scipy 0.11.0 to 0.12.0 changes a linear scipy.interpolate.interp1d, breaks my constantly updated interpolator

Question

I have been playing around with a package that uses a linear scipy.interpolate.interp1d to create a history function for the ode solver in scipy, described here.

The relevant bit of code goes something like

def update(self, ti, Y):
        """ Add one new (ti, yi) to the interpolator """
        self.itpr.x = np.hstack([self.itpr.x,  [ti]])
        yi = np.array([Y]).T
        self.itpr.y = np.hstack([self.itpr.y,  yi])
        #self.itpr._y = np.hstack([self.itpr.y,  yi])
        self.itpr.fill_value = Y

Where "self.itpr" is initialized in __init__:

def __init__(self, g, tc=0):
    """ g(t) = expression of Y(t) for t<tc """

    self.g = g
    self.tc = tc
    # We must fill the interpolator with 2 points minimum
    self.itpr = scipy.interpolate.interp1d(
        np.array([tc-1, tc]),  # X
        np.array([self.g(tc), self.g(tc)]).T,  # Y
        kind='linear',  bounds_error=False, 
        fill_value = self.g(tc))

Where g is some function that returns an array of values that are solutions to a set of differential equations and tc is the current time.

This seems nice to me because a new interpolator object doesn't have to be re-created every time I want to update the ranges of values (which happens at each explicit time step during a simulation). This method of updating the interpolator works well under scipy v 0.11.0. However, after updating to v 0.12.0 I ran into issues. I see that the new interpolator now includes an array _y that seems to just be another copy of the original. Is it safe and/or sane to just update _y as outlined above as well? Is there a simpler, more pythonic way to address this that would hopefully be more robust to future updates in scipy? Again, in v 0.11 everything works well and expected results are produced, and in v 0.12 I get an IndexError when _y is referenced as it isn't updated in my function while y itself is.

Any help/pointers would be appreciated!

Answer 1

It looks like _y is just a copy of y that has been reshaped by interp1d._reshape_yi(). It should therefore be safe to just update it using:

  self.itpr._y = self.itpr._reshape_yi(self.itpr.y)

In fact, as far as I can tell it's only _y that gets used internally by the interpolator, so I think you could get away without actually updating y at all.

A much more elegant solution would be to make _y a property of the interpolator that returns a suitably reshaped copy of y. It's possible to achieve this by monkey-patching your specific instance of interp1d after it has been created (see Alex Martelli's answer here for more explanation):

x = np.arange(100)
y = np.random.randn(100)
itpr = interp1d(x,y)

# method to get self._y from self.y
def get_y(self):
    return self._reshape_yi(self.y)
meth = property(get_y,doc='reshaped version of self.y')

# make this a method of this interp1d instance only
basecls = type(itpr)
cls = type(basecls.__name__, (basecls,), {})
setattr(cls, '_y', meth)
itpr.__class__ = cls

# itpr._y is just a reshaped version of itpr.y
print itpr.y.shape,itpr._y.shape
>>> (100,) (100, 1)

Now itpr._y gets updated when you update itpr.y

itpr.x = np.arange(110)
itpr.y = np.random.randn(110)
print itpr._y.shape
>>> (110,) (110, 1)

This is all quite fiddly and not very Pythonic - it's much easier to fix the scipy source code (scipy/interpolate/interpolate.py). All you'd need to do is remove the last line from interp1d.__init__() where it sets:

self._y = y

and add these lines:

@property
def _y(self):
    return self._reshape_yi(self.y)