Multiple scipy.integrate.ode instances

Question

I would like to use scipy.integrate.ode (or scipy.integrate.odeint) instances in multiple threads (one for each CPU core) in order to solve multiple IVPs at a time. However the documentation says: "This integrator is not re-entrant. You cannot have two ode instances using the “vode” integrator at the same time."

(Also odeint causes internal errors if instantiated multiple times although the documentation does not say so.)

Any idea what can be done?

Have you seen [this comment](http://stackoverflow.com/questions/23574146/solving-two-uncoupled-odes-within-a-loop-using-python-and-scipy-integrate-ode#comment36180664_23574361) and the corresponding posts? Also, odeint [might no be susceptible](http://comments.gmane.org/gmane.comp.python.scientific.user/36151) to the same problem. — Andras Deak -- Слава Україні, Dec 15 '15 at 16:29
@AndrasDeak: Thanks for looking at it! The link to the .pdf did not work anymore. However I rather need an implicit solver than an explicit one and Runge Kutta is an explicit one, I think. I also got internal errors when I tried to use mutithreaded odeint rather than ode. I think `ode.set_integrator('lsoda')` is the same implementation as odeint. — Matthias, Dec 16 '15 at 12:47

score 11 · Accepted Answer · answered Dec 15 '15 at 23:54

One option is to use multiprocessing (i.e. use processes instead of threads). Here's an example that uses the map function of the multiprocessing.Pool class.

The function solve takes a set of initial conditions and returns a solution generated by odeint. The "serial" version of the code in the main section calls solve repeatedly, once for each set of initial conditions in ics. The "multiprocessing" version uses the map function of a multiprocessing.Pool instance to run several processes simultaneously, each calling solve. The map function takes care of doling out the arguments to solve.

My computer has four cores, and as I increase num_processes, the speedup maxes out at about 3.6.

from __future__ import division, print_function

import sys
import time
import multiprocessing as mp
import numpy as np
from scipy.integrate import odeint



def lorenz(q, t, sigma, rho, beta):
    x, y, z = q
    return [sigma*(y - x), x*(rho - z) - y, x*y - beta*z]


def solve(ic):
    t = np.linspace(0, 200, 801)
    sigma = 10.0
    rho = 28.0
    beta = 8/3
    sol = odeint(lorenz, ic, t, args=(sigma, rho, beta), rtol=1e-10, atol=1e-12)
    return sol


if __name__ == "__main__":
    ics = np.random.randn(100, 3)

    print("multiprocessing:", end='')
    tstart = time.time()
    num_processes = 5
    p = mp.Pool(num_processes)
    mp_solutions = p.map(solve, ics)
    tend = time.time()
    tmp = tend - tstart
    print(" %8.3f seconds" % tmp)

    print("serial:         ", end='')
    sys.stdout.flush()
    tstart = time.time()
    serial_solutions = [solve(ic) for ic in ics]
    tend = time.time()
    tserial = tend - tstart
    print(" %8.3f seconds" % tserial)

    print("num_processes = %i, speedup = %.2f" % (num_processes, tserial/tmp))

    check = [(sol1 == sol2).all()
             for sol1, sol2 in zip(serial_solutions, mp_solutions)]
    if not all(check):
        print("There was at least one discrepancy in the solutions.")

On my computer, the output is:

multiprocessing:    6.904 seconds
serial:            24.756 seconds
num_processes = 5, speedup = 3.59

Do you know why the CPU load always corresponds with the number of processes (up to the number of cores) rather than with the set size of the initial conditions? — Matthias, Dec 16 '15 at 12:38
That sounds reasonable. After all, how could the CPU load exceed the number of cores? — Warren Weckesser, Dec 16 '15 at 14:36
E.g. if the set size of the initial conditions was 2 (i.e. there are two separate ode calculations) and the number of processes 4 (on a 4 core machine) the CPU load would always be 100 % although only 1 core could be used for each ode calculation. I would the CPU load have expected to remain at about 50 % and 2 processes to remain idle. — Matthias, Dec 16 '15 at 18:41
Each call to `solve` is run in a single process. With two sets of initial conditions, the best you can expect is to use two CPUs at 100%. That's what I see if make ics have shape (2, 3), and I set `num_processes=4`, but I have to make the `t` array large to clearly see it using the Activity Monitor on Mac OS X. In this case, multiprocessing is only useful if each solution takes a long time, otherwise the multiprocessing overhead time (e.g time to start the processes, time to copy the data between processes, etc) will dominate the overall time. — Warren Weckesser, Dec 16 '15 at 22:36
If you have just two sets of initial conditions, but you see four cores running at 100% for an extended period, then I don't know what is happening. — Warren Weckesser, Dec 16 '15 at 22:37

score 1 · Answer 2 · answered Oct 13 '18 at 18:17

SciPy.integrate.ode appears to use the LLNL SUNDIALS solvers, although SciPy doesn't say so explicitly, but they should, in my opinion.

The current version of the CVODE ode solver, 3.2.2, is re-entrant, which means that it can be used to solve multiple problems concurrently. The relevant information appears in User Documentation for CVODE v3.2.0 (SUNDIALS v3.2.0).

All state information used by cvode to solve a given problem is saved in a structure, and a pointer to that structure is returned to the user. There is no global data in the cvode package, and so, in this respect, it is reentrant. State information specific to the linear solver is saved in a separate structure, a pointer to which resides in the cvode memory structure. The reentrancy of cvode was motivated by the anticipated multicomputer extension, but is also essential in a uniprocessor setting where two or more problems are solved by intermixed calls to the package from within a single user program.

But I don't know whether SciPy.integrate.ode, or other ode solvers like scikits.odes.ode, support this concurrency.

Multiple scipy.integrate.ode instances

2 Answers2

Linked