From e1a4356b1239169d567807adb05a95c3a9e56040 Mon Sep 17 00:00:00 2001 From: "Gary M. Stump" Date: Tue, 19 Apr 2022 09:13:58 -0400 Subject: [PATCH 1/4] added econforge interpolation into LinearInterp and scipy interpolation into CubicInterp --- HARK/fast_interpolation.py | 16 ++++++++++++++++ HARK/interpolation.py | 16 ++++++++++++++++ 2 files changed, 32 insertions(+) create mode 100644 HARK/fast_interpolation.py diff --git a/HARK/fast_interpolation.py b/HARK/fast_interpolation.py new file mode 100644 index 000000000..a98f8c0c0 --- /dev/null +++ b/HARK/fast_interpolation.py @@ -0,0 +1,16 @@ +import numpy as np +from interpolation import interp # econ forge +from scipy.interpolate import CubicSpline # scipy + +# econforge interp and then uses the hark extrapolation for values outside of the range +def fast_linear_1d_interp(x, obj): + result_temp = interp(obj.x_list, obj.y_list, x) + above_upper_bound = x > obj.x_list[-1] + i = len(obj.x_list) - 1 + alpha = (x[above_upper_bound] - obj.x_list[i - 1]) / (obj.x_list[i] - obj.x_list[i - 1]) + result_temp[above_upper_bound] = (1.0 - alpha) * obj.y_list[i - 1] + alpha * obj.y_list[i] + return [result_temp] + +# gets the scipy cubic spline object +def get_cubic_spline(obj): + return CubicSpline(np.array(obj.x_list), np.array(obj.y_list)) diff --git a/HARK/interpolation.py b/HARK/interpolation.py index fc18d7931..94f89f463 100644 --- a/HARK/interpolation.py +++ b/HARK/interpolation.py @@ -15,6 +15,8 @@ from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP from .core import MetricObject +from HARK.fast_interpolation import fast_linear_1d_interp, get_cubic_spline + def _isscalar(x): """ @@ -814,6 +816,10 @@ def _evalOrDer(self, x, _eval, _Der): A list including the level and/or derivative of the interpolated function where requested. """ + # econforge interpolation check + if _eval and not _Der and not self.decay_extrap and x[0] >= self.x_list[0]: + return fast_linear_1d_interp(x, self) + i = np.maximum(np.searchsorted(self.x_list[:-1], x), 1) alpha = (x - self.x_list[i - 1]) / (self.x_list[i] - self.x_list[i - 1]) @@ -989,6 +995,11 @@ def __init__( self.coeffs.append(temp) self.coeffs = np.array(self.coeffs) + # cache scipy cubic spline if all x values are increasing + self.cs = None + if not np.any(np.diff(self.x_list) <= 0): + self.cs = get_cubic_spline(self) + def _evaluate(self, x): """ Returns the level of the interpolated function at each value in x. Only @@ -1014,6 +1025,11 @@ def _evaluate(self, x): - self.coeffs[pos, 2] * np.exp(alpha * self.coeffs[pos, 3]) ) else: + + # use cubic scipy spline if cached + if self.cs is not None: + return self.cs(x) + m = len(x) pos = np.searchsorted(self.x_list, x) y = np.zeros(m) From f33d8ccabb60df033c1c3616de283cd637ed85ff Mon Sep 17 00:00:00 2001 From: "Gary M. Stump" Date: Mon, 25 Apr 2022 12:34:33 -0400 Subject: [PATCH 2/4] updated linear interpolation to catch empty array and temporarily removed cubic spline interpolation until tolerance issues are resolved --- HARK/interpolation.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/HARK/interpolation.py b/HARK/interpolation.py index 94f89f463..7647c0303 100644 --- a/HARK/interpolation.py +++ b/HARK/interpolation.py @@ -816,9 +816,14 @@ def _evalOrDer(self, x, _eval, _Der): A list including the level and/or derivative of the interpolated function where requested. """ - # econforge interpolation check - if _eval and not _Der and not self.decay_extrap and x[0] >= self.x_list[0]: - return fast_linear_1d_interp(x, self) + # ========== fast interpolation check for econ hark library ========== + # sometimes x is empty causing the below to error, so added a try except + # tested with an extra check to see if x is empty, but this added to the run time + try: + if _eval and not _Der and not self.decay_extrap and x[0] >= self.x_list[0]: + return fast_linear_1d_interp(x, self) + except: + print("empty x sent to LinearInterp") i = np.maximum(np.searchsorted(self.x_list[:-1], x), 1) alpha = (x - self.x_list[i - 1]) / (self.x_list[i] - self.x_list[i - 1]) @@ -995,11 +1000,6 @@ def __init__( self.coeffs.append(temp) self.coeffs = np.array(self.coeffs) - # cache scipy cubic spline if all x values are increasing - self.cs = None - if not np.any(np.diff(self.x_list) <= 0): - self.cs = get_cubic_spline(self) - def _evaluate(self, x): """ Returns the level of the interpolated function at each value in x. Only @@ -1026,10 +1026,6 @@ def _evaluate(self, x): ) else: - # use cubic scipy spline if cached - if self.cs is not None: - return self.cs(x) - m = len(x) pos = np.searchsorted(self.x_list, x) y = np.zeros(m) From 289ff4f7f89780064fef1208b83b7638586a5566 Mon Sep 17 00:00:00 2001 From: "Gary M. Stump" Date: Tue, 26 Apr 2022 09:28:47 -0400 Subject: [PATCH 3/4] inserted econforge interpolation to bilinear and trilinear hark interpolation methods --- HARK/interpolation.py | 165 ++++++++++++++++++++++-------------------- 1 file changed, 88 insertions(+), 77 deletions(-) diff --git a/HARK/interpolation.py b/HARK/interpolation.py index 7647c0303..e997fd55e 100644 --- a/HARK/interpolation.py +++ b/HARK/interpolation.py @@ -16,6 +16,7 @@ from .core import MetricObject from HARK.fast_interpolation import fast_linear_1d_interp, get_cubic_spline +from interpolation.splines import CGrid, eval_linear def _isscalar(x): @@ -1386,34 +1387,39 @@ def __init__(self, f_values, x_list, y_list, xSearchFunc=None, ySearchFunc=None) self.xSearchFunc = xSearchFunc self.ySearchFunc = ySearchFunc + self.customgrid = CGrid(self.x_list, self.y_list) + def _evaluate(self, x, y): - """ - Returns the level of the interpolated function at each value in x,y. - Only called internally by HARKinterpolator2D.__call__ (etc). - """ - if _isscalar(x): - x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) - y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) - else: - x_pos = self.xSearchFunc(self.x_list, x) - x_pos[x_pos < 1] = 1 - x_pos[x_pos > self.x_n - 1] = self.x_n - 1 - y_pos = self.ySearchFunc(self.y_list, y) - y_pos[y_pos < 1] = 1 - y_pos[y_pos > self.y_n - 1] = self.y_n - 1 - alpha = (x - self.x_list[x_pos - 1]) / ( - self.x_list[x_pos] - self.x_list[x_pos - 1] - ) - beta = (y - self.y_list[y_pos - 1]) / ( - self.y_list[y_pos] - self.y_list[y_pos - 1] - ) - f = ( - (1 - alpha) * (1 - beta) * self.f_values[x_pos - 1, y_pos - 1] - + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos] - + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1] - + alpha * beta * self.f_values[x_pos, y_pos] - ) - return f + + return eval_linear(self.customgrid, self.f_values, np.column_stack((x,y))) + + # """ + # Returns the level of the interpolated function at each value in x,y. + # Only called internally by HARKinterpolator2D.__call__ (etc). + # """ + # if _isscalar(x): + # x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) + # y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) + # else: + # x_pos = self.xSearchFunc(self.x_list, x) + # x_pos[x_pos < 1] = 1 + # x_pos[x_pos > self.x_n - 1] = self.x_n - 1 + # y_pos = self.ySearchFunc(self.y_list, y) + # y_pos[y_pos < 1] = 1 + # y_pos[y_pos > self.y_n - 1] = self.y_n - 1 + # alpha = (x - self.x_list[x_pos - 1]) / ( + # self.x_list[x_pos] - self.x_list[x_pos - 1] + # ) + # beta = (y - self.y_list[y_pos - 1]) / ( + # self.y_list[y_pos] - self.y_list[y_pos - 1] + # ) + # f = ( + # (1 - alpha) * (1 - beta) * self.f_values[x_pos - 1, y_pos - 1] + # + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos] + # + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1] + # + alpha * beta * self.f_values[x_pos, y_pos] + # ) + # return f def _derX(self, x, y): """ @@ -1544,57 +1550,62 @@ def __init__( self.ySearchFunc = ySearchFunc self.zSearchFunc = zSearchFunc + self.customgrid = CGrid(self.x_list, self.y_list, self.z_list) + def _evaluate(self, x, y, z): - """ - Returns the level of the interpolated function at each value in x,y,z. - Only called internally by HARKinterpolator3D.__call__ (etc). - """ - if _isscalar(x): - x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) - y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) - z_pos = max(min(self.zSearchFunc(self.z_list, z), self.z_n - 1), 1) - else: - x_pos = self.xSearchFunc(self.x_list, x) - x_pos[x_pos < 1] = 1 - x_pos[x_pos > self.x_n - 1] = self.x_n - 1 - y_pos = self.ySearchFunc(self.y_list, y) - y_pos[y_pos < 1] = 1 - y_pos[y_pos > self.y_n - 1] = self.y_n - 1 - z_pos = self.zSearchFunc(self.z_list, z) - z_pos[z_pos < 1] = 1 - z_pos[z_pos > self.z_n - 1] = self.z_n - 1 - alpha = (x - self.x_list[x_pos - 1]) / ( - self.x_list[x_pos] - self.x_list[x_pos - 1] - ) - beta = (y - self.y_list[y_pos - 1]) / ( - self.y_list[y_pos] - self.y_list[y_pos - 1] - ) - gamma = (z - self.z_list[z_pos - 1]) / ( - self.z_list[z_pos] - self.z_list[z_pos - 1] - ) - f = ( - (1 - alpha) - * (1 - beta) - * (1 - gamma) - * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1] - + (1 - alpha) - * (1 - beta) - * gamma - * self.f_values[x_pos - 1, y_pos - 1, z_pos] - + (1 - alpha) - * beta - * (1 - gamma) - * self.f_values[x_pos - 1, y_pos, z_pos - 1] - + (1 - alpha) * beta * gamma * self.f_values[x_pos - 1, y_pos, z_pos] - + alpha - * (1 - beta) - * (1 - gamma) - * self.f_values[x_pos, y_pos - 1, z_pos - 1] - + alpha * (1 - beta) * gamma * self.f_values[x_pos, y_pos - 1, z_pos] - + alpha * beta * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1] - + alpha * beta * gamma * self.f_values[x_pos, y_pos, z_pos] - ) - return f + + return eval_linear(self.customgrid, self.f_values, np.column_stack((x,y,z))) + + # """ + # Returns the level of the interpolated function at each value in x,y,z. + # Only called internally by HARKinterpolator3D.__call__ (etc). + # """ + # if _isscalar(x): + # x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) + # y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) + # z_pos = max(min(self.zSearchFunc(self.z_list, z), self.z_n - 1), 1) + # else: + # x_pos = self.xSearchFunc(self.x_list, x) + # x_pos[x_pos < 1] = 1 + # x_pos[x_pos > self.x_n - 1] = self.x_n - 1 + # y_pos = self.ySearchFunc(self.y_list, y) + # y_pos[y_pos < 1] = 1 + # y_pos[y_pos > self.y_n - 1] = self.y_n - 1 + # z_pos = self.zSearchFunc(self.z_list, z) + # z_pos[z_pos < 1] = 1 + # z_pos[z_pos > self.z_n - 1] = self.z_n - 1 + # alpha = (x - self.x_list[x_pos - 1]) / ( + # self.x_list[x_pos] - self.x_list[x_pos - 1] + # ) + # beta = (y - self.y_list[y_pos - 1]) / ( + # self.y_list[y_pos] - self.y_list[y_pos - 1] + # ) + # gamma = (z - self.z_list[z_pos - 1]) / ( + # self.z_list[z_pos] - self.z_list[z_pos - 1] + # ) + # f = ( + # (1 - alpha) + # * (1 - beta) + # * (1 - gamma) + # * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1] + # + (1 - alpha) + # * (1 - beta) + # * gamma + # * self.f_values[x_pos - 1, y_pos - 1, z_pos] + # + (1 - alpha) + # * beta + # * (1 - gamma) + # * self.f_values[x_pos - 1, y_pos, z_pos - 1] + # + (1 - alpha) * beta * gamma * self.f_values[x_pos - 1, y_pos, z_pos] + # + alpha + # * (1 - beta) + # * (1 - gamma) + # * self.f_values[x_pos, y_pos - 1, z_pos - 1] + # + alpha * (1 - beta) * gamma * self.f_values[x_pos, y_pos - 1, z_pos] + # + alpha * beta * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1] + # + alpha * beta * gamma * self.f_values[x_pos, y_pos, z_pos] + # ) + # return f def _derX(self, x, y, z): """ From 679c43781d48489ddc83d7fccb5bfc72633b9e9b Mon Sep 17 00:00:00 2001 From: "Gary M. Stump" Date: Wed, 27 Apr 2022 09:22:34 -0400 Subject: [PATCH 4/4] added try except around Bilinear and Trilinear interpolation calls to account for an outlier case where x and y test values are one point with integers which causes an error in the expected data format for python 3.8 and up --- HARK/interpolation.py | 168 ++++++++++++++++++++++-------------------- 1 file changed, 87 insertions(+), 81 deletions(-) diff --git a/HARK/interpolation.py b/HARK/interpolation.py index e997fd55e..e636fabac 100644 --- a/HARK/interpolation.py +++ b/HARK/interpolation.py @@ -1391,35 +1391,38 @@ def __init__(self, f_values, x_list, y_list, xSearchFunc=None, ySearchFunc=None) def _evaluate(self, x, y): - return eval_linear(self.customgrid, self.f_values, np.column_stack((x,y))) - - # """ - # Returns the level of the interpolated function at each value in x,y. - # Only called internally by HARKinterpolator2D.__call__ (etc). - # """ - # if _isscalar(x): - # x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) - # y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) - # else: - # x_pos = self.xSearchFunc(self.x_list, x) - # x_pos[x_pos < 1] = 1 - # x_pos[x_pos > self.x_n - 1] = self.x_n - 1 - # y_pos = self.ySearchFunc(self.y_list, y) - # y_pos[y_pos < 1] = 1 - # y_pos[y_pos > self.y_n - 1] = self.y_n - 1 - # alpha = (x - self.x_list[x_pos - 1]) / ( - # self.x_list[x_pos] - self.x_list[x_pos - 1] - # ) - # beta = (y - self.y_list[y_pos - 1]) / ( - # self.y_list[y_pos] - self.y_list[y_pos - 1] - # ) - # f = ( - # (1 - alpha) * (1 - beta) * self.f_values[x_pos - 1, y_pos - 1] - # + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos] - # + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1] - # + alpha * beta * self.f_values[x_pos, y_pos] - # ) - # return f + try: + return eval_linear(self.customgrid, self.f_values, np.column_stack((x,y))) + except: + print("econforge point format not supported : ",np.column_stack((x,y))) + + """ + Returns the level of the interpolated function at each value in x,y. + Only called internally by HARKinterpolator2D.__call__ (etc). + """ + if _isscalar(x): + x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) + y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) + else: + x_pos = self.xSearchFunc(self.x_list, x) + x_pos[x_pos < 1] = 1 + x_pos[x_pos > self.x_n - 1] = self.x_n - 1 + y_pos = self.ySearchFunc(self.y_list, y) + y_pos[y_pos < 1] = 1 + y_pos[y_pos > self.y_n - 1] = self.y_n - 1 + alpha = (x - self.x_list[x_pos - 1]) / ( + self.x_list[x_pos] - self.x_list[x_pos - 1] + ) + beta = (y - self.y_list[y_pos - 1]) / ( + self.y_list[y_pos] - self.y_list[y_pos - 1] + ) + f = ( + (1 - alpha) * (1 - beta) * self.f_values[x_pos - 1, y_pos - 1] + + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos] + + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1] + + alpha * beta * self.f_values[x_pos, y_pos] + ) + return f def _derX(self, x, y): """ @@ -1554,58 +1557,61 @@ def __init__( def _evaluate(self, x, y, z): - return eval_linear(self.customgrid, self.f_values, np.column_stack((x,y,z))) - - # """ - # Returns the level of the interpolated function at each value in x,y,z. - # Only called internally by HARKinterpolator3D.__call__ (etc). - # """ - # if _isscalar(x): - # x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) - # y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) - # z_pos = max(min(self.zSearchFunc(self.z_list, z), self.z_n - 1), 1) - # else: - # x_pos = self.xSearchFunc(self.x_list, x) - # x_pos[x_pos < 1] = 1 - # x_pos[x_pos > self.x_n - 1] = self.x_n - 1 - # y_pos = self.ySearchFunc(self.y_list, y) - # y_pos[y_pos < 1] = 1 - # y_pos[y_pos > self.y_n - 1] = self.y_n - 1 - # z_pos = self.zSearchFunc(self.z_list, z) - # z_pos[z_pos < 1] = 1 - # z_pos[z_pos > self.z_n - 1] = self.z_n - 1 - # alpha = (x - self.x_list[x_pos - 1]) / ( - # self.x_list[x_pos] - self.x_list[x_pos - 1] - # ) - # beta = (y - self.y_list[y_pos - 1]) / ( - # self.y_list[y_pos] - self.y_list[y_pos - 1] - # ) - # gamma = (z - self.z_list[z_pos - 1]) / ( - # self.z_list[z_pos] - self.z_list[z_pos - 1] - # ) - # f = ( - # (1 - alpha) - # * (1 - beta) - # * (1 - gamma) - # * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1] - # + (1 - alpha) - # * (1 - beta) - # * gamma - # * self.f_values[x_pos - 1, y_pos - 1, z_pos] - # + (1 - alpha) - # * beta - # * (1 - gamma) - # * self.f_values[x_pos - 1, y_pos, z_pos - 1] - # + (1 - alpha) * beta * gamma * self.f_values[x_pos - 1, y_pos, z_pos] - # + alpha - # * (1 - beta) - # * (1 - gamma) - # * self.f_values[x_pos, y_pos - 1, z_pos - 1] - # + alpha * (1 - beta) * gamma * self.f_values[x_pos, y_pos - 1, z_pos] - # + alpha * beta * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1] - # + alpha * beta * gamma * self.f_values[x_pos, y_pos, z_pos] - # ) - # return f + try: + return eval_linear(self.customgrid, self.f_values, np.column_stack((x,y,z))) + except: + print("econforge point format not supported : ",np.column_stack((x,y,z))) + + """ + Returns the level of the interpolated function at each value in x,y,z. + Only called internally by HARKinterpolator3D.__call__ (etc). + """ + if _isscalar(x): + x_pos = max(min(self.xSearchFunc(self.x_list, x), self.x_n - 1), 1) + y_pos = max(min(self.ySearchFunc(self.y_list, y), self.y_n - 1), 1) + z_pos = max(min(self.zSearchFunc(self.z_list, z), self.z_n - 1), 1) + else: + x_pos = self.xSearchFunc(self.x_list, x) + x_pos[x_pos < 1] = 1 + x_pos[x_pos > self.x_n - 1] = self.x_n - 1 + y_pos = self.ySearchFunc(self.y_list, y) + y_pos[y_pos < 1] = 1 + y_pos[y_pos > self.y_n - 1] = self.y_n - 1 + z_pos = self.zSearchFunc(self.z_list, z) + z_pos[z_pos < 1] = 1 + z_pos[z_pos > self.z_n - 1] = self.z_n - 1 + alpha = (x - self.x_list[x_pos - 1]) / ( + self.x_list[x_pos] - self.x_list[x_pos - 1] + ) + beta = (y - self.y_list[y_pos - 1]) / ( + self.y_list[y_pos] - self.y_list[y_pos - 1] + ) + gamma = (z - self.z_list[z_pos - 1]) / ( + self.z_list[z_pos] - self.z_list[z_pos - 1] + ) + f = ( + (1 - alpha) + * (1 - beta) + * (1 - gamma) + * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1] + + (1 - alpha) + * (1 - beta) + * gamma + * self.f_values[x_pos - 1, y_pos - 1, z_pos] + + (1 - alpha) + * beta + * (1 - gamma) + * self.f_values[x_pos - 1, y_pos, z_pos - 1] + + (1 - alpha) * beta * gamma * self.f_values[x_pos - 1, y_pos, z_pos] + + alpha + * (1 - beta) + * (1 - gamma) + * self.f_values[x_pos, y_pos - 1, z_pos - 1] + + alpha * (1 - beta) * gamma * self.f_values[x_pos, y_pos - 1, z_pos] + + alpha * beta * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1] + + alpha * beta * gamma * self.f_values[x_pos, y_pos, z_pos] + ) + return f def _derX(self, x, y, z): """