Fast table interpolation on regular time data for ODE

Question

I am using scipy.integrate.odeint to simulate the reaction of a system with known input signals via integration. The simplified code below illustrates what I'm doing. It simulates the response of a number of first order systems to sinusoidal inputs by first generating the data and interpolating to the data in the integration function.

I would like to speed up the simulations by exploiting the fact that the time points for the data table are regularly spaced and that we are doing interpolation between the same rows of the table for all the input variables. Profiling has shown that the actual code I'm working on is indeed spending a significant amount of time in numpy.interp.

The exact question then is if I should just roll my own interpolation function which handles this special case or if there exists an off-the shelf solution which handles this case efficiently. I have searched for this but come up empty handed.

import numpy
import scipy.integrate

import matplotlib.pyplot as plt

Ntime = 400
Ninputs = 4
time = numpy.arange(Ntime)

# Generate some out-of-phase sinusoids for input data
DATA = numpy.array([numpy.sin(time/10 - i) for i in range(Ninputs)]).T

# Time constants
taus = numpy.ones(Ninputs)*10

def intfun(x, t):
    # I would like to speed up this lookup:
    inputs = numpy.array([numpy.interp(t, time, DATA[:, i]) 
                          for i in range(Ninputs)])

    return -(1/taus)*x + inputs  # first order

x0 = numpy.ones(Ninputs)

x = scipy.integrate.odeint(intfun, x0, time)

plt.plot(time, x)
plt.show()

Answer 1

I see the following inefficiencies in your approach:

• For every call of intfun, numpy.interp has to start from zero searching the point from DATA used for interpolation every time it is called. Though this uses a compiled binary search, it takes some time.

• For every call of intfun, you have a considerable overhead from numpy.interp which was made for interpolating more than one data point.

• You do not exploit the data being equidistant at all.

• You have to loop over the input dimension (Ninputs), as numpy.interp only works with one-dimensional data.

So, what can you do?

1. To avoid the search, you can exploit the regular sampling and only use the relevant part of the array for interpolation:

   def intfun(x,t):
       j = int(t)
       inputs = numpy.array([
                   numpy.interp( t, time[j:j+1], DATA[j:j+1,i] )
                   for i in range(Ninputs)
                   ])
       return -(1/taus)*x + inputs
   

   Empirically, this gives a slight speed boost.

2. You can exploit your knowledge of the regular sampling and do the interpolation yourself. This way you get rid of the overhead from numpy.interp and you can do a multidimensional interpolation:

   def intfun(x,t):
       j = int(t)
       t -= j
       inputs = DATA[j]*(1-t) + DATA[j+1]*t
       return -(1/taus)*x + inputs
   

   Empirically, this gives a better speed boost than the above.

3. You use a fixed step size for the integration (instead of an adaptive integrator like odeint). This allows you to perform the interpolation beforehand in a vectorised (and hence efficient) manner. However, you may need more steps and lose time this way. Also you have to go through the trouble of finding a good step size and you have to find or implement a fixed-step-size integrator for Python.

4. You switch to a compiled language with static types. This way, you get rid of the Python overhead from the interpolation operations in suggestion 2 and you can make use of an adaptive step size at the same time.