Merge pull request #33 from EconForge/albop/extrap

FIX: interp extrapolates with closest value on the grid.

Author: albop

interpolation/multilinear/fungen.py

@@ -20,6 +20,9 @@
 def clamp(x,a,b):
     return min(b,max(a,x))
 
+
+
+
 # returns the index of a 1d point along a 1d dimension
 @generated_jit(nopython=True)
 def get_index(gc, x):
@@ -195,3 +198,20 @@ def extract_row(a, n, tup):
     s = "def extract_row(a, n, tup): return ({},)".format(str.join(', ', [f"a[n,{i}]" for i in range(d)]))
     eval(compile(ast.parse(s),'<string>','exec'), dd)
     return dd['extract_row']
+
+
+
+# find closest point inside the grid domain
+@generated_jit
+def project(grid, point):
+    s = "def __project(grid, point):\n"
+    d = len(grid.types)
+    for i in range(d):
+        if isinstance(grid.types[i], numba.types.Array):
+            s += f"    x_{i} = min(max(point[{i}], grid[{i}][0]), grid[{i}][grid[{i}].shape[0]-1])\n"
+        else:
+            s += f"    x_{i} = min(max(point[{i}], grid[{i}][0]), grid[{i}][1])\n"
+    s += f"    return ({str.join(', ', ['x_{}'.format(i) for i in range(d)])},)"
+    d = {}
+    eval(compile(ast.parse(s),'<string>','exec'), d)
+    return d['__project']

interpolation/multilinear/mlinterp.py

@@ -18,7 +18,7 @@
 # Actual interpolation function #
 #################################
 
-from .fungen import fmap, funzip, get_coeffs, tensor_reduction, get_index, extract_row
+from .fungen import fmap, funzip, get_coeffs, tensor_reduction, get_index, extract_row, project
 
 from numba import njit
 from typing import Tuple
@@ -136,7 +136,8 @@ def {funname}(*args):
     grid = {grid_s}
     C = args[{it.d}]
     point = {point_s}
-    res = mlinterp(grid, C, point)
+    ppoint = project(grid, point)
+    res = mlinterp(grid, C, ppoint)
     return res
     """
         return source
@@ -152,7 +153,8 @@ def {funname}(*args):
     res = zeros(N)
     # return res
     for n in range(N):
-        res[n] = mlinterp(grid, C, {p_s})
+        ppoint = project(grid, {p_s})
+        res[n] = mlinterp(grid, C, ppoint)
     return res
 """
         return source
@@ -168,7 +170,8 @@ def {funname}(*args):
     N = points_x.shape[0]
     res = zeros(N)
     for n in range(N):
-        res[n] = mlinterp(grid, C, (points_x[n],))
+        ppoint = project(grid,(points_x[n],))
+        res[n] = mlinterp(grid, C, ppoint)
     return res
 """
         elif it.d==2:
@@ -184,7 +187,8 @@ def {funname}(*args):
     res = zeros((N,M))
     for n in range(N):
         for m in range(M):
-            res[n,m] = mlinterp(grid, C, (points_x[n], points_y[m]))
+            ppoint = project(grid,(points_x[n], points_y[m]))
+            res[n,m] = mlinterp(grid, C, ppoint)
     return res
 """
         else: