Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Switch backend from PyCall to PythonCall and improve package management #118

Merged
merged 15 commits into from
Oct 10, 2023

Conversation

LilithHafner
Copy link
Member

@LilithHafner LilithHafner commented Sep 19, 2023

This PR

  • Changes the backend from PyCall to PythonCall
  • Removes most of the glue code that makes Py[thon]Call and DifferentialEquaitons compatible (that code has been upstreamed into package extensions Add PyCall extension SciMLBase.jl#502 & Add PythonCall Extension SciMLBase.jl#519)
  • Makes the install script perform a single Pkg operation instead of several to reduce noise and the risk of re-precompilation
  • Stops using the user's default environment (uses @diffeqpy instead)
  • Automatically installs Julia & packages if Julia is missing when loading de or ode (previously this only happened when loading de)
  • Updates all examples in the readme so that they still work.

This PR is breaking because

  • from julia import Main no longer works. This can be replaced with from juliacall import Main, but julia.Main (from PyCall.jl) and juliacall.Main (from PythonCall.jl) have slightly different semantics around conversion and eval/seval. For the examples in the README, all that is needed is to replace
from julia import Main
Main.eval("1+1")

with

de.seval("1+1")
  • Automatic conversion semantics are different. For calls into the SciML ecosystem, this is mostly patched, but if folks are heavily utilizing passing data back and forth between Python and Julia in their own code and do not have a direct dependency on julia/PyCall.jl, this will break their code.
  • Differential equation results are not automatically converted when passed back to Python. Specifically, sol.u is often Vector{Vector{T}} instead of a numpy array. This necessitates a bit of reshaping before passing to matplotlib.pyplot. de.stack is available for converting that Vector{Vector{T}} into a Matrix{T} which can be plotted, though even then it is transposed compared to what diffeqpy previously returned.
  • It breaks some numba interop. (TODO)

See changes to the README to see the impact of the breaking changes on downstream code.

@LilithHafner LilithHafner marked this pull request as draft September 19, 2023 17:31
README.md Outdated
Comment on lines 102 to 107
return -u
def f(du,u,p,t):
du[0] = -u[0]

u0 = 0.5
u0 = np.array([0.5])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason why scalars are no longer supported?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With PythonCall, a python float is passed to Julia as an Py with a pytype of float. This object works in many places where one would expect a float, but not all. For example, it's not a Number and is not iterable (iteration falls back to python's iteration, and float is not iterable in python) so the check any(isnan, x) will error if we pass in a python float.

It's possible to use scalars, but they must be converted when they are passed into a differential equations problem (not too hard) and also whenever they are returned from a function that computes derivatives (somewhat harder).

Similarly, any use of non-mutating derivative functions requires the passed object to be converted from a python object to a Julia object, and the PythonCall semantics aren't quite what I'd like for that. For example, if a python function returns x = [1,2,3] we can call pyconvert(Any, x) which will convert it to a PyList{Any}, and the eltype of Any will cause errors down the line.

These conversion issues can be overcome, but I worry they would be brittle and could more than double the size of the diffeqpy codebase.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it would be too bad to add a PythonCall extension where at the DEProblem constructors and such we can wrap the f to "do the right thing" and convert the u0/p types.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, I'll give it a shot.

diffeqpy/setup.jl Outdated Show resolved Hide resolved
@LilithHafner
Copy link
Member Author

CI fails because of #379. There might be more errors left to fix even after that, though.

@LilithHafner LilithHafner changed the title Switch backend from PyCall to PythonCall Switch backend from PyCall to PythonCall and improve package management Oct 10, 2023
@ChrisRackauckas ChrisRackauckas marked this pull request as ready for review October 10, 2023 14:30
@ChrisRackauckas
Copy link
Member

Oops I accidentally hit ready to review 😅

@LilithHafner
Copy link
Member Author

Aside from this breaking some interactions with numba, this is ready.

@ChrisRackauckas
Copy link
Member

Does this interact well with ModelingToolkit? Try using the modelingtoolkitize function for translating it to a symbolic MTK model and building the function in Julia, i.e. https://docs.sciml.ai/ModelingToolkit/stable/tutorials/modelingtoolkitize/

@LilithHafner
Copy link
Member Author

It doesn't come bundled with ModelningToolkit, so users who want to use modelingtoolkitize would need to install that package separately, but aside from that, yes:

>>> from diffeqpy import de
>>> rober = de.seval("""function rober(du, u, p, t)
...     y₁, y₂, y₃ = u
...     k₁, k₂, k₃ = p
...     du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
...     du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
...     du[3] = k₂ * y₂^2
...     nothing
... end""")
>>> 
>>> prob = de.ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))
>>> sys = de.modelingtoolkitize(prob)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/x/.julia/packages/PythonCall/qTEA1/src/jlwrap/any.jl", line 195, in __getattr__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_getattr)), k)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: Julia: UndefVarError: `modelingtoolkitize` not defined
>>> from juliacall import Main
>>> Main.seval("import Pkg; Pkg.add(\"ModelingToolkit\")")
[...]
>>> Main.seval("using ModelingToolkit")
[ Info: Precompiling ModelingToolkit [961ee093-0014-501f-94e3-6117800e7a78]
>>> mt = Main.seval("ModelingToolkit")
>>> sys = mt.modelingtoolkitize(prob)
>>> sys
Julia:
Model ##MTKizedODE#292 with 3 equations
States (3):
  x₁(t) [defaults to 1.0]
  x₂(t) [defaults to 0.0]
  x₃(t) [defaults to 0.0]
Parameters (3):
  α₁ [defaults to 0.04]
  α₂ [defaults to 3.0e7]
  α₃ [defaults to 10000.0]

@ChrisRackauckas
Copy link
Member

Awesome. Can you ship ModelingToolkit with it? We do that with diffeqr because it's a big acceleration. See: https://github.com/SciML/diffeqr#performance-enhancements . I can update the docs for that, removing numba and replacing it with ModelingToolkit. Just double check that the simulation with ModelingToolkit is unimpeded (i.e. it's similar in time to the optimal version in just Julia) and I think this is good to merge!

@LilithHafner
Copy link
Member Author

I can confirm that modelingtoolkitize can accelerate an ODEProblem defined with a Python function to within a factor of 2 of pure Julia speeds.

Methodology

import time
from diffeqpy import de
from juliacall import Main
mtk = Main.seval('import Pkg; Pkg.activate("mtk", shared=true); using ModelingToolkit; ModelingToolkit')
def rober(u,p,t):
    du0 = p[0]*(u[1]-u[0])
    du1 = u[0]*(p[1]-u[2]) - u[1]
    du2 = u[0]*u[1] - p[2]*u[2]
    return [du0, du1, du2]

rober_jl = de.seval("""
function rober_jl(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""")

prob = de.ODEProblem(rober, [1.0,0.0,0.0], (0.0,1e5), (0.04, 3e7, 1e4))
sys_mtk = mtk.modelingtoolkitize(prob)
prob_mtk = de.ODEProblem(sys_mtk)
prob_jl = de.ODEProblem(rober_jl, [1.0,0.0,0.0], (0.0,1e5), (0.04, 3e7, 1e4))

def timeit(f, str, n):
    start = time.time()
    for i in range(n):
        f()
    print(str, "mean of", n, "trials", (time.time() - start)/n)

timeit(lambda: de.solve(prob), "python", 1)
timeit(lambda: de.solve(prob), "python", 100)
timeit(lambda: de.solve(prob_mtk), "python mtk", 1)
timeit(lambda: de.solve(prob_mtk), "python mtk", 100)
timeit(lambda: de.solve(prob_jl), "julia", 1)
timeit(lambda: de.solve(prob_jl), "julia", 100)

Results

python mean of 1 trials 2.3302979469299316
python mean of 100 trials 0.013549079895019531
python mtk mean of 1 trials 1.6208293437957764
python mtk mean of 100 trials 0.0003797006607055664
julia mean of 1 trials 0.9586193561553955
julia mean of 100 trials 0.00020106077194213866

I'm thinking of adding a function like diffeqr's jtioptimize_ode called de.jit that takes in any type of problem, calls modelingtoolkitize on it and then converts it back to its original problem type. ode won't get ModelingToolkit or the jit because it's whole point is lowered compilation times.

@ChrisRackauckas
Copy link
Member

I'm thinking of adding a function like diffeqr's jtioptimize_ode called de.jit that takes in any type of problem, calls modelingtoolkitize on it and then converts it back to its original problem type. ode won't get ModelingToolkit or the jit because it's whole point is lowered compilation times.

Yes that would be awesome! Let's do that as a follow up.

It would be cool to try and see if all of the GPU ensemble stuff works as well:

https://github.com/SciML/diffeqr#gpu-accelerated-ode-solving-of-ensembles

I assume that there isn't anything that's required to make it work? But's lets follow this up as new PRs.

@ChrisRackauckas ChrisRackauckas merged commit 6770470 into SciML:master Oct 10, 2023
2 checks passed
@LilithHafner LilithHafner deleted the pythoncall branch October 10, 2023 20:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants