Cubic Interpolation on lower bound

[alanlujan91]

This is related to issue #971.

An issue with CubicInterp is that it previously did not work as expected for the lower bound of the range where the interpolation is defined. As an example:

from HARK.interpolation import CubicInterp
import numpy as np
f = lambda x: x**2
df = lambda x: 2*x
x = np.linspace(-1,1,5)
y = f(x)
dy = df(x)
interp = CubicInterp(x,y,dy)
interp(x)

Out[1]: array([ nan, 0.25, 0.  , 0.25, 1.  ])

Clearly, the first element should be 1.. This is fixed when you set lower_extrap=True.

interp = CubicInterp(x,y,dy, lower_extrap=True)
interp(x)

Out[2]: array([1.  , 0.25, 0.  , 0.25, 1.  ])

However, x=-1 is not extrapolation, as it is one of the data points given, so it should return a value without setting lower_extrap=True.

A previous attempt to resolving this was the following line in `CubicInterp':

#972

...
y[x == self.x_list.min()] = self.y_list.min()
...

but this only works for a particular set of cases. A better solution is

#1008

...
y[x == self.x_list[0]] = self.y_list[0]
... 

interp = CubicInterp(x,y,dy)
interp(x)

Out[2]: array([1.  , 0.25, 0.  , 0.25, 1.  ])

which explicitly plugs in the lower bound. At some point, this correction was overwritten again, so HARK code is back to #972.

Evidently, #972 does not break any tests, as it creeped back in without notice. The way I realized this is back on HARK is because it created a conflict on #1011.

[mnwhite]

For background, this was a known bug / intended behavior when I wrote the original code. I intentionally didn't put the line `y[x == self.x_list[0]] = self.y_list[0]` because that's a really costly operation to fix a knife-edge case. It can probably be fixed by changing the "sided-ness" of the searchsorted operation that gets the segment (and then adjusting other parts of the algorithm).

…

On Thu, Aug 19, 2021 at 1:39 PM alanlujan91 ***@***.***> wrote: This is related to issue #971 <#971>. An issue with CubicInterp is that it previously did not work as expected for the lower bound of the range where the interpolation is defined. As an example: from HARK.interpolation import CubicInterp import numpy as np f = lambda x: x**2 df = lambda x: 2*x x = np.linspace(-1,1,5) y = f(x) dy = df(x) interp = CubicInterp(x,y,dy) interp(x) Out[1]: array([ nan, 0.25, 0. , 0.25, 1. ]) Clearly, the first element should be 1.. This is fixed when you set lower_extrap=True. interp = CubicInterp(x,y,dy, lower_extrap=True) interp(x) Out[2]: array([1. , 0.25, 0. , 0.25, 1. ]) However, x=-1 is not extrapolation, as it is one of the data points given, so it should return a value without setting lower_extrap=True. A previous attempt to resolving this was the following line in `CubicInterp': #972 <#972> ... y[x == self.x_list.min()] = self.y_list.min() ... but this only works for a particular set of cases. A better solution is #1008 <#1008> ... y[x == self.x_list[0]] = self.y_list[0] ... interp = CubicInterp(x,y,dy) interp(x) Out[2]: array([1. , 0.25, 0. , 0.25, 1. ]) which explicitly plugs in the lower bound. At some point, this correction was overwritten again, so HARK code is back to #972 <#972>. Evidently, #972 <#972> does not break any tests, as it creeped back in without notice. The way I realized this is back on HARK is because it created a conflict on #1011 <#1011>. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1060>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADKRAFMKFVGGC5PWMK5UAKTT5U6T3ANCNFSM5COWFZPQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .

[llorracc]

Well, it needs to be fixed one way or another (I wasted a fair amount of time running it down a few months ago).

My preference would be just to substituting the scipy version universally.