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()