Issue with speed and precision.

[HiggsBoson3310]

Hi,
I was trying to use the package and encountered some weird behavior. I was just trying to compare with Scipy integrators, and I have found very low speeds and very jagged solutions.

I can provide the script if required, but here is the call to the integrator
to=time.time()
solver = nbkode.RungeKutta45(rhs, t0=ts[0],y0=y0,params=[25.,0.5])
ti = time.time()
print("Time of initialization: %.8f s"%(ti-to))
ts, ys = solver.run(ts)
print(type(ys))
tf = time.time()
print("Numba time %.8f s, time of evaluation %.8f"%(tf-to, tf-ti))

The timing is:

Scipy time 0.02010775 s
Time of initialization: 4.81732988 s
C:\..\envs\nbode\lib\site-packages\nbkode\core.py:532: NumbaExperimentalFeatureWarning: First-class function type feature is experimental
if not step(t_bound, rhs, cache, *args): <class 'numpy.ndarray'>
Numba time 12.74941874 s, time of evaluation 7.93208885

and the output is very jagged

Is this an issue or am I calling the integrator incorrectly? I can provide used code if necessary. I am using Windows 10, with python version 3.8.12 and my numba versions are
numba 0.54.1 py38hf11a4ad_0 numbakit-ode 0.4 pypi_0 pypi

[hgrecco]

Numbakit-ode is tested against scipy integrators here so it should provide exactly the same result if the parameters are the same. Can you provide information a full script?

Ps. As a side comment, do not use time.time() for timing. If you do not want to use the timeit module, choose time.monotoic() or time.perf_counter().

[HiggsBoson3310]

Hi!
Thank you for the response.

Here is my script:

import numpy as np
import numba
import nbkode 
import scipy.integrate as integ
import matplotlib.pyplot as plt
import time
NS = 5
np.random.seed(5)
t0 = 0
y0 = np.random.random(NS)
Mi = 5*(np.random.random(size=(NS,NS))-0.5)
for i in range(NS):
    Mi[i,i] = -abs(Mi[i,i])
ts = np.linspace(0,50,5000)

def prof(t,to,gam):
    return  0.5*(1+np.tanh(gam*(t-to)))*np.exp(-gam**5 * ((t-to)**2)) 

def rhs(t,y,params):
    MM = Mi*0.5*(1+np.tanh(params[1]*(t-params[0])))*np.exp(-params[1]**5 * ((t-params[0])**2)) 
    return np.dot(MM,y)


to=time.time()
solver_sci = integ.solve_ivp(rhs, [ts[0],ts[-1]],y0,args=[[25.,0.5]], t_eval=ts)
tss = solver_sci.t
yss = solver_sci.y
tf = time.time()
print("Scipy time %.8f s"%(tf-to))

to=time.time()
solver = nbkode.RungeKutta45(rhs, t0=ts[0],y0=y0,params=[25.,0.5])
ti = time.time()
print("Time of initialization: %.8f s"%(ti-to))
ts, ys = solver.run(ts)
print(type(ys))
tf = time.time()
print("Numba time %.8f s, time of evaluation %.8f"%(tf-to, tf-ti))



fig, ax = plt.subplots(2,1)
ax[0].grid(True)
ax[0].set_title("Numba ode tools")
ax[1].grid(True)
ax[1].set_title("Scipy")
#plt.plot(ts, prof(ts, 25. , 0.5))
for i in range(NS):
    ax[0].plot(ts,ys[:,i],linestyle='-')
    ax[1].plot(tss,yss[i,:],linestyle="-")
#plt.plot(ts,ys[:,1],'r-')

plt.show()

and thanks for the tips on code timing!

[hgrecco]

With the current main I am getting this plot:

which looks the same.

[hgrecco]

Explaining the timing takes a little longer. Due to a problem with numba we cannot (yet) cache the compiled integrator and applied onto a generic function. This means that initialization takes a long time every single time.

Until this bug in numba is fixed, the performance improvement is visible only when running long simulations. But you can actually make it visible running a fixed number of steps:

to=time.time()
solver = nbkode.RungeKutta45(rhs, t0=ts[0], y0=y0,params=[25.,0.5])
tsn, ysn = solver.step(n=32)
tf = time.time()
print("Numba time %.8f s"%(tf-to))
print(f"{tsn[-1]}, {len(tsn)} / {len(ts)}")

solver = nbkode.RungeKutta45(rhs, t0=ts[0], y0=y0,params=[25.,0.5])
tsn, ysn = solver.step(n=1)
to=time.time()
tsn, ysn = solver.step(n=31)
tf = time.time()
print("Numba time %.8f s"%(tf-to))
print(f"{tsn[-1]}, {len(tsn)} / {len(ts)}")

which prints

Numba time 12.14614391 s
239.18818934278295, 32 / 5000
Numba time 0.00111103 s
239.18818934278295, 31 / 5000

[HiggsBoson3310]

Thanks for the explanation. I wonder if the "jaggedness" I'm seeing has to do with my installation. I'll try resinstalling and using a clean environment.

[hgrecco]

try installing directly from github to get the latest version

[HiggsBoson3310]

Update: I was not able to get the installation with windows working. It gave me this cryptic error

  Resolved https://github.com/hgrecco/numbakit-ode to commit 3b2ce017f7f614d57229ac97bf5f0fe330672f52
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
  Preparing metadata (pyproject.toml) ... error
  error: subprocess-exited-with-error

  × Preparing metadata (pyproject.toml) did not run successfully.
  │ exit code: 1
  ╰─> [45 lines of output]

with a bunch of lines.

I tried installing it in WSL and it worked like a charm. I can replicate the plot you obtained.

If you think the windows installation problem is important I could open another issue.

Thank you for your help!

[hgrecco]

Please do open an issue regarding windows installation.