Python vs. Julia ODE Solver

General Usage Performance
22

It is. But you are not running exactly the “my” version, because I have much less allocations than you.

Just to be clear, it was not possible to do this in the first version because of the rand(Float64,(N,1)) definition I mentioned above:

julia> v = rand(Float64,(5,1));

julia> y = rand(5);

julia> z = zeros(5);

julia> z[1:5] .= v .+ y
ERROR: DimensionMismatch("cannot broadcast array to have fewer dimensions")

The problem is that v here is a two-dimension array (even though the second dimension is 1):

julia> v
5×1 Array{Float64,2}:
 0.23460063971299383
 0.977599317555971
 0.11413962191900962
 0.9259610441668678
 0.0253857101609829

That is solved if v is a vector:

julia> v = rand(Float64,5);

julia> z[1:5] .= v .+ y
5-element view(::Array{Float64,1}, 1:5) with eltype Float64:
 1.7707818802990465
 1.5662031379445822
 1.1979284388950557
 0.7127751475367243
 1.1919315206952428

It might be a parallel discussion why/if the first broadcasting should work.

The final code I am running is here, It seems to be faster than what you are running there still.

Code 
using OrdinaryDiffEq, BenchmarkTools, LSODA, Plots, Profile
using LinearAlgebra

function forward_RHS(du, u, p, t)
    v_s = view(u, 1:p[2])
    u_s = view(u, p[2]+1:2*p[2])
    w_s = view(u, 2*p[2]+1:3*p[2])
    s = view(u, 3*p[2]+1:4*p[2])

    p[26] .= t .* view(p[23],:,1)
    p[26] .= sin.(p[26])
    mul!(p[24],p[20],p[26])

    du[1:p[2]] .= ((p[3].*(v_s .- p[4]).^2 .+ p[5]).*(v_s .- 1.0)
                  .- p[12].*u_s.*(v_s .- p[14]) .- p[13].*w_s.*(v_s .- p[15])
                  .+ p[22]
                  .+ p[24]
                  .+ mul!(p[25],p[21],s))./p[16]
    du[p[2]+1:2*p[2]] .= (exp.(p[6].*(v_s .- p[7])) .+ p[8] .- u_s)./p[17]
    du[2*p[2]+1:3*p[2]] .= (p[9].*(v_s .- p[10]).^2 .+ p[11] .- w_s)./p[18]
    du[3*p[2]+1:4*p[2]] .= (-1 .* s .+ (v_s .> 0.25).*(v_s .< 0.3).*(view(du,1:p[2]) .> 0.0).*view(du,1:p[2])./0.05)./p[19]

    nothing

end

function runprob(N)

 parameters = (2, N,
              -500, 0.2, 0.0, 100.0, 0.7, 0.0, 15.0, 0.24, 0.02, 1700.0, 8.0, 0.19, -0.5,
              25.0, 6.0, 12.0, 10.0,
              2.0*1.4 .* (rand(Float64,N,2) .- 0.5), #20 
              0.5 .* (rand(Float64,N,N) .- 0.8), #21
              0.14 .* (rand(Float64,N) .- 0.5), #22
              rand(Float64,2,4),  #23
              zeros(N), #24
              zeros(N), #25 
              zeros(2)) #26

  prob = ODEProblem(forward_RHS, zeros(4*N), (0.0, 500.0), parameters)

  solve(prob, lsoda(), saveat=0.1, reltol=1e-8, abstol=1e-8)

end

#@time solution = runprob(N)

I get:

julia> @time runprob(20);
  0.132728 seconds (195.93 k allocations: 18.890 MiB)

julia> @time runprob(90);
  1.365390 seconds (512.74 k allocations: 54.033 MiB)

julia> @time runprob(100);
  5.902604 seconds (563.84 k allocations: 59.490 MiB, 0.33% gc time)

julia> @time runprob(120);
  6.846018 seconds (559.88 k allocations: 62.312 MiB, 0.26% gc time)

julia> @time runprob(200);
 16.374404 seconds (837.08 k allocations: 96.083 MiB, 0.10% gc time)

Perhaps there are indeed some “jumps” in the performance?

3 Likes