How to monitor the process of SciPy.odeint?

Question

SciPy can solve ode equations by scipy.integrate.odeint or other packages, but it gives result after the function has been solved completely. However, if the ode function is very complex, the program will take a lot of time(one or two days) to give the whole result. So how can I mointor the step it solve the equations(print out result when the equation hasn't been solved completely)?

Answer 1

When I was googling for an answer, I couldn't find a satisfactory one. So I made a simple gist with a proof-of-concept solution using the tqdm project. Hope that helps you.

Edit: Moderators asked me to give an explanation of what is going on in the link above.

First of all, I am using scipy's OOP version of odeint (solve_ivp) but you could adapt it back to odeint. Say you want to integrate from time T0 to T1 and you want to show progress for every 0.1% of progress. You can modify your ode function to take two extra parameters, a pbar (progress bar) and a state (current state of integration). Like so:

def fun(t, y, omega, pbar, state):
    # state is a list containing last updated time t:
    # state = [last_t, dt]
    # I used a list because its values can be carried between function
    # calls throughout the ODE integration
    last_t, dt = state
    
    # let's subdivide t_span into 1000 parts
    # call update(n) here where n = (t - last_t) / dt
    time.sleep(0.1)
    n = int((t - last_t)/dt)
    pbar.update(n)
    
    # we need this to take into account that n is a rounded number.
    state[0] = last_t + dt * n
    
    # YOUR CODE HERE
    dydt = 1j * y * omega
    return dydt

This is necessary because the function itself must know where it is located, but scipy's odeint doesn't really give this context to the function. Then, you can integrate fun with the following code:

T0 = 0
T1 = 1
t_span = (T0, T1)
omega = 20
y0 = np.array([1], dtype=np.complex)
t_eval = np.arange(*t_span, 0.25/omega)

with tqdm(total=1000, unit="‰") as pbar:
    sol = solve_ivp(
        fun,
        t_span,
        y0,
        t_eval=t_eval,
        args=[omega, pbar, [T0, (T1-T0)/1000]],
    )

Note that anything mutable (like a list) in the args is instantiated once and can be changed from within the function. I recommend doing this rather than using a global variable.

This will show a progress bar that looks like this:

100%|█████████▉| 999/1000 [00:13<00:00, 71.69‰/s]

Answer 2

You could split the integration domain and integrate the segments, taking the last value of the previous as initial condition of the next segment. In-between, print out whatever you want. Use numpy.concatenate to assemble the pieces if necessary.

In a standard example of a 3-body solar system simulation, replacing the code

u0 = solsys.getState0();
t = np.arange(0, 100*365.242*day, 0.5*day);
%timeit u_res = odeint(lambda u,t: solsys.getDerivs(u), u0, t, atol = 1e11*1e-8, rtol = 1e-9)

output: 1 loop, best of 3: 5.53 s per loop

with a progress-reporting code

def progressive(t,N):
    nk = [ int(n+0.5) for n in np.linspace(0,len(t),N+1) ]
    u0 = solsys.getState0();
    u_seg = [ np.array([u0]) ];
    for k in range(N):
        u_seg.append( odeint(lambda u,t: solsys.getDerivs(u), u0, t[nk[k]:nk[k+1]], atol = 1e11*1e-8, rtol = 1e-9)[1:] )
        print t[nk[k]]/day
        for b in solsys.bodies: print("%10s %s"%(b.name,b.x))
    return np.concatenate(u_seg)
%timeit u_res = progressive(t,20)

output: 1 loop, best of 3: 5.96 s per loop

shows only a slight 8% overhead for console printing. With a more substantive ODE function, the fraction of the reporting overhead will reduce significantly.

---

That said, python, at least with its standard packages, is not the tool for industrial-scale number-crunching. Always use compiled versions with strong typing of variables to reduce interpretative overhead as much as possible.

• Use some heavily developed and tested package like Sundials or the julia-lang framework differentialequations.jl directly coding the ODE function in an appropriate compiled language. Use the higher-order methods for larger step sizes, thus smaller steps. Test if using implicit or exponential/Rosenbrock methods reduces the number of steps or ODE function evaluations per fixed interval further. The difference can be a factor of 10 to 100 in speedup.

• Use a python wrapper of the above with some acceleration-friendly implementation of your ODE function.

• Use the quasi-source-translating tool JITcode to translate your python ODE function to a spaghetti list of C instruction that then give a compiled function that can be (almost) directly called from the compiled FORTRAN kernel of odeint.

Answer 3

Simple and Clear. If you want to integrate an ODE from T0 to T1: In the last line of the code, before return, you can use print((t/T1)*100,end='') Then use a sys.stdout.flush() to keep the same line of printing.

Here is an example. My integrating time [0 0.2]

ddt[-2]=(beta/(Ap2*(L-x)))*(-Qgap+Ap*u)
ddt[-1]=(beta/(Ap2*(L+x)))*(Qgap-Ap*u)
print("\rCompletion percentage "+str(format(((t/0.2)*100),".4f")),end='')
sys.stdout.flush()
return ddt

It slows a bit the solving process by fraction of seconds, but it serves perfectly the purpose rather than creating new functions.