ENH: added mlinterp

.travis.yml

@@ -1,8 +1,7 @@
 language: python
 
 python:
-    - "2.7"
-    - "3.5"
+    - "3.6"
 
 sudo: false
 
@@ -11,11 +10,7 @@ install:
   # conda line below will keep everything up-to-date.  We do this
   # conditionally because it saves us some downloading if the version is
   # the same.
-  - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
-      wget http://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh;
-    else
-      wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
-    fi
+  - wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
   - bash miniconda.sh -b -p $HOME/miniconda
   - export PATH="$HOME/miniconda/bin:$PATH"
   - hash -r

interpolation/__init__.py

@@ -1 +1 @@
-from .multilinear.mlinterp import interp
+from .multilinear.mlinterp import interp, mlinterp

interpolation/multilinear/fungen.py

@@ -169,3 +169,11 @@ def fun(C,l): return (1-l[0])*((1-l[1])*C[0][0]+l[1]*C[0][1]) + l[0]*((1-l[1])*C
     else:
         print("Not implemented")
     return fun# funzip(((1,2), (2,3), (4,3)))n
+
+@generated_jit
+def extract_row(a, n, tup):
+    d = len(tup.types)
+    dd = {}
+    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']

interpolation/multilinear/mlinterp.py

@@ -18,26 +18,45 @@
 # Actual interpolation function #
 #################################
 
-from .fungen import fmap, funzip, get_coeffs, tensor_reduction, get_index
+from .fungen import fmap, funzip, get_coeffs, tensor_reduction, get_index, extract_row
 
 from numba import njit
 from typing import Tuple
 
 from numba.types import UniTuple, Array, float64
 from numba.types import Float
 
+import numpy as np
+from numba import generated_jit
 
 #logic of multilinear interpolation
 
-@njit
-def _mlinterp(grid: Tuple, c, u: Tuple)->float:
 
-    # get indices and barycentric coordinates
-    tmp = fmap(get_index, grid, u)
-    indices, barycenters = funzip(tmp)
-    coeffs = get_coeffs(c, indices)
-    v = tensor_reduction(coeffs, barycenters)
-    return v
+@generated_jit
+def mlinterp(grid, c, u):
+    if isinstance(u, UniTuple):
+        def mlininterp(grid: Tuple, c: Array, u: Tuple)->float:
+            # get indices and barycentric coordinates
+            tmp = fmap(get_index, grid, u)
+            indices, barycenters = funzip(tmp)
+            coeffs = get_coeffs(c, indices)
+            v = tensor_reduction(coeffs, barycenters)
+            return v
+    elif isinstance(u, Array) and u.ndim==2:
+        def mlininterp(grid: Tuple, c: Array, u: Array)->float:
+            N = u.shape[0]
+            res = np.zeros(N)
+            for n in range(N):
+                uu = extract_row(u, n, grid)
+                # get indices and barycentric coordinates
+                tmp = fmap(get_index, grid, uu)
+                indices, barycenters = funzip(tmp)
+                coeffs = get_coeffs(c, indices)
+                res[n] = tensor_reduction(coeffs, barycenters)
+            return res
+    else:
+        mlininterp = None
+    return mlininterp
 
 ### The rest of this file constrcts function `interp`
 
@@ -90,6 +109,7 @@ def make_mlinterp(it, funname):
 
     if it.values =='vector':
         return None
+
     if it.eval in ('float', 'tuple') and it.values =='vector':
         # raise Exception("Non supported. (return type unknown)")
         return None
@@ -110,7 +130,7 @@ def {funname}(*args):
     grid = {grid_s}
     C = args[{it.d}]
     point = {point_s}
-    res = _mlinterp(grid, C, point)
+    res = mlinterp(grid, C, point)
     return res
     """
         return source
@@ -126,7 +146,7 @@ def {funname}(*args):
     res = zeros(N)
     # return res
     for n in range(N):
-        res[n] = _mlinterp(grid, C, {p_s})
+        res[n] = mlinterp(grid, C, {p_s})
     return res
 """
         return source
@@ -146,7 +166,7 @@ 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]))
+            res[n,m] = mlinterp(grid, C, (points_x[n], points_y[m]))
     return res
 """
         return source

interpolation/multilinear/tests/test_multilinear.py

@@ -50,20 +50,45 @@
 tests = [a2, a3, a4, c1, c2, c3, c4]
 tests_failing = [b1, b2, b3, b4, b5, d4, d5, d1, d2, d3, d4, d5]
 
-from interpolation.multilinear.mlinterp import mlinterp
+from interpolation.multilinear.mlinterp import mlinterp, interp
+
+
+def test_mlinterp():
+
+    # simple multilinear interpolation api
+
+    import numpy as np
+    from interpolation import mlinterp
+
+    # from interpolation.multilinear.mlinterp import mlininterp, mlininterp_vec
+    x1 = np.linspace(0,1,10)
+    x2 = np.linspace(0,1,20)
+    y = np.random.random((10,20))
+
+    z1 = np.linspace(0,1,30)
+    z2 = np.linspace(0,1,30)
+
+    pp = np.random.random((2000,2))
+
+    res0 = mlinterp((x1,x2), y, pp)
+    res0 = mlinterp((x1,x2), y, (0.1, 0.2))
+
 
 def test_multilinear():
 
+    # flat flexible api
+
     for t in tests:
 
         tt = [typeof(e) for e in t]
-        rr = mlinterp(*t)
+        print(tt)
+        rr = interp(*t)
+
         try:
             print(f"{tt}: {rr.shape}")
         except:
             print(f"{tt}: OK")
 
-
 #
 # exit()
 #