nicholaswogan / numbalsoda Goto Github PK

View Code? Open in Web Editor NEW

91.0 91.0 9.0 297 KB

Python wrapper of LSODA (solving ODEs) which can be called from within numba functions.

License: MIT License

CMake 0.33% Python 3.99% Jupyter Notebook 61.27% C++ 17.05% Julia 0.44% Fortran 16.92%

numba ode-solver odepack python

numbalsoda's People

Contributors

Stargazers

Watchers

Forkers

chrisrackauckas geraintpratten tiagomiguelrs vgurev tbridel joeyshuttleworth tobias-fischer ahsoka genesofeve

numbalsoda's Issues

elementwise vector operation and passing size of the system

Hello,

I am new to numblsoda so apologies if I missed something obvious. But I am trying to solve a system of N equations with N being an user input. Is there any way to pass N to the step function(rhs in the docs) written below? Also, is there any way to do vector operations without writing for loops?

Here is my minimal working code in scipy:

N=1000
km=1+np.log(2)/3
kp=1.1*km
kbm=10**-6
T=24
t=np.arange(0,T+1,1)

d=np.array([-((kp*k/N + kbm)*(N - k) + km*k) for k in range(0,N+1)])
u=np.array([km*(k + 1) for k in range(0,N)])
l=np.array([(kp*(k - 1)/N + kbm)*(N - k + 1) for k in range(1,N+1)])

pi=np.zeros(N+1, dtype=np.float64)
pi[N]=1
dp=np.zeros(N+1)

def step(p, t):
      dp[0]=d[0]*p[0]+u[0]*p[1]
      dp[1:N]=d[1:N]*p[1:N]+u[1:]*p[2:]+l[:N-1]*p[:N-1]
      dp[N]=d[N]*p[N]+l[N-1]*p[N-1]
      return dp

sol = odeint(step, pi, t, mxstep=2000)

Thanks.

Easier passing data to rhs

Thanks for the package, some 400x faster than plain Python for the system I'm looking at. I just wanted to point out that Numba functions close over their local scope so passing arbitrary (constant?) data can be much easier than in your example. Mine looks like this:

def make_lsoda_func(c,R,dstar,E,F,N):
    @cfunc(lsoda_sig)
    def rhs(t, x, du, p):
        codim3(t,x,c,R,dstar,E,F,N,du)
    return rhs

Cannot import dop853 from numbalsoda

Hi!
I've just installed your library and got this error when running the code shown in the readme.
ImportError: cannot import name 'dop853' from 'numbalsoda'
Could I be missing a library? Thanks in advance

Numba cache and iterative solving

Hi,

I am trying to use NumbaLSODA to solve many systems iteratively.
However, it seems that I can't out the function in the numba cache and so I have a big overhead at each iteration.
I was wondering if you have an advice about it.

This is my code.
First I put all numba related function inside a file

numba_func.py

NOPYTHON = True
CACHE = True
PARALLEL = True

@nb.cfunc(lsoda_sig, nopython=NOPYTHON, cache=CACHE)
def rhs(x, i, di, p):
    di[0] = ...
    di[1] = ...
    di[2] = ...

funcptr = rhs.address

@nb.njit('(float64)(float64, float64, float64, int16)', nopython=NOPYTHON, parallel=PARALLEL, cache=CACHE)
def solve(a, b, c, funcptr):

    w0 = np.array([a, b, ...], dtype=np.float64)
    p  = np.array([c], dtype=np.float64)

    x = np.linspace(0, 100, 500)

    usol, success = lsoda(funcptr, w0, x, data=p)
    
    return usol[-1][1]

Then I use another file to solve the systems one after the others

from numba_func import solve, funcptr

gs = []
for a, b, c in zip(as, bs, cs):
    gs = np.append(gs, solve(a, b, c, funcptr))

I got the following warning:

NumbaWarning: Cannot cache compiled function "solve" as it uses dynamic globals (such as ctypes pointers and large global arrays)

I guess the idea is to pass the variable funcptr correctly so that numba is happy but so far I failed...

Installation via pip produces error when importing numbalsoda

Hi Nick, thank you very much for building this Python package, it definitely appears very promising!

I installed the latest version of numbalsoda (0.3.4) via pip and the installation was successful. I am using Python 3.9.7 and conda version 22.9.0.

However, when trying to import it, I get the following error:

'The procedure entry point ZSt28__throw_bad_array_new_lengthv could not be located in the dynamic link library C:\Users...\numbalsoda\liblsoda.dll'

Installing numbalsoda via conda (conda install -c conda-forge numbalsoda) solves the issue and the package is successfully imported into the Python script. Conda installs version 0.3.2, and thus I was wondering if the difference in the versions may play a role in this issue.

Thank you very much!

Make NumPy types for u0 and usol configurable

First of all, awesome code! I got a speedup of an order of magnitude basically for free!

But, as far as I can tell, u0 and usol are required to be of type np.float64. However, often times np.float32 is more than enough.

Maybe I just don't see the option to change it, but if it doesn't exist, it'd be nice if we could change it.

Thank you!

Timeout for long integration

Is it possible to add some timeout for the integration? I've been solving a stiff problem and the solving takes hours without stopping. I think it would be better to halt it if 10 or 20s are exceeded.

What about people being stuck with windows at work 😓

Hi,

Thanks for this package which make things so much faster!
Is there a way to make your pip install available to Windows user?
For instance by making a conda package from it?

Unexpected behavior of NumbaLSODA and nb.cufnc in parallel mode

Hi There,
Thanks for sharing this package. I am seeing up to 20x speedup in my code. As I am simulating the evolution of multiple initial conditions, I want to speed up the code further by using parallel=True flag in Numba. I adapted the example you provided in Stackoverflow for my code to use NumbaLSODA with parallel=True flag in Numba. However, I observed that the speedup was marginal even with 8 cores . In some cases, my parallel code was actually slower (!).

After debugging on the 2 state example using the code below, I narrowed the reason down to the way the parameters in the rhs function are defined. if the parameters are either global variables or passed using the data argument in lsoda, the parallel = True speeds up the code almost linearly versus the number of cores. However, if the parameters are defined locally as in f_local below, the speed up is marginal. In fact, the parallel version using f_local is slower than the series version using f_global or f_param.

I thought that local declarations of variables is a good coding practice and it speeds up code execution, but this does not seem to be the case with Numba and NumbaLSODA. I do not know if this behavior is caused by cfunc in Numba or NumbaLSODA. It was definitely unexpected and I thought I will bring it to your notice. Do you know the reason for this behavior? Are local constants not compiled just as global constants? I am unable to dig deeper into the reason and was wondering if you could help.
Best
Amar

from NumbaLSODA import lsoda_sig, lsoda
from matplotlib import pyplot as plt
import numpy as np
import numba as nb
import time

a_glob=np.array([1.5,1.5])

@nb.cfunc(lsoda_sig)
def f_global(t, u_, du, p): # variable a is global
    u = nb.carray(u_, (2,))
    du[0] = u[0]-u[0]*u[1]*a_glob[0]
    du[1] = u[0]*u[1]-u[1]*a_glob[1]

    
@nb.cfunc(lsoda_sig)
def f_local(t, u_, du, p): # variable a is local
    u = nb.carray(u_, (2,))
    a = np.array([1.5,1.5]) 
    du[0] = u[0]-u[0]*u[1]*a[0]
    du[1] = u[0]*u[1]-u[1]*a[1]

    
@nb.cfunc(lsoda_sig)
def f_param(t, u_, du, p): # pass in a as a paameter
    u = nb.carray(u_, (2,))
    du[0] = u[0]-u[0]*u[1]*p[0]
    du[1] = u[0]*u[1]-u[1]*p[1]

funcptr_glob = f_global.address
funcptr_local = f_local.address
funcptr_param = f_param.address
t_eval = np.linspace(0.0,20.0,201)
np.random.seed(0)
a = np.array([1.5,1.5])

@nb.njit(parallel=True)
def main(n, flag):
#     a = np.array([1.5,1.5])
    u1 = np.empty((n,len(t_eval)), np.float64)
    u2 = np.empty((n,len(t_eval)), np.float64)
    for i in nb.prange(n):
        u0 = np.empty((2,), np.float64)
        u0[0] = np.random.uniform(4.5,5.5)
        u0[1] = np.random.uniform(0.7,0.9)
        if flag ==1: # global
            usol, success = lsoda(funcptr_glob, u0, t_eval, rtol = 1e-8, atol = 1e-8)
        if flag ==2: # local
            usol, success = lsoda(funcptr_local, u0, t_eval, rtol = 1e-8, atol = 1e-8)
        if flag ==3: # param
            usol, success = lsoda(funcptr_param, u0, t_eval, data = a, rtol = 1e-8, atol = 1e-8)
        u1[i] = usol[:,0]
        u2[i] = usol[:,1]
    return u1, u2

@nb.njit(parallel=False)
def main_series(n, flag): # same function as above but with parallel flag = False
#     a = np.array([1.5,1.5])
u1 = np.empty((n,len(t_eval)), np.float64)
u2 = np.empty((n,len(t_eval)), np.float64)
for i in nb.prange(n):
    u0 = np.empty((2,), np.float64)
    u0[0] = np.random.uniform(4.5,5.5)
    u0[1] = np.random.uniform(0.7,0.9)
    if flag ==1: # global
        usol, success = lsoda(funcptr_glob, u0, t_eval, rtol = 1e-8, atol = 1e-8)
    if flag ==2: # local
        usol, success = lsoda(funcptr_local, u0, t_eval, rtol = 1e-8, atol = 1e-8)
    if flag ==3: # param
        usol, success = lsoda(funcptr_param, u0, t_eval, data = a, rtol = 1e-8, atol = 1e-8)
    u1[i] = usol[:,0]
    u2[i] = usol[:,1]
return u1, u2

n = 100
# calling first time for JIT compiling
u1, u2 = main(n,1)
u1, u2 = main(n,2)
u1, u2 = main(n,3)

u1, u2 = main_series(n,1)
u1, u2 = main_series(n,1)
u1, u2 = main_series(n,1)

# Running code for large number 
n = 10000
a1 = time.time()
u1, u2 = main(n,1) # global
a2 = time.time()
print(a2 - a1) # this is fast


a1 = time.time()
u1, u2 = main(n,2) # local
a2 = time.time()
print(a2 - a1) # this is slow

a1 = time.time()
u1, u2 = main(n,3) # param
a2 = time.time()
print(a2 - a1) # this is fast and almost identical performance as global

a1 = time.time()
u1, u2 = main_series(n,1) # global
a2 = time.time()
print(a2 - a1) # this is faster than local + parallel

a1 = time.time()
u1, u2 = main_series(n,2) # local
a2 = time.time()
print(a2 - a1) # this is slow

a1 = time.time()
u1, u2 = main_series(n,3) # param
a2 = time.time()
print(a2 - a1) # this is fast and almost identical performance as global

[lsoda] 500 steps taken before reaching tout

I am running a simulation and if I understand correctly the maximum number of steps (500) is performed before the integration finishes. Is there a way to change this parameter?

Thanks

Backward in time using numbalsoda

I am using numblsoda to solve couples of first order differential equations. I need to solve those equations backward in time for some particular situations. I could not figure out how I can do that. I tried to use decreasing time step but it does not work as it works for solve_ivp and odeint.

for instance in this code how can I go backward in time?

from numbalsoda import lsoda_sig, lsoda, dop853
from numba import njit, cfunc
import numpy as np

@cfunc(lsoda_sig)
def rhs(t, u, du, p):
    du[0] = u[0]-u[0]*u[1]
    du[1] = u[0]*u[1]-u[1]*p[0]

funcptr = rhs.address # address to ODE function
u0 = np.array([5.,0.8]) # Initial conditions
data = np.array([1.0]) # data you want to pass to rhs (data == p in the rhs).
t_eval = np.arrange(0.0,50.0,0.1) # times to evaluate solution


usol, success = lsoda(funcptr, u0, t_eval, data = data)

I tried t_eval = np.arrange(50, 0, -0.1) but it does not work. It would be great if someone could help me to fix this issue.

solve_ivp is missing or not added in conda package

In installed numbalsoda using 'conda install' but importing solve_ivp failed. Comparing the content of the packages, installed via conda and pip, discovered that file /driver_solve_ivp.py is simply missing
Traceback (most recent call last):
File "<...>\test_numbalsoda.py", line 1, in
from numbalsoda import lsoda_sig, lsoda, dop853, solve_ivp
ImportError: cannot import name 'solve_ivp' from 'numbalsoda' (<...>\envs\conda_python310\lib\site-packages\numbalsoda_init.py)_

I encountered this issue while trying install the package via conda install to bypass the issue #29.
Here is the folder with conda package:

Request for PyPI package

First of all, thanks for creating this great library! I'm developing a framework for cosmological simulations, and NumbaLSODA speeds up our ODE integration by about 20x compared to using scipy.integrate.odeint() with pure Python, and even more compared to scipy.integrate.solve_ivp().

Would it be possible for you to push NumbaLSODA wheels to PyPI? This would simplify our installation process, as the support for Git dependencies has been dropped from setuptools.

Could not find module

After running the following in Jupiter Notebook:

pip install numbalsoda

Collecting numbalsoda
  Downloading numbalsoda-0.3.4-cp39-cp39-win_amd64.whl (151 kB)
     -------------------------------------- 151.6/151.6 kB 4.4 MB/s eta 0:00:00
Requirement already satisfied: numpy in c:\users\[user]\anaconda3\lib\site-packages (from numbalsoda) (1.21.5)
Requirement already satisfied: numba in c:\users\[user]\anaconda3\lib\site-packages (from numbalsoda) (0.55.1)
Requirement already satisfied: llvmlite<0.39,>=0.38.0rc1 in c:\users\[user]\anaconda3\lib\site-packages (from numba->numbalsoda) (0.38.0)
Requirement already satisfied: setuptools in c:\users\[user]\anaconda3\lib\site-packages (from numba->numbalsoda) (63.4.1)
Installing collected packages: numbalsoda
Successfully installed numbalsoda-0.3.4
Note: you may need to restart the kernel to use updated packages.

It is not possible for me to import the package, and I get the following error message:

from numbalsoda import lsoda_sig, lsoda, dop853

FileNotFoundError: Could not find module 'C:\Users\[user]\[folder]\Gitlab\numbalsoda\numbalsoda\liblsoda.dll' (or one of its dependencies). Try using the full path with constructor syntax.

Can I use 2d-array of 'u0'?

I want to use the 2d-shape initial values of the multiple-particles system. But, I think I cannot use the slicing or indexing.

def swing_lsoda(network):
    """
    rhs = swing_cover(network)
    funcptr = rhs.address
    ...
    data = np.array([1.0])
    usol, success = lsoda(funcptr, u0, t_eval, data = data)
    """
    @nb.cfunc(lsoda_sig) # network
    def rhs(t, u, du, p): # p0=m, p1=gamma, p2=P, p3=K
        # u is 2d-array 
        Interaction = p[3] * SinIntE(network, u[0])
        du[0] = u[1] #Thetas of nodes
        du[1] = 1/p[0]*(p[2] - p[1]*u[1] - Interaction) #Omega

How can I try it?

do you have a tip how to workaround the global variable?

Hi,

thanks for the great code, it's very fast.

I'm struggling though with a need for a global variable, or maybe an argument that is updated by the ODE itself

@nb.cfunc(lsoda_sig)
def particle_ode(t, y_, dydt, p_):
    y = nb.carray(y_,(3,))
    p = nb.carray(p_,(11,))
    rho1, rho, nu , Cd, g, rhop, d, zu, zl, Vc0, trec = p
   
    global tzl # <---- here 

                           
    zp   = y[0]                # particle position      [m]
    V    = y[1]                # particle velocity      [m/s]
    tzl =  y[2] 
    Cam = 0.5                  # added mass coefficient [-]

    # determine Vc
    if (y[0]<= zl):
        Vc = Vc0
        tzl = t     # < ----- here I need to adjust this as we move through with y[0]
    else:
        Vc = Vc0 * np.exp(-1*((t-tzl)/trec))
    FS = (rho(zp) - rho1)*Vc*g

    # force balance on the particle 
    dzpdt = V
    dVdt  = ( (rhop-rho(zp)) * Vp(d) * g \
            - 0.5 *Cd(V*d/nu(zp))* rho(zp) * Ap(d) * np.abs(V) * V \
            - FS \
            ) / (rhop*Vp(d) + Cam*rho(zp)*Vp(d))

    dydt = [dzpdt, dVdt, 0.0 ]
funcptr = particle_ode.address

Support for other integrators

Hi, thanks for this great package!

I know this package is only supposed to deal with the lsoda integrator, but I was wondering if you could at least add support for an explicit integrator as well, maybe rk23. The reason I ask is that I want to solve a very large system of ODEs (>10^6), and because of that, implicit integrators are not really suitable (since they require a matrix solve Ax=b at each step), so I'm stuck with scipy's rk23 or rk45, which are fine, except that because the integration involves too many time steps, I think solve_ivp's explicit while loop is becoming the bottle-neck.

Adaptive timestepping and hijacking the "step" method

Hi again @Nicholaswogan !

I was wondering if you ever thought about the adaptive timestepping capability that exists in Scipy LSODA ? Is anything like this possible with your implementation ?

Thanks !
Best,
T. Bridel-Bertomeu

Strange results when t_eval has repeat values

Thanks for the great package!

Summary:
I'm trying to switch over from using scipy's odeint to NumbaLSODA as it seems much better! Sometimes the t_eval array I supply has repeated times at the end and this produces some strange output. This was okay with odeint though, would it be possible to allow it for NumbaLSODA too? I think it is trying to integrate an infinitely small timestep rather than just evaluating the same timestep multiple times.

Details:
I'm trying to add this to LEGWORK so that we can speed up how we evolve the orbits of binary stars due to gravitational wave emission. However, if the evolution gets too close to the merger then LSODA produces a bunch of warnings about the timesteps being too small (since the derivatives get so large). To avoid this, I set any timesteps that are too close to the merger equal to the last allowable timestep - this leads to repeated values and the weird results that I see. (I don't just trim the array since I do this in a 2D array with a collection of binaries and so each individual binary might have a different length array which would wouldn't work in a 2D array I think.)

MWE:
Evolve binary evolution due to gravitational wave emission using Peters (1964)

from numba import cfunc
import numba as nb
import numpy as np
from NumbaLSODA import lsoda_sig, lsoda

@cfunc(lsoda_sig)
def de_dt(t, e, de, data):
    e_ = nb.carray(e, (1,))
    beta, c_0 = nb.carray(data, (2,))
    de[0] = -19 / 12 * beta / c_0**4 * (e_[0]**(-29.0 / 19.0) * (1 - e_[0]**2)**(3.0/2.0)) \
        / (1 + (121./304.) * e_[0]**2)**(1181./2299.)

ecc = 0.5
beta = 2.47e22
c_0 = 1.677e10
t_merge = 1.81e17

funcptr = de_dt.address
u0 = np.array([ecc])
data = np.array([beta, c_0])
t_eval = np.array([0.0, 0.1, 0.2, 0.5, 0.8, 0.9, 0.9, 0.9]) * t_merge

ecc_evol, _ = lsoda(funcptr, u0, t_eval=t_eval, data=data)
print(ecc_evol)

This outputs

[[ 5.00000000e-01]
 [ 4.83244343e-01]
 [ 4.64840518e-01]
 [ 3.95592532e-01]
 [ 2.82404708e-01]
 [ 2.15245963e-01]
 [-1.36311572e+57]
 [-1.36311572e+57]]

So you can see that as soon as the timesteps repeat, you get this hugely small number (and eccentricity should be between 0 and 1)

Potential interpolant issue

Hi @Nicholaswogan,

Great work on this repository, it looks like you managed a consistent and quite high speed-up compared to SciPy. I have been working with LSODA recently using solve_ivp, and I have encountered an issue relative to the interpolant (PR at SciPy). I am wondering if this implementation of LSODA also has a similar issue, as it looks to me that the root of the problem lies within the FORTRAN implementation. In particular, LSODA::intdy seems very similar to me to DINTDY in ODEPACK, which could well mean that the issue is present. Let me know if you can have a look.

Using numbalsoda solve_ivp

Hi,

I would like to solve a stiff ODE with the solve_ivp wrapper by numbalsoda - as you did here. Unfortunatly i can't load the module solve_ivp from numbalsoda (not available in the package).

Could your provide the module solve_ivp?

Discrepancy between lsoda and standard scipy

Hello,

I got some weird behavior using LSODA.
I compared them to scipy ode and scipy looked right while LSODA is off.
I wrote a test code to demonstrate my point:

# This Python file uses the following encoding: utf-8

import numba as nb
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
from NumbaLSODA import lsoda_sig, lsoda



@nb.cfunc(lsoda_sig)
def rhs(x: float,
        i: tuple,
        di: tuple,
        p: tuple) -> None:

    # Reconstruct the complex number
    i0 = i[0] + 1j*i[1]
    i1 = i[2] + 1j*i[3]
    i2 = i[4] + 1j*i[5]

    kp     = p[0]
    ks     = p[1]
    ki     = p[2]
    ep     = p[3]
    xi     = p[4]

    delta = kp-ks-ki

    di0 = 1j*kp*ep*i1*i2*            np.exp(-1j*delta*x)/4 + 1j*kp*xi*i0/8*(  abs(i0)**2 + 2*abs(i1)**2 + 2*abs(i2)**2)
    di1 = 1j*ks*ep*i0*i2.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ks*xi*i1/8*(2*abs(i0)**2 +   abs(i1)**2 + 2*abs(i2)**2)
    di2 = 1j*ki*ep*i0*i1.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ki*xi*i2/8*(2*abs(i0)**2 + 2*abs(i1)**2 +   abs(i2)**2)


    di[0] = di0.real
    di[1] = di0.imag
    di[2] = di1.real
    di[3] = di1.imag
    di[4] = di2.real
    di[5] = di2.imag
funcptr = rhs.address

def vectorfield(x: float,
                i: tuple,
                p: tuple) -> None:

    # Reconstruct the complex number
    i0 = i[0]
    i1 = i[1]
    i2 = i[2]

    kp     = p[0]
    ks     = p[1]
    ki     = p[2]
    ep     = p[3]
    xi     = p[4]

    delta = kp-ks-ki

    di0 = 1j*kp*ep*i1*i2*            np.exp(-1j*delta*x)/4 + 1j*kp*xi*i0/8*(  abs(i0)**2 + 2*abs(i1)**2 + 2*abs(i2)**2)
    di1 = 1j*ks*ep*i0*i2.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ks*xi*i1/8*(2*abs(i0)**2 +   abs(i1)**2 + 2*abs(i2)**2)
    di2 = 1j*ki*ep*i0*i1.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ki*xi*i2/8*(2*abs(i0)**2 + 2*abs(i1)**2 +   abs(i2)**2)

    r = np.zeros(3, dtype=np.complex128)
    r[0] = di0
    r[1] = di1
    r[2] = di2

    return r

def currents_lsoda():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0.real, i0.imag,
                   i1.real, i1.imag,
                   i2.real, i2.imag], dtype=np.float64)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 10e-3, 100)

    usol, success = lsoda(funcptr, w0, x, data=p)

    return (usol[:,0] + 1j*usol[:,1],
            usol[:,2] + 1j*usol[:,3],
            usol[:,4] + 1j*usol[:,5])


def currents_scipy():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0,
                   i1,
                   i2], dtype=np.complex128)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 75e-3, 100)
    ys = np.array([])
    xs = np.array([])

    r = ode(vectorfield).set_integrator('zvode').set_f_params(p).set_initial_value(w0, t=x[0])

    dx = x[1]-x[0]
    while r.successful() and r.t-dx <=x[-1]:
        r.integrate(r.t+dx)
        xs = np.append(xs, r.t)
        ys = np.append(ys, r.y)
    ys = np.reshape(ys, (int(len(ys)/len(w0)), len(w0))).T

    return xs, ys



i0_l, i1_l, i2_l = currents_lsoda()
xs, (i0_s, i1_s, i2_s) = currents_scipy()


x = np.linspace(0, 75e-3, len(i0_l))


fig, ax = plt.subplots(3, 2)

ax[0][0].set_title('lsoda')
ax[0][0].plot(x*1e3, np.abs(i0_l)*1e6)
ax[0][0].set_ylabel('i0')

ax[1][0].plot(x*1e3, np.abs(i1_l)*1e3)
ax[1][0].set_ylabel('i1')

ax[2][0].plot(x*1e3, np.abs(i2_l)*1e3)
ax[2][0].set_ylabel('i2')

ax[-1][0].set_xlabel('distance (mm)')


ax[0][1].set_title('scipy')
ax[0][1].plot(xs*1e3, np.abs(i0_s)*1e6)
ax[0][1].set_ylabel('i0')

ax[1][1].plot(xs*1e3, np.abs(i1_s)*1e3)
ax[1][1].set_ylabel('i1')

ax[2][1].plot(xs*1e3, np.abs(i2_s)*1e3)
ax[2][1].set_ylabel('i2')

ax[-1][1].set_xlabel('distance (mm)')

for i in fig.get_axes():
    i.grid('both')
fig.tight_layout()
plt.show()

Which gives:

The right column is the expected results.

Is there something wrong in the way I am using LSODA?

Installation with pip fails

Good evening,

I am working on a project involving the long-time integration of an ODE, and the results you show definitely look very promising, but I cannot install the package ...

I am working on mac os X v12.4 with a M1 chip.
Python v3.10.8, pip v22.3.1, setuptools v65.5.1.

To really test only the install of numbalsoda and remove all possible side effects from pre-installed packages, I put myself in a new virtualenv.
🔴 I am getting the following error when doing pip install numbalsoda:

pip install numbalsoda
Collecting numbalsoda
  Downloading numbalsoda-0.3.4.tar.gz (241 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 241.3/241.3 kB 2.8 MB/s eta 0:00:00
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
  Preparing metadata (pyproject.toml) ... done
Collecting numpy
  Downloading numpy-1.23.5-cp310-cp310-macosx_11_0_arm64.whl (13.4 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 13.4/13.4 MB 9.5 MB/s eta 0:00:00
Collecting numba
  Downloading numba-0.56.4-cp310-cp310-macosx_11_0_arm64.whl (2.4 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.4/2.4 MB 9.4 MB/s eta 0:00:00
Requirement already satisfied: setuptools in ./py3env/lib/python3.10/site-packages (from numba->numbalsoda) (65.5.1)
Collecting llvmlite<0.40,>=0.39.0dev0
  Downloading llvmlite-0.39.1-cp310-cp310-macosx_11_0_arm64.whl (23.1 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 23.1/23.1 MB 10.3 MB/s eta 0:00:00
Building wheels for collected packages: numbalsoda
  Building wheel for numbalsoda (pyproject.toml) ... error
  error: subprocess-exited-with-error
  
  × Building wheel for numbalsoda (pyproject.toml) did not run successfully.
  │ exit code: 1
  ╰─> [73 lines of output]
      /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/setuptools/dist.py:770: UserWarning: Usage of dash-separated 'description-file' will not be supported in future versions. Please use the underscore name 'description_file' instead
        warnings.warn(
      /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/setuptools/config/setupcfg.py:508: SetuptoolsDeprecationWarning: The license_file parameter is deprecated, use license_files instead.
        warnings.warn(msg, warning_class)
      
      
      --------------------------------------------------------------------------------
      -- Trying "Ninja" generator
      --------------------------------
      ---------------------------
      ----------------------
      -----------------
      ------------
      -------
      --
      Not searching for unused variables given on the command line.
      -- The C compiler identification is AppleClang 13.1.6.13160021
      -- Detecting C compiler ABI info
      -- Detecting C compiler ABI info - done
      -- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang - skipped
      -- Detecting C compile features
      -- Detecting C compile features - done
      -- The CXX compiler identification is AppleClang 13.1.6.13160021
      -- Detecting CXX compiler ABI info
      -- Detecting CXX compiler ABI info - done
      -- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang++ - skipped
      -- Detecting CXX compile features
      -- Detecting CXX compile features - done
      -- Configuring done
      -- Generating done
      -- Build files have been written to: /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57/_cmake_test_compile/build
      --
      -------
      ------------
      -----------------
      ----------------------
      ---------------------------
      --------------------------------
      -- Trying "Ninja" generator - success
      --------------------------------------------------------------------------------
      
      Configuring Project
        Working directory:
          /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57/_skbuild/macosx-12.0-arm64-3.10/cmake-build
        Command:
          cmake /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57 -G Ninja -DCMAKE_INSTALL_PREFIX:PATH=/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57/_skbuild/macosx-12.0-arm64-3.10/cmake-install -DPYTHON_VERSION_STRING:STRING=3.10.8 -DSKBUILD:INTERNAL=TRUE -DCMAKE_MODULE_PATH:PATH=/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/skbuild/resources/cmake -DPYTHON_EXECUTABLE:PATH=/Users/tbridel/Documents/1-CODES/py3env/bin/python -DPYTHON_INCLUDE_DIR:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/include/python3.10 -DPYTHON_LIBRARY:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/lib/libpython3.10.dylib -DPython_EXECUTABLE:PATH=/Users/tbridel/Documents/1-CODES/py3env/bin/python -DPython_ROOT_DIR:PATH=/Users/tbridel/Documents/1-CODES/py3env -DPython_INCLUDE_DIR:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/include/python3.10 -DPython_FIND_REGISTRY:STRING=NEVER -DPython3_EXECUTABLE:PATH=/Users/tbridel/Documents/1-CODES/py3env/bin/python -DPython3_ROOT_DIR:PATH=/Users/tbridel/Documents/1-CODES/py3env -DPython3_INCLUDE_DIR:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/include/python3.10 -DPython3_FIND_REGISTRY:STRING=NEVER -DCMAKE_MAKE_PROGRAM:FILEPATH=/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/ninja/data/bin/ninja -DSKBUILD=ON -DCMAKE_BUILD_TYPE:STRING=Release -DCMAKE_OSX_DEPLOYMENT_TARGET:STRING=12.0 -DCMAKE_OSX_ARCHITECTURES:STRING=arm64
      
      CMake Error at /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/cmake/data/share/cmake-3.25/Modules/CMakeDetermineFortranCompiler.cmake:33 (message):
        Could not find compiler set in environment variable FC:
      
        /usr/local/bin/gfortran.
      Call Stack (most recent call first):
        CMakeLists.txt:3 (project)
      
      
      CMake Error: CMAKE_Fortran_COMPILER not set, after EnableLanguage
      CMake Error: CMAKE_CXX_COMPILER not set, after EnableLanguage
      -- Configuring incomplete, errors occurred!
      See also "/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57/_skbuild/macosx-12.0-arm64-3.10/cmake-build/CMakeFiles/CMakeOutput.log".
      Traceback (most recent call last):
        File "/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/skbuild/setuptools_wrap.py", line 632, in setup
          env = cmkr.configure(
        File "/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/skbuild/cmaker.py", line 334, in configure
          raise SKBuildError(
      
      An error occurred while configuring with CMake.
        Command:
          cmake /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57 -G Ninja -DCMAKE_INSTALL_PREFIX:PATH=/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57/_skbuild/macosx-12.0-arm64-3.10/cmake-install -DPYTHON_VERSION_STRING:STRING=3.10.8 -DSKBUILD:INTERNAL=TRUE -DCMAKE_MODULE_PATH:PATH=/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/skbuild/resources/cmake -DPYTHON_EXECUTABLE:PATH=/Users/tbridel/Documents/1-CODES/py3env/bin/python -DPYTHON_INCLUDE_DIR:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/include/python3.10 -DPYTHON_LIBRARY:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/lib/libpython3.10.dylib -DPython_EXECUTABLE:PATH=/Users/tbridel/Documents/1-CODES/py3env/bin/python -DPython_ROOT_DIR:PATH=/Users/tbridel/Documents/1-CODES/py3env -DPython_INCLUDE_DIR:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/include/python3.10 -DPython_FIND_REGISTRY:STRING=NEVER -DPython3_EXECUTABLE:PATH=/Users/tbridel/Documents/1-CODES/py3env/bin/python -DPython3_ROOT_DIR:PATH=/Users/tbridel/Documents/1-CODES/py3env -DPython3_INCLUDE_DIR:PATH=/opt/homebrew/opt/[email protected]/Frameworks/Python.framework/Versions/3.10/include/python3.10 -DPython3_FIND_REGISTRY:STRING=NEVER -DCMAKE_MAKE_PROGRAM:FILEPATH=/private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-build-env-tb3aezox/overlay/lib/python3.10/site-packages/ninja/data/bin/ninja -DSKBUILD=ON -DCMAKE_BUILD_TYPE:STRING=Release -DCMAKE_OSX_DEPLOYMENT_TARGET:STRING=12.0 -DCMAKE_OSX_ARCHITECTURES:STRING=arm64
        Source directory:
          /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57
        Working directory:
          /private/var/folders/8_/vb42yyg57lj90txl41w_m4ph0000gn/T/pip-install-s6p_6_f6/numbalsoda_aa6a60982d73407590bdf13f5c02eb57/_skbuild/macosx-12.0-arm64-3.10/cmake-build
      Please see CMake's output for more information.
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
  ERROR: Failed building wheel for numbalsoda
Failed to build numbalsoda
ERROR: Could not build wheels for numbalsoda, which is required to install pyproject.toml-based projects

🔴 This error also occurs with numbalsoda v0.3.3.

🟢 The pip install runs fine for numbalsoda v0.2.1 however.

Could you please advise @Nicholaswogan ?

Thank you very much in advance !
Best,
T. Bridel-Bertomeu

Is it possible to realize some stop condition?

It looks like yours solution is fastest in Python ODE tasks.

In solve_ivp() we have events which allow to implement stop condition during ODE process solving.
Is it possible to find/implement something similar in NumbaLSODA?

It could be useful if I dont know real parameter's range of integration but know where it should be stopped during entire integration process.

Thanks.

Adding other Runge-Kutta time integrators

Hye again !

Could you please share the steps one would have to take to add a new Runge-Kutta integrator to the project ?

I'm happy to implement it push a PR at some point, but I wanted to know if you had maybe a process that would make it easier to add a new Butcher-tableau based RK integrator ?

Thanks !
T. Bridel-Bertomeu

Jacobian matrix?

Hello, this package looks quite interesting & promising. I took a quick look at the examples and the wrapper. It doesn't look like NumbaLSODA is taking Jacobian matrix?

I did a few comparisons with scipy's solve_ivp, with BDF method and jacobian matrix provided. For small sized (say, with a few hundreds of ODEs), NumbaLSODA is indeed much faster but the performance of solve_ivp (BDF+jacobian) isn't totally unacceptable. For larger sized problem (say, with thousands of ODEs) NumbaLSODA is either marginally faster or slower - okay this isn't really fair, since in those cases with thousands of ODEs, with/without jacobian makes a huge difference in solve_ivp. NumbaLSODE would probably still outperform if solve_ivp is with LSODA and without Jacobian.

But I'd imagine jacobian may further improve the performance?

Thank you again.