ValueError second time .run() is called

[astrojuanlu]

This code:

import numpy as np
import nbkode

def pop(t, z, p):
    x, y = z
    a, b, c, d = p
    return np.array([a*x - b*x*y, -c*y + d*x*y])

z0 = np.array([10, 5], dtype=float)
t0 = 0
p = np.array([1.5, 1, 3, 1])

solver8 = nbkode.DOP853(pop, t0, z0, params = p)

ts3 = np.linspace(0, 15, 300)
solver8.run(ts3)
solver8.run(ts3)  # Again

Has this effect:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-14-065a3f1d5d78> in <module>
     15 ts3 = np.linspace(0, 15, 300)
     16 solver8.run(ts3)
---> 17 solver8.run(ts3)  # Again

~/.pyenv/versions/poliastro38/lib/python3.8/site-packages/nbkode/core.py in run(self, t)
    323 
    324         if t[0] < self._ts[0]:
--> 325             raise ValueError(
    326                 f"Cannot interpolate at t={t[0]} as it is smaller "
    327                 f"than the current smallest value in history ({self._ts[0]})"

ValueError: Cannot interpolate at t=0.0 as it is smaller than the current smallest value in history (14.742948523125582)

This is Python 3.8, nbkode 0.3.

[hgrecco]

Yes. This is not a bug but a design choice (which might be a bad one :-) but at least it was a conscious one). The API is more similar to a class based API than a function based API. The solver has a state defined by N (t, y). N is 2 by default, but it can be larger if the method requires it for calculation (i.e. multistep methods). These (t, y) are ALWAYS the ones PROVIDED by the algorithm (e.g. in fixed step methods tn = t0 + n dt )

When the user requires a t

• smaller than the minimum timepoint stored int the history (i.e. in the past), an exception is raised.
• between the minimum and maximum timepoints stored int the history the value for y is interpolated.
• larger than the minimum timepoint stored int the history (i.e. in the future), the algorith is stepped forward until the result can be provided.

Is this choice not ergonomic in your application?

[astrojuanlu]

Thanks for the explanation @hgrecco - I see that this is actually mentioned in the docs:

A solver instance remember and therefore the following command will not recalculate the what has happened between 0 and 10. [...]
You can do this as many times as you want, as long as you move forward in time.

Regarding your question:

Is this choice not ergonomic in your application?

It's not a big deal, I just had to RTFM :) I was thinking if I could propose an improvement to the API or the error message, but now that I understand what run does, it's quite evident.

Closing!

[summereasy]

I just started using nbkode and I encountered the same ValueError message which leads me to here :D

By reading the upper posts(and the documents), I understand that this is a design choice, not a bug. it is clear.

But my question is: Does that mean I have to compile the solver every-each time if I change the initial conditions y0 (or input params)?

Say, put a scenario like this:
I have got Sun, Earth and some 100,000 massless particles spread in the middle of space. I want to test their motion especially if their positions are in the neighborhood of Lagrangian Points. In this case, if the solver is a "jitted-function", then I can redo it with different y0s, even in parallel mode I think. In nbkode, do I need to compile the solver for each particle? (coz the compile time is relatively long ... ) Thanks!

Since I have very limited knowledge of numba, forgive me if this question is a "naive" one :)

BTW, nbkode is really great! it is easy to use, flexible as scipy does, and got some 200x runtime speedup in my test, and that is literally freaking me out! It really did help a lot. Thanks to @hgrecco and @maurosilber and all the contributors.

[maurosilber]

If you need to change the RHS function each time, yes (at least, until some of the issues linked in #5 are addressed in numba).

But if it is the same function, and only the initial conditions y0 change, it shouldn't recompile.

Note that you need to @numba.jit the function beforehand. While the solver jits the function for you inside __init__, it results in a "different" function, and it needs to recompile part of the solver. This jitting inside __init__ happens in two cases:

1. when it is not jitted
2. when it depends on parameters (rhs(t, y, p)), which can be passed with the params argument, as a closure (rhs(t, y)) is generated inside __init__.

Here is a minimal example:

import contextlib
import time

import numba
import nbkode

@contextlib.contextmanager
def timeit(message):
    start = time.perf_counter()
    yield
    end = time.perf_counter()
    total = end - start
    print(f"{total:.2e} seconds", message)

@numba.njit  # this is important
def func(t, y):
    return -y

for i in range(3):
    print("Loop", i)

    with timeit("Initialize"):
        solver = nbkode.RungeKutta23(func, 0, 1)

    with timeit("First run"):
        solver.run(5)

    with timeit("Second run"):
        solver.run(10)

which prints:

Loop 0
7.81e+00 seconds Initialize
../numbakit-ode/nbkode/core.py:617: NumbaExperimentalFeatureWarning: First-class function type feature is experimental
1.37e+01 seconds First run
1.06e-03 seconds Second run
Loop 1
1.05e-03 seconds Initialize
1.03e-03 seconds First run
1.10e-03 seconds Second run
Loop 2
1.04e-03 seconds Initialize
1.08e-03 seconds First run
1.05e-03 seconds Second run

[summereasy]

@maurosilber Thanks for the example code!
But with RHS functions coming with params, even those params are fixed values, the solver will be recompiled:

import contextlib
import time

import numba
import nbkode

@contextlib.contextmanager
def timeit(message):
    start = time.perf_counter()
    yield
    end = time.perf_counter()
    total = end - start
    print(f"{total:.2e} seconds", message)

@numba.njit  # this is important
def func(t, y, p):
    return -y * p

for i in range(3):
    print("Loop", i)

    with timeit("Initialize"):
        solver = nbkode.DOP853(func, 0, 1, params=2)

    with timeit("First run"):
        solver.run(5)

    with timeit("Second run"):
        solver.run(10)

It prints:

Loop 0
2.74e+00 seconds Initialize
.../nbkode/core.py:617: NumbaExperimentalFeatureWarning: First-class function type feature is experimental
8.38e+00 seconds First run
3.74e-04 seconds Second run
Loop 1
2.25e+00 seconds Initialize
.../nbkode/core.py:617: NumbaExperimentalFeatureWarning: First-class function type feature is experimental
8.37e+00 seconds First run
4.33e-04 seconds Second run
Loop 2
2.21e+00 seconds Initialize
.../nbkode/core.py:617: NumbaExperimentalFeatureWarning: First-class function type feature is experimental
8.53e+00 seconds First run
3.49e-04 seconds Second run

[maurosilber]

Yes, that's what I mentioned in item 2. But, you can do that closure outside the __init__.

For instance,

@numba.njit
def func_with_params(t, y, p):
   return -p * y

@numba.njit
def func(t, y):
    return func_with_params(t, y, 2)

Or, more generally:

def param_closure(func, params):
    @numba.njit
    def rhs(t, y):
        return func(t, y, params)
    return rhs

rhs = param_closure(func_with_params, 2)

[summereasy]

Oh! Right! It works! Finally, I got what you mean!
Great Thanks for the example, really counting on it!
Thanks again, nbkode rocks!