pyjulia slower than expected for some operations

[mjohnson541]

I've been using pyjulia to wrap a julia package, but the wrapper seems slower than I would expect. The pyjulia code using diffeqpy involves some setup to generate the react object which has a julia ODEProblem with a julia function attached, but the main operation is a solve like this:

sol = de.solve(react.ode,de.CVODE_BDF(),abstol=1e-16,reltol=1e-6)

In pure julia at first evaluation this takes about 20 seconds, and after that it takes about 6 seconds
In pyjulia at first evaluation this takes about 30 seconds and after that it takes about 18 seconds.

I get the same result if I interface directly with DifferentialEquations without diffeqpy.

I can't think of any reason it should be that different. Does some of the jit compilation have to happen every evaluation in pyjulia? I'm using the python-jl fix if that affects things.

[tkf]

If you use PyJulia for passing around Julia functions and values, it should be as fast as pure Julia solution (if the main computation is done in Julia). My guess is that there are some quirks in how PyCall wrap/unwrap things (or usage of it) and the ODE solver accidentally calls Julia functions via Python. I'd run Julia profiler to see if/where PyCall shows up in the hot loop.

[mjohnson541]

Thanks! That was a very convoluted debug XD. It turned out to not be a performance issue. It ultimately stemmed from the fact that typeof(obj).name in pyjulia returns as Module.type while julia proper it returns as just type. This was causing it to think it couldn't do certain optimizations. Probably doesn't affect many people.

[tkf]

Are you defining the derivative for the ODE as a Python function? What is the dimension of the system? If it's low dimensional and the derivative is defined as a Python function, there is not much we can do to improve the performance.

For diffeqpy specific questions, it's best to ask @ChrisRackauckas

[ChrisRackauckas]

Yeah this is something I noticed way back in this blog post (https://www.stochasticlifestyle.com/why-numba-and-cython-are-not-substitutes-for-julia/) but never really figured out how to get around it. I think the issue is that even if we eval a Julia function to define the derivative, the resulting function is a Python function that points to a Julia function, and so each call has that indirection in there. Somehow we want to wrap that so that way in the problem constructors or from the Julia side it unwraps itself and directly calls the function in there. I've wondered if there's a way to do like Main.function_eval(...) that would then tag it to know that it's a function and skip that step in higher order functions.

[mjohnson541]

No, so it was defined as a Julia function and I've now achieved performance pairity between the python wrapper and the julia code.

What actually was happening is that I essentially had a check that looks to see if I have a specialized vectorized object for an arbitrary type, so if my object is some type Arrhenius I essentially check:

Symbol(string(typeof(obj).name) * "vec") in definedfcnlist

Where definedfcnlist contains :Arrheniusvec. In julia proper string(typeof(obj).name) == Arrhenius, but in pyjulia string(typeof(obj).name) == ReactionMechanismSimulator.Arrhenius this meant in pyjulia it would search for :ReactionMechanismSimulator.Arrheniusvec instead, not find it and then not vectorize that calculation slowing it down. Once I figured out the output of typeof(obj).name was different in pyjulia and julia It was easy enough to resolve.

[tkf]

If the ODE function is implemented in Julia, there is, in principle, no need to call back in Python. So, the performance difference is still a mystery to me. It might be some quirks in PyCall.jl or some bad interactions with how DifferentialEquations.jl works. I need to debug this myself to find out exactly what's going on (though I don't have enough bandwidth to handle this right now). If you can share an MWE that would be great.

[ChrisRackauckas]

It looks like the difference has really dropped. It could've been a WSL problem before.

from julia.api import Julia
jl = Julia(compiled_modules=False)
from diffeqpy import de
from julia import 
import timeit

import numpy as np
jul_f = Main.eval("""
function f(du,u,p,t)
  x, y, z = u[1], u[2], u[3]
  sigma, rho, beta = p[1], p[2], p[3]
  du[1] = sigma * (y - x)
  du[2] = x * (rho - z) - y
  du[3] = x * y - beta * z
end""")
u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
t = np.linspace(0, 100, 1001)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(jul_f, u0, tspan, p)
sol = de.solve(prob,saveat=t,abstol=1e-8,reltol=1e-8)
 
def time_func():
    sol = de.solve(prob,saveat=t,abstol=1e-8,reltol=1e-8)
 
time_func()
timeit.Timer(time_func).timeit(number=10000)/10000 # 0.003481317869999998 seconds

vs

using DifferentialEquations, BenchmarkTools
function lorenz!(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end
u0 = [1.0;0.0;0.0]
p = [10.0,28.0,8/3]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)
@btime solve(prob,Tsit5(),saveat=0.1,reltol=1e-8,abstol=1e-8) # 3.283 ms (1055 allocations: 154.44 KiB)

That last amount could just be simple setup stuff, so this can be closed.