From 1d1ab4bf95d3604b7c1871717c34aea3e87fe0e5 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 12 Sep 2023 14:27:09 -0500 Subject: [PATCH 01/14] WME --- diffeqpy/de.py | 8 +++++-- diffeqpy/install.jl | 6 ++--- diffeqpy/ode.py | 8 +++++-- diffeqpy/setup.jl | 46 ++++++++++++++++++++++++++++++++------ diffeqpy/tests/test_dde.py | 2 +- setup.py | 2 +- 6 files changed, 56 insertions(+), 16 deletions(-) diff --git a/diffeqpy/de.py b/diffeqpy/de.py index a0567ce..bd868dc 100644 --- a/diffeqpy/de.py +++ b/diffeqpy/de.py @@ -7,10 +7,14 @@ _ensure_installed() # PyJulia have to be loaded after `_ensure_installed()` -from julia import Main +from juliacall import Main script_dir = os.path.dirname(os.path.realpath(__file__)) Main.include(os.path.join(script_dir, "setup.jl")) -from julia import DifferentialEquations +# TODO: find a better way to do this or upstream this function into juliacall +def load_julia_package(name): + return Main.seval(f"using {name}; {name}") + +DifferentialEquations = load_julia_package("DifferentialEquations") sys.modules[__name__] = DifferentialEquations # mutate myself diff --git a/diffeqpy/install.jl b/diffeqpy/install.jl index 00dfa35..6157c9a 100644 --- a/diffeqpy/install.jl +++ b/diffeqpy/install.jl @@ -2,7 +2,7 @@ using Pkg Pkg.add("DifferentialEquations") Pkg.add("OrdinaryDiffEq") Pkg.add("DiffEqBase") -Pkg.add("PyCall") -Pkg.build("PyCall") +Pkg.develop("PythonCall") +Pkg.build("PythonCall") using DifferentialEquations -using PyCall +using PythonCall diff --git a/diffeqpy/ode.py b/diffeqpy/ode.py index e0adb50..38b17f2 100644 --- a/diffeqpy/ode.py +++ b/diffeqpy/ode.py @@ -1,10 +1,14 @@ import os import sys -from julia import Main +from juliacall import Main script_dir = os.path.dirname(os.path.realpath(__file__)) Main.include(os.path.join(script_dir, "setup.jl")) -from julia import OrdinaryDiffEq +# TODO: find a better way to do this or upstream this function into juliacall +def load_julia_package(name): + return Main.seval(f"using {name}; {name}") + +OrdinaryDiffEq = load_julia_package("OrdinaryDiffEq") sys.modules[__name__] = OrdinaryDiffEq # mutate myself diff --git a/diffeqpy/setup.jl b/diffeqpy/setup.jl index e58cde6..055f1f1 100644 --- a/diffeqpy/setup.jl +++ b/diffeqpy/setup.jl @@ -6,16 +6,48 @@ catch err rethrow() end -import PyCall +import PythonCall -PyCall.PyObject(::typeof(DiffEqBase.solve)) = - PyCall.pyfunctionret(DiffEqBase.solve,Any,Vararg{PyCall.PyAny}) +# PythonCall does not implicitly convert the return values of Python functions to Julia +# values. If a user writes a function in Python and passes it to a solver, that function +# will return objects of type Py, but the solver will expect Julia objects. To solve this, +# we add explicit conversion any time a Py is passed into the DifferentialEquations.jl +# ecosystem as a function. -function DiffEqBase.numargs(f::PyCall.PyObject) - inspect = PyCall.pyimport("inspect") - PyCall.hasproperty(f,:py_func) ? _f = f.py_func : _f = f +# TODO: upstream this into a package extension between PythonCall and DiffEqBase +wrap(f::PythonCall.Py) = (args...; kws...) -> pyconvert(Any, f(args...; kws...)) +using DiffEqBase +DiffEqBase.DAEFunction(f::Py) = DiffEqBase.DAEFunction(wrap(f)) +DiffEqBase.ODEFunction{T, U}(f::Py) where {T, U} = DiffEqBase.ODEFunction{T, U}(wrap(f)) +# ... +# TODO: make this cover all entrypoints or use some other mechanism to ensure we get every +# function + +# function wrap(f::PythonCall.Py) +# # preserve the number of function arguments because DifferentialEquations uses that to +# # determine if the function is mutating +# inspect = PythonCall.pyimport("inspect") +# PythonCall.hasproperty(f,:py_func) ? _f = f.py_func : _f = f +# nargs = length(first(inspect.getfullargspec(_f))) +# if PythonCall.pyconvert(Bool, inspect.ismethod(_f)) +# # `f` is a bound method (i.e., `self.f`) but `getfullargspec` +# # includes `self` in the `args` list. So, we are subtracting +# # 1 manually here: +# nargs -= 1 +# end +# @eval (args::Vararg{Any, $nargs}; kw...) -> pyconvert(Any, $f(args...; kw...)) +# end + + +# PyCall.PyObject(::typeof(DiffEqBase.solve)) = + # PyCall.pyfunctionret(DiffEqBase.solve,Any,Vararg{PyCall.PyAny}) + +# TODO: upstream this into a package extension between PythonCall and DiffEqBase +function DiffEqBase.numargs(f::PythonCall.Py) + inspect = PythonCall.pyimport("inspect") + PythonCall.hasproperty(f,:py_func) ? _f = f.py_func : _f = f nargs = length(first(inspect.getfullargspec(_f))) - if inspect.ismethod(_f) + if PythonCall.pyconvert(Bool, inspect.ismethod(_f)) # `f` is a bound method (i.e., `self.f`) but `getfullargspec` # includes `self` in the `args` list. So, we are subtracting # 1 manually here: diff --git a/diffeqpy/tests/test_dde.py b/diffeqpy/tests/test_dde.py index 0b4b1ea..315e400 100644 --- a/diffeqpy/tests/test_dde.py +++ b/diffeqpy/tests/test_dde.py @@ -1,4 +1,4 @@ -from julia import Main +from juliacall import Main from .. import de diff --git a/setup.py b/setup.py index e72f068..807d8c8 100644 --- a/setup.py +++ b/setup.py @@ -24,6 +24,6 @@ def readme(): author_email='contact@juliadiffeq.org', license='MIT', packages=['diffeqpy','diffeqpy.tests'], - install_requires=['julia>=0.2', 'jill'], + install_requires=['juliacall>=0.9.14', 'jill'], # TODO: specify looser compat bound for juliacall if safe. include_package_data=True, zip_safe=False) From 6767b16937c6662a0cfc102c5bad9345e4318b61 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 12 Sep 2023 14:36:24 -0500 Subject: [PATCH 02/14] update README.md to reflect BREAKING CHANGES --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 7053153..35a492b 100644 --- a/README.md +++ b/README.md @@ -164,12 +164,12 @@ sol = de.solve(prob) ``` Additionally, you can directly define the functions in Julia. This will allow -for more specialization and could be helpful to increase the efficiency over -the Numba version for repeat or long calls. This is done via `julia.Main.eval`: +for more specialization and could be helpful to increase the efficiency over the +Numba version for repeat or long calls. This is done via `juliacall.Main.seval`: ```py -from julia import Main -jul_f = Main.eval("(u,p,t)->-u") # Define the anonymous function in Julia +from juliacall import Main +jul_f = Main.seval("(u,p,t)->-u") # Define the anonymous function in Julia prob = de.ODEProblem(jul_f, u0, tspan) sol = de.solve(prob) ``` From fe6ec4fc019882987df929c3efcdc8b026487a7d Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 19 Sep 2023 12:23:42 -0500 Subject: [PATCH 03/14] clean up and update examples in README to work --- README.md | 37 +++++++++++-------- diffeqpy/setup.jl | 93 ++++++++++++++++++++++++++++++----------------- 2 files changed, 81 insertions(+), 49 deletions(-) diff --git a/README.md b/README.md index 35a492b..7e9dfff 100644 --- a/README.md +++ b/README.md @@ -75,6 +75,8 @@ translate these docs to Python code. Python does not allow `!` in function names, so this is also [a limitation of pyjulia](https://pyjulia.readthedocs.io/en/latest/limitations.html#mismatch-in-valid-set-of-identifiers) To use functions which on the Julia side have a `!`, like `step!`, replace `!` by `_b`, for example: +TODO: + ```julia from diffeqpy import de def f(u,p,t): @@ -95,11 +97,12 @@ is valid Python code for using [the integrator interface](https://diffeq.sciml.a ```py from diffeqpy import de +import numpy as np -def f(u,p,t): - return -u +def f(du,u,p,t): + du[0] = -u[0] -u0 = 0.5 +u0 = np.array([0.5]) tspan = (0., 1.) prob = de.ODEProblem(f, u0, tspan) sol = de.solve(prob) @@ -116,7 +119,7 @@ We can plot the solution values using matplotlib: ```py import matplotlib.pyplot as plt -plt.plot(sol.t,sol.u) +plt.plot(sol.t,[u[0] for u in sol.u]) plt.show() ``` @@ -128,7 +131,7 @@ We can utilize the interpolation to get a finer solution: import numpy t = numpy.linspace(0,1,100) u = sol(t) -plt.plot(t,u) +plt.plot(t,[ui[0] for ui in u]) plt.show() ``` @@ -165,12 +168,11 @@ sol = de.solve(prob) Additionally, you can directly define the functions in Julia. This will allow for more specialization and could be helpful to increase the efficiency over the -Numba version for repeat or long calls. This is done via `juliacall.Main.seval`: +Numba version for repeat or long calls. This is done via `de.seval` or `ode.seval`: ```py -from juliacall import Main -jul_f = Main.seval("(u,p,t)->-u") # Define the anonymous function in Julia -prob = de.ODEProblem(jul_f, u0, tspan) +jul_f = de.seval("(u,p,t)->-u") # Define the anonymous function in Julia +prob = de.ODEProblem(jul_f, u0[0], tspan) sol = de.solve(prob) ``` @@ -182,18 +184,21 @@ To solve systems of ODEs, simply use an array as your initial condition and define `f` as an array function: ```py -def f(u,p,t): +def f(du,u,p,t): x, y, z = u sigma, rho, beta = p - return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z] + du[0] = sigma * (y - x) + du[1] = x * (rho - z) - y + du[2] = x * y - beta * z -u0 = [1.0,0.0,0.0] +u0 = np.array([1.0,0.0,0.0]) tspan = (0., 100.) -p = [10.0,28.0,8/3] +p = (10.0,28.0,8/3) prob = de.ODEProblem(f, u0, tspan, p) sol = de.solve(prob,saveat=0.01) +u = np.array([list(ui) for ui in sol.u]) -plt.plot(sol.t,sol.u) +plt.plot(sol.t,u) plt.show() ``` @@ -202,7 +207,7 @@ plt.show() or we can draw the phase plot: ```py -ut = numpy.transpose(sol.u) +ut = numpy.transpose(u) from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') @@ -214,6 +219,8 @@ plt.show() ### In-Place Mutating Form +TODO + When dealing with systems of equations, in many cases it's helpful to reduce memory allocations by using mutating functions. In diffeqpy, the mutating form adds the mutating vector to the front. Let's make a fast version of the diff --git a/diffeqpy/setup.jl b/diffeqpy/setup.jl index 055f1f1..a25d8a2 100644 --- a/diffeqpy/setup.jl +++ b/diffeqpy/setup.jl @@ -8,40 +8,6 @@ end import PythonCall -# PythonCall does not implicitly convert the return values of Python functions to Julia -# values. If a user writes a function in Python and passes it to a solver, that function -# will return objects of type Py, but the solver will expect Julia objects. To solve this, -# we add explicit conversion any time a Py is passed into the DifferentialEquations.jl -# ecosystem as a function. - -# TODO: upstream this into a package extension between PythonCall and DiffEqBase -wrap(f::PythonCall.Py) = (args...; kws...) -> pyconvert(Any, f(args...; kws...)) -using DiffEqBase -DiffEqBase.DAEFunction(f::Py) = DiffEqBase.DAEFunction(wrap(f)) -DiffEqBase.ODEFunction{T, U}(f::Py) where {T, U} = DiffEqBase.ODEFunction{T, U}(wrap(f)) -# ... -# TODO: make this cover all entrypoints or use some other mechanism to ensure we get every -# function - -# function wrap(f::PythonCall.Py) -# # preserve the number of function arguments because DifferentialEquations uses that to -# # determine if the function is mutating -# inspect = PythonCall.pyimport("inspect") -# PythonCall.hasproperty(f,:py_func) ? _f = f.py_func : _f = f -# nargs = length(first(inspect.getfullargspec(_f))) -# if PythonCall.pyconvert(Bool, inspect.ismethod(_f)) -# # `f` is a bound method (i.e., `self.f`) but `getfullargspec` -# # includes `self` in the `args` list. So, we are subtracting -# # 1 manually here: -# nargs -= 1 -# end -# @eval (args::Vararg{Any, $nargs}; kw...) -> pyconvert(Any, $f(args...; kw...)) -# end - - -# PyCall.PyObject(::typeof(DiffEqBase.solve)) = - # PyCall.pyfunctionret(DiffEqBase.solve,Any,Vararg{PyCall.PyAny}) - # TODO: upstream this into a package extension between PythonCall and DiffEqBase function DiffEqBase.numargs(f::PythonCall.Py) inspect = PythonCall.pyimport("inspect") @@ -64,3 +30,62 @@ catch err @error "Failed to import DifferentialEquations.jl" exception = (err, catch_backtrace()) rethrow() end + +#= + +from diffeqpy import de +import numpy as np + +def f(du,u,p,t): + du[0] = -u[0] + +u0 = np.array([0.5]) +tspan = (0., 1.) +prob = de.ODEProblem(f, u0, tspan) +sol = de.solve(prob) + + + +import matplotlib.pyplot as plt +plt.plot(sol.t,[u[0] for u in sol.u]) +plt.show() + + + +import numpy +t = numpy.linspace(0,1,100) +u = sol(t) +plt.plot(t,[ui[0] for ui in u]) +plt.show() + + +jul_f = de.seval("(u,p,t)->-u") # Define the anonymous function in Julia +prob = de.ODEProblem(jul_f, u0[0], tspan) +sol = de.solve(prob) + + + +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 = np.array([1.0,0.0,0.0]) +tspan = (0., 100.) +p = (10.0,28.0,8/3) +prob = de.ODEProblem(f, u0, tspan, p) +sol = de.solve(prob,saveat=0.01) +u = np.array([list(ui) for ui in sol.u]) + +plt.plot(sol.t,u) +plt.show() + +ut = numpy.transpose(u) +from mpl_toolkits.mplot3d import Axes3D +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.plot(ut[0,:],ut[1,:],ut[2,:]) +plt.show() +=# \ No newline at end of file From bd4facc2f946b953e8fe259c164234857e057800 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 19 Sep 2023 12:24:44 -0500 Subject: [PATCH 04/14] cleanup --- diffeqpy/setup.jl | 59 ----------------------------------------------- 1 file changed, 59 deletions(-) diff --git a/diffeqpy/setup.jl b/diffeqpy/setup.jl index a25d8a2..1d92213 100644 --- a/diffeqpy/setup.jl +++ b/diffeqpy/setup.jl @@ -30,62 +30,3 @@ catch err @error "Failed to import DifferentialEquations.jl" exception = (err, catch_backtrace()) rethrow() end - -#= - -from diffeqpy import de -import numpy as np - -def f(du,u,p,t): - du[0] = -u[0] - -u0 = np.array([0.5]) -tspan = (0., 1.) -prob = de.ODEProblem(f, u0, tspan) -sol = de.solve(prob) - - - -import matplotlib.pyplot as plt -plt.plot(sol.t,[u[0] for u in sol.u]) -plt.show() - - - -import numpy -t = numpy.linspace(0,1,100) -u = sol(t) -plt.plot(t,[ui[0] for ui in u]) -plt.show() - - -jul_f = de.seval("(u,p,t)->-u") # Define the anonymous function in Julia -prob = de.ODEProblem(jul_f, u0[0], tspan) -sol = de.solve(prob) - - - -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 = np.array([1.0,0.0,0.0]) -tspan = (0., 100.) -p = (10.0,28.0,8/3) -prob = de.ODEProblem(f, u0, tspan, p) -sol = de.solve(prob,saveat=0.01) -u = np.array([list(ui) for ui in sol.u]) - -plt.plot(sol.t,u) -plt.show() - -ut = numpy.transpose(u) -from mpl_toolkits.mplot3d import Axes3D -fig = plt.figure() -ax = fig.add_subplot(111, projection='3d') -ax.plot(ut[0,:],ut[1,:],ut[2,:]) -plt.show() -=# \ No newline at end of file From 48a1a7f7fc1f6a4efb42629ac88268fe8a197318 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Mon, 9 Oct 2023 18:40:12 -0500 Subject: [PATCH 05/14] remove setup.jl --- diffeqpy/de.py | 5 +---- diffeqpy/install.jl | 3 +-- diffeqpy/ode.py | 5 +---- diffeqpy/setup.jl | 32 -------------------------------- 4 files changed, 3 insertions(+), 42 deletions(-) delete mode 100644 diffeqpy/setup.jl diff --git a/diffeqpy/de.py b/diffeqpy/de.py index bd868dc..4b666d0 100644 --- a/diffeqpy/de.py +++ b/diffeqpy/de.py @@ -9,12 +9,9 @@ # PyJulia have to be loaded after `_ensure_installed()` from juliacall import Main -script_dir = os.path.dirname(os.path.realpath(__file__)) -Main.include(os.path.join(script_dir, "setup.jl")) - # TODO: find a better way to do this or upstream this function into juliacall def load_julia_package(name): - return Main.seval(f"using {name}; {name}") + return Main.seval(f"using {name}: {name}; {name}") DifferentialEquations = load_julia_package("DifferentialEquations") sys.modules[__name__] = DifferentialEquations # mutate myself diff --git a/diffeqpy/install.jl b/diffeqpy/install.jl index 6157c9a..bb02417 100644 --- a/diffeqpy/install.jl +++ b/diffeqpy/install.jl @@ -1,8 +1,7 @@ using Pkg Pkg.add("DifferentialEquations") Pkg.add("OrdinaryDiffEq") -Pkg.add("DiffEqBase") -Pkg.develop("PythonCall") +Pkg.add("PythonCall") Pkg.build("PythonCall") using DifferentialEquations using PythonCall diff --git a/diffeqpy/ode.py b/diffeqpy/ode.py index 38b17f2..c3ab6e7 100644 --- a/diffeqpy/ode.py +++ b/diffeqpy/ode.py @@ -3,12 +3,9 @@ from juliacall import Main -script_dir = os.path.dirname(os.path.realpath(__file__)) -Main.include(os.path.join(script_dir, "setup.jl")) - # TODO: find a better way to do this or upstream this function into juliacall def load_julia_package(name): - return Main.seval(f"using {name}; {name}") + return Main.seval(f"using {name}: {name}; {name}") OrdinaryDiffEq = load_julia_package("OrdinaryDiffEq") sys.modules[__name__] = OrdinaryDiffEq # mutate myself diff --git a/diffeqpy/setup.jl b/diffeqpy/setup.jl deleted file mode 100644 index 1d92213..0000000 --- a/diffeqpy/setup.jl +++ /dev/null @@ -1,32 +0,0 @@ -@debug "Importing DiffEqBase.jl...." -try - import DiffEqBase -catch err - @error "Failed to import DiffEqBase.jl" exception = (err, catch_backtrace()) - rethrow() -end - -import PythonCall - -# TODO: upstream this into a package extension between PythonCall and DiffEqBase -function DiffEqBase.numargs(f::PythonCall.Py) - inspect = PythonCall.pyimport("inspect") - PythonCall.hasproperty(f,:py_func) ? _f = f.py_func : _f = f - nargs = length(first(inspect.getfullargspec(_f))) - if PythonCall.pyconvert(Bool, inspect.ismethod(_f)) - # `f` is a bound method (i.e., `self.f`) but `getfullargspec` - # includes `self` in the `args` list. So, we are subtracting - # 1 manually here: - return nargs - 1 - else - return nargs - end -end - -@debug "Importing DifferentialEquationsjl...." -try - import DifferentialEquations -catch err - @error "Failed to import DifferentialEquations.jl" exception = (err, catch_backtrace()) - rethrow() -end From 65980e3fe4321498e50f21d68f0fcf358fefd613 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Mon, 9 Oct 2023 19:19:01 -0500 Subject: [PATCH 06/14] get to working state (squash with last) --- diffeqpy/__init__.py | 9 +++++++++ diffeqpy/de.py | 18 ++---------------- diffeqpy/install.jl | 8 ++------ diffeqpy/ode.py | 12 ++---------- 4 files changed, 15 insertions(+), 32 deletions(-) diff --git a/diffeqpy/__init__.py b/diffeqpy/__init__.py index 511ae03..95ccf31 100644 --- a/diffeqpy/__init__.py +++ b/diffeqpy/__init__.py @@ -35,3 +35,12 @@ def _ensure_installed(*kwargs): if not _find_julia(): # TODO: this should probably ensure that packages are installed too install(*kwargs) + +# TODO: upstream this function or an alternative into juliacall +def load_julia_package(name): + # This is terrifying to many people. However, it seems SciML takes pragmatic approach. + _ensure_installed() + + # Must be loaded after `_ensure_installed()` + from juliacall import Main + return Main.seval(f"using {name}: {name}; {name}") diff --git a/diffeqpy/de.py b/diffeqpy/de.py index 4b666d0..42f9d63 100644 --- a/diffeqpy/de.py +++ b/diffeqpy/de.py @@ -1,17 +1,3 @@ -import os import sys - -from . import _ensure_installed - -# This is terrifying to many people. However, it seems SciML takes pragmatic approach. -_ensure_installed() - -# PyJulia have to be loaded after `_ensure_installed()` -from juliacall import Main - -# TODO: find a better way to do this or upstream this function into juliacall -def load_julia_package(name): - return Main.seval(f"using {name}: {name}; {name}") - -DifferentialEquations = load_julia_package("DifferentialEquations") -sys.modules[__name__] = DifferentialEquations # mutate myself +from . import load_julia_package +sys.modules[__name__] = load_julia_package("DifferentialEquations") # mutate myself diff --git a/diffeqpy/install.jl b/diffeqpy/install.jl index bb02417..9272ea8 100644 --- a/diffeqpy/install.jl +++ b/diffeqpy/install.jl @@ -1,7 +1,3 @@ using Pkg -Pkg.add("DifferentialEquations") -Pkg.add("OrdinaryDiffEq") -Pkg.add("PythonCall") -Pkg.build("PythonCall") -using DifferentialEquations -using PythonCall +Pkg.add(["DifferentialEquations", "OrdinaryDiffEq"]) +using DifferentialEquations, OrdinaryDiffEq # Precompile diff --git a/diffeqpy/ode.py b/diffeqpy/ode.py index c3ab6e7..9557672 100644 --- a/diffeqpy/ode.py +++ b/diffeqpy/ode.py @@ -1,11 +1,3 @@ -import os import sys - -from juliacall import Main - -# TODO: find a better way to do this or upstream this function into juliacall -def load_julia_package(name): - return Main.seval(f"using {name}: {name}; {name}") - -OrdinaryDiffEq = load_julia_package("OrdinaryDiffEq") -sys.modules[__name__] = OrdinaryDiffEq # mutate myself +from . import load_julia_package +sys.modules[__name__] = load_julia_package("OrdinaryDiffEq") # mutate myself From 22c18529e91b1d5153f1aeb16873e05442a918ed Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Mon, 9 Oct 2023 19:19:30 -0500 Subject: [PATCH 07/14] revert all README changes --- README.md | 41 +++++++++++++++++------------------------ 1 file changed, 17 insertions(+), 24 deletions(-) diff --git a/README.md b/README.md index 7bb80a5..4c256d2 100644 --- a/README.md +++ b/README.md @@ -75,9 +75,7 @@ translate these docs to Python code. Python does not allow `!` in function names, so this is also [a limitation of pyjulia](https://pyjulia.readthedocs.io/en/latest/limitations.html#mismatch-in-valid-set-of-identifiers) To use functions which on the Julia side have a `!`, like `step!`, replace `!` by `_b`, for example: -TODO: - -```julia +```py from diffeqpy import de def f(u,p,t): @@ -99,12 +97,11 @@ is valid Python code for using [the integrator interface](https://diffeq.sciml.a ```py from diffeqpy import de -import numpy as np -def f(du,u,p,t): - du[0] = -u[0] +def f(u,p,t): + return -u -u0 = np.array([0.5]) +u0 = 0.5 tspan = (0., 1.) prob = de.ODEProblem(f, u0, tspan) sol = de.solve(prob) @@ -121,7 +118,7 @@ We can plot the solution values using matplotlib: ```py import matplotlib.pyplot as plt -plt.plot(sol.t,[u[0] for u in sol.u]) +plt.plot(sol.t,sol.u) plt.show() ``` @@ -133,7 +130,7 @@ We can utilize the interpolation to get a finer solution: import numpy t = numpy.linspace(0,1,100) u = sol(t) -plt.plot(t,[ui[0] for ui in u]) +plt.plot(t,u) plt.show() ``` @@ -169,12 +166,13 @@ sol = de.solve(prob) ``` Additionally, you can directly define the functions in Julia. This will allow -for more specialization and could be helpful to increase the efficiency over the -Numba version for repeat or long calls. This is done via `de.seval` or `ode.seval`: +for more specialization and could be helpful to increase the efficiency over +the Numba version for repeat or long calls. This is done via `julia.Main.eval`: ```py -jul_f = de.seval("(u,p,t)->-u") # Define the anonymous function in Julia -prob = de.ODEProblem(jul_f, u0[0], tspan) +from julia import Main +jul_f = Main.eval("(u,p,t)->-u") # Define the anonymous function in Julia +prob = de.ODEProblem(jul_f, u0, tspan) sol = de.solve(prob) ``` @@ -186,21 +184,18 @@ To solve systems of ODEs, simply use an array as your initial condition and define `f` as an array function: ```py -def f(du,u,p,t): +def f(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 + return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z] -u0 = np.array([1.0,0.0,0.0]) +u0 = [1.0,0.0,0.0] tspan = (0., 100.) -p = (10.0,28.0,8/3) +p = [10.0,28.0,8/3] prob = de.ODEProblem(f, u0, tspan, p) sol = de.solve(prob,saveat=0.01) -u = np.array([list(ui) for ui in sol.u]) -plt.plot(sol.t,u) +plt.plot(sol.t,sol.u) plt.show() ``` @@ -209,7 +204,7 @@ plt.show() or we can draw the phase plot: ```py -ut = numpy.transpose(u) +ut = numpy.transpose(sol.u) from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') @@ -221,8 +216,6 @@ plt.show() ### In-Place Mutating Form -TODO - When dealing with systems of equations, in many cases it's helpful to reduce memory allocations by using mutating functions. In diffeqpy, the mutating form adds the mutating vector to the front. Let's make a fast version of the From 10d5b1888eae523cd21ac9df58a9d5d2128b8235 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Mon, 9 Oct 2023 20:41:40 -0500 Subject: [PATCH 08/14] remove TODO --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 807d8c8..01fef8b 100644 --- a/setup.py +++ b/setup.py @@ -24,6 +24,6 @@ def readme(): author_email='contact@juliadiffeq.org', license='MIT', packages=['diffeqpy','diffeqpy.tests'], - install_requires=['juliacall>=0.9.14', 'jill'], # TODO: specify looser compat bound for juliacall if safe. + install_requires=['juliacall>=0.9.14', 'jill'], include_package_data=True, zip_safe=False) From 64171fa50ab06852133e2394bcc2b4a064b35ffd Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Mon, 9 Oct 2023 20:35:51 -0500 Subject: [PATCH 09/14] Make changes to README so it still works (BREAKING CHANGES) --- README.md | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index 4c256d2..730da11 100644 --- a/README.md +++ b/README.md @@ -167,11 +167,10 @@ sol = de.solve(prob) Additionally, you can directly define the functions in Julia. This will allow for more specialization and could be helpful to increase the efficiency over -the Numba version for repeat or long calls. This is done via `julia.Main.eval`: +the Numba version for repeat or long calls. This is done via `seval`: ```py -from julia import Main -jul_f = Main.eval("(u,p,t)->-u") # Define the anonymous function in Julia +jul_f = de.seval("(u,p,t)->-u") # Define the anonymous function in Julia prob = de.ODEProblem(jul_f, u0, tspan) sol = de.solve(prob) ``` @@ -195,7 +194,7 @@ p = [10.0,28.0,8/3] prob = de.ODEProblem(f, u0, tspan, p) sol = de.solve(prob,saveat=0.01) -plt.plot(sol.t,sol.u) +plt.plot(sol.t,de.transpose(de.stack(sol.u))) plt.show() ``` @@ -204,11 +203,11 @@ plt.show() or we can draw the phase plot: ```py -ut = numpy.transpose(sol.u) +us = de.stack(sol.u) from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') -ax.plot(ut[0,:],ut[1,:],ut[2,:]) +ax.plot(us[0,:],us[1,:],us[2,:]) plt.show() ``` @@ -240,7 +239,7 @@ sol = de.solve(prob) or using a Julia function: ```py -jul_f = Main.eval(""" +jul_f = de.seval(""" function f(du,u,p,t) x, y, z = u sigma, rho, beta = p @@ -275,7 +274,7 @@ tspan = (0.0,1.0) prob = de.SDEProblem(f,g,u0,tspan) sol = de.solve(prob,reltol=1e-3,abstol=1e-3) -plt.plot(sol.t,sol.u) +plt.plot(sol.t,de.stack(sol.u)) plt.show() ``` @@ -430,14 +429,14 @@ the solver accuracy by accurately stepping at the points of discontinuity. Together this is: ```py -f = Main.eval(""" +f = de.seval(""" function f(du, u, h, p, t) du[1] = 1.1/(1 + sqrt(10)*(h(p, t-20)[1])^(5/4)) - 10*u[1]/(1 + 40*u[2]) du[2] = 100*u[1]/(1 + 40*u[2]) - 2.43*u[2] end""") u0 = [1.05767027/3, 1.030713491/3] -h = Main.eval(""" +h = de.seval(""" function h(p,t) [1.05767027/3, 1.030713491/3] end From 80af3ddbe2a58719cd4c954bbca0fd0fbc8e2730 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Mon, 9 Oct 2023 21:07:06 -0500 Subject: [PATCH 10/14] fix a test failure BREAKING --- diffeqpy/tests/test_dde.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/diffeqpy/tests/test_dde.py b/diffeqpy/tests/test_dde.py index 315e400..b6e9a6d 100644 --- a/diffeqpy/tests/test_dde.py +++ b/diffeqpy/tests/test_dde.py @@ -4,14 +4,14 @@ def test(): - f = Main.eval(""" + f = de.seval(""" function f(du, u, h, p, t) du[1] = 1.1/(1 + sqrt(10)*(h(p, t-20)[1])^(5/4)) - 10*u[1]/(1 + 40*u[2]) du[2] = 100*u[1]/(1 + 40*u[2]) - 2.43*u[2] end""") u0 = [1.05767027/3, 1.030713491/3] - h = Main.eval(""" + h = de.seval(""" function h(p,t) [1.05767027/3, 1.030713491/3] end From f9d66533dc1c10fe433fddcb29c8ed893cbb7f80 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 10 Oct 2023 07:47:13 -0500 Subject: [PATCH 11/14] use @diffeqpy instead of @pyjuliapkg for enviornmnet and add and load PythonCall in install script to precompile extensions --- diffeqpy/__init__.py | 2 +- diffeqpy/install.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/diffeqpy/__init__.py b/diffeqpy/__init__.py index 95ccf31..ddc6be6 100644 --- a/diffeqpy/__init__.py +++ b/diffeqpy/__init__.py @@ -43,4 +43,4 @@ def load_julia_package(name): # Must be loaded after `_ensure_installed()` from juliacall import Main - return Main.seval(f"using {name}: {name}; {name}") + return Main.seval(f"import Pkg; Pkg.activate(\"diffeqpy\", shared=true); import {name}; {name}") diff --git a/diffeqpy/install.jl b/diffeqpy/install.jl index 9272ea8..03ddb76 100644 --- a/diffeqpy/install.jl +++ b/diffeqpy/install.jl @@ -1,3 +1,4 @@ using Pkg -Pkg.add(["DifferentialEquations", "OrdinaryDiffEq"]) -using DifferentialEquations, OrdinaryDiffEq # Precompile +Pkg.activate("diffeqpy", shared=true) +Pkg.add(["DifferentialEquations", "OrdinaryDiffEq", "PythonCall"]) +using DifferentialEquations, OrdinaryDiffEq, PythonCall # Precompile From 9358a7a84c3757b3ca081686b05ce8af3af4445f Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 10 Oct 2023 07:48:09 -0500 Subject: [PATCH 12/14] remove now-unnecessary line in tests --- diffeqpy/tests/test_dde.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/diffeqpy/tests/test_dde.py b/diffeqpy/tests/test_dde.py index b6e9a6d..aab6262 100644 --- a/diffeqpy/tests/test_dde.py +++ b/diffeqpy/tests/test_dde.py @@ -1,5 +1,3 @@ -from juliacall import Main - from .. import de From 155099e5f2b2c472b96c121fcdfbb755994e0395 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 10 Oct 2023 07:50:31 -0500 Subject: [PATCH 13/14] remove trailing whitespace in test to work around https://github.com/JuliaPy/PythonCall.jl/issues/379 --- diffeqpy/tests/test_dde.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/diffeqpy/tests/test_dde.py b/diffeqpy/tests/test_dde.py index aab6262..62d467a 100644 --- a/diffeqpy/tests/test_dde.py +++ b/diffeqpy/tests/test_dde.py @@ -12,8 +12,7 @@ def test(): h = de.seval(""" function h(p,t) [1.05767027/3, 1.030713491/3] - end - """) + end""") tspan = (0.0, 100.0) constant_lags = [20.0] From f03075da8cde8a5cfe7c4bced08192152f7aa706 Mon Sep 17 00:00:00 2001 From: Lilith Hafner Date: Tue, 10 Oct 2023 08:20:13 -0500 Subject: [PATCH 14/14] fix numba examples where easy and advertize errors where it's not --- README.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 730da11..fd37d78 100644 --- a/README.md +++ b/README.md @@ -162,7 +162,7 @@ import numba numba_f = numba.jit(f) prob = de.ODEProblem(numba_f, u0, tspan) -sol = de.solve(prob) +sol = de.solve(prob) # ERROR ``` Additionally, you can directly define the functions in Julia. This will allow @@ -309,11 +309,11 @@ sol = de.solve(prob) # Now let's draw a phase plot -ut = numpy.transpose(sol.u) +us = de.stack(sol.u) from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') -ax.plot(ut[0,:],ut[1,:],ut[2,:]) +ax.plot(us[0,:],us[1,:],us[2,:]) plt.show() ``` @@ -358,11 +358,11 @@ sol = de.solve(prob,saveat=0.005) # Now let's draw a phase plot -ut = numpy.transpose(sol.u) +us = de.stack(sol.u) from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') -ax.plot(ut[0,:],ut[1,:],ut[2,:]) +ax.plot(us[0,:],us[1,:],us[2,:]) plt.show() ``` @@ -411,7 +411,7 @@ def f(resid,du,u,p,t): numba_f = numba.jit(f) prob = de.DAEProblem(numba_f,du0,u0,tspan,differential_vars=differential_vars) -sol = de.solve(prob) +sol = de.solve(prob) # ERROR ``` ## Delay Differential Equations