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 https://github.com/econ-ark/HARK/issues/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: 2x 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 https://github.com/econ-ark/HARK/pull/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 https://github.com/econ-ark/HARK/pull/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 https://github.com/econ-ark/HARK/pull/972.

Evidently, #972 https://github.com/econ-ark/HARK/pull/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 https://github.com/econ-ark/HARK/pull/1011.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/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.