performance

[rvignolo]

Is there any way we can improve the performance when using a Python function for an ODEProblem?

I know that I can use something like:

jul_f = de.eval("""
function f(du,u,p,t)
  x, y, z = u
  sigma, rho, beta = p
  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.)
p = [10.0,28.0,2.66]
prob = de.ODEProblem(jul_f, u0, tspan, p)
sol = de.solve(prob,saveat=t,abstol=1e-8,reltol=1e-8)

instead of

def f(du,u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    du[0] = sigma * (y - x)
    du[1] = x * (rho - z) - y
    du[2] = x * y - beta * z

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,2.66]
prob = de.ODEProblem(f, u0, tspan, p)
sol = de.solve(prob,saveat=t,abstol=1e-8,reltol=1e-8)

but I cannot translate my function to Julia code (it is way more complicated than this example and it relies in many Python objects and auxiliary functions).
Maybe what is making the code run slow is that PyCall.jl needs to convert python objects to julia objects and vice versa. So... is there any way I can let Julia or Python know the types of du, u and t.

Thanks! If there is no easy way to hack this it is fine, but I just wanted to know if I could make my code run faster.

[tkf]

If you have a sympy expression, it should be easy to compile a Julia function from it. Here is an example that is useful for low-dimensional system #32. This example also contains a way to "trace" simple Python functions like you quoted using sympy (see if __name__ == "__main__": of #32).

[rvignolo]

Hi, thanks for your reply!

Yeah, I saw that issue some days ago and it is great for simple examples. However, for more complicated cases it is not. For example, imagine your f function depends on the evaluation at time t of interpolators that are instantiated in Python (such as SciPy's interp1d or interp2d objects). In this case, your jl_lambdify won't work, because you will end with a Julia function that depend on object that does not exist.

On the other hand, if the f expression depend on a functional such as a Derivative I would have to compute it first on Python and then migrate the resulting expression, but if the argument of the Derivative has interpolated data I would have to use a finite difference scheme with a predefined step h.

Lastly and more important, if the f expression depend on an Integral of interpolated data, I won't be able to solve it.

Any thoughts?

[tkf]

If your function calls complex functions in SciPy, wouldn't they be the dominating factor in the performance?

If your problem size is somewhat intermediate so that using sympy does not make sense while Python/Julia interface still is not negligible, maybe you can warp your Python function with something like this (untested):

fjl = Main.eval("f -> (args...) -> pycall(f, Nothing, args...)")(f)

This would eliminate the dynamic conversion of the returned value (which is not used as you are using the in-place version).

[rvignolo]

I see, thanks for the tip! I will work on that and show you some performance results.

On the other hand, is it necessary (and, in that case, possible) to specify the type of args? I believe that it is not necessary because when f is evaluated from Julia, we already know the arguments types and dynamic conversion is not needed. But I want to be sure.

Thanks!

[tkf]

Yeah, in Julia, argument types are specialized automatically. So, you don't have to do anything for performance. Argument types are specified mainly for dispatching.

[ChrisRackauckas]

If you work out a nice example, it would be worth adding it to the docs.

[ChrisRackauckas]

Use the new fast_prob = de.jit(prob)