Switch backend from PyCall to PythonCall and improve package management

[LilithHafner]

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.

[ChrisRackauckas]

Is there a reason why scalars are no longer supported?

[LilithHafner]

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.

[ChrisRackauckas]

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.

[LilithHafner]

Okay, I'll give it a shot.

[LilithHafner]

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

[ChrisRackauckas]

Oops I accidentally hit ready to review 😅

[LilithHafner]

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

[ChrisRackauckas]

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]

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]

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]

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]

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.