Fixes to object API + updated README.

Author: albop

README.md

@@ -2,56 +2,118 @@
 
 [![Build Status](https://travis-ci.org/EconForge/interpolation.py.svg?branch=master)](https://travis-ci.org/EconForge/interpolation.py)
 
-This library provides fast numba-accelerated interpolation routines
-for the following cases:
+The library contains:
+- splines.*: fast numba-compatible multilinear and cubic interpolation
+- multilinear.*: fast numba-compatible multilinear interpolation (alternative implementation)
+- smolyak.*: smolyak polinomials
+- complete polynomials
 
-## multilinear and cubic splines in any dimension
+## multilinear and cubic interpolation
+
+Fast numba-accelerated interpolation routines
+for multilinear and cubic interpolation, with any number of dimensions.
+Several interfaces are provided.
+
+### eval_linear
+
+Preferred interface for multilinear interpolation. It can interpolate on uniform
+and nonuniform cartesian grids. Several extrapolation options are available.
 
-Here is an example in 3 dimensions:
 
 ```python
-from interpolation.splines import LinearSpline, CubicSpline
-a = np.array([0.0,0.0,0.0])         # lower boundaries
-b = np.array([1.0,1.0,1.0])         # upper boundaries
-orders = np.array([50,50,50])       # 10 points along each dimension
-values = np.random.random(orders)   # values at each node of the grid
-S = np.random.random((10^6,3))    # coordinates at which to evaluate the splines
+import numpy as np
 
-# multilinear
-lin = LinearSpline(a,b,orders,values)
-V = lin(S)
-# cubic
-spline = CubicSpline(a,b,orders,values) # filter the coefficients
-V = spline(S)                       # interpolates -> (100000,) array
+from interpolation.splines import UCGrid, CGrid, nodes
 
-```
+# we interpolate function
+f = lambda x,y: np.sin(np.sqrt(x**2+y**2+0.00001))/np.sqrt(x**2+y**2+0.00001)
 
-Unfair timings: (from `misc/speed_comparison.py`)
-```
-# interpolate 10^6 points on a 50x50x50 grid.
-Cubic: 0.11488723754882812
-Linear: 0.03426337242126465
-Scipy (linear): 0.6502540111541748
-```
+# uniform cartesian grid
+grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
 
-More details are available as an example [notebook](https://github.com/EconForge/interpolation.py/blob/master/examples/cubic_splines_python.ipynb)
+# get grid points
+gp = nodes(grid)   # 100x2 matrix
+
+# compute values on grid points
+values = f(gp[:,0], gp[:,1]).reshape((10,10))
+
+from interpolation.splines import eval_linear
+# interpolate at one point
+point = np.array([0.1,0.45]) # 1d array
+val = eval_linear(grid, values, point)  # float
+
+# interpolate at many points:
+points = np.random.random((10000,2))
+eval_linear(grid, values, points) # 10000 vector
+
+# output can be preallocated
+out = np.zeros(10000)
+eval_linear(grid, values, points, out) # 10000 vector
+
+
+# default calls extrapolate data by using the nearest value inside the grid
+# other extrapolation options can be chosen among NEAREST, LINEAR, CONSTANT
+
+from interpolation.splines import extrap_options as xto
+points = np.random.random((100,2))*3-1
+eval_linear(grid, values, points, xto.NEAREST)  # 100
+eval_linear(grid, values, points, xto.LINEAR)   # 10000 vector
+eval_linear(grid, values, points, xto.CONSTANT) # 10000 vector
+
+
+# one can also approximate on nonuniform cartesian grids
+
+grid = CGrid(np.linspace(-1,1,100)**3, (-1.0, 1.0, 10))
+
+points = np.random.random((10000,2))
+eval_linear(grid, values, points) # 10000 vector
+
+
+# it is also possible to interpolate vector-valued functions in the following way
+
+f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
+g = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
+grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
+gp = nodes(grid)   # 100x2 matrix
+mvalues = np.concatenate([
+   f(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None],
+   g(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None]
+],axis=2) # 10x10x2 array
+points = np.random.random((1000,2))
+eval_linear(grid, mvalues, points[:,1]) # 2 elements vector
+eval_linear(grid, mvalues, points)      # 1000x2 matrix      
+out = np.zeros((1000,2))
+eval_linear(grid, mvalues, points, out) # 1000x2 matrix
 
-Other available features are:
-- linear extrapolation
-- computation of derivatives
-- interpolate many functions at once (multi-splines or vector valued splines)
 
-Experimental
-- evaluation on the GPU (with numba.cuda)
-- parallel evaluation (with guvectorize)
 
-In the near future:
-- JIT classes for all interpolation objects
+# finally, the same syntax can be used to interpolate using cubic splines
+# one just needs to prefilter the coefficients first
+# the same set of options apply but nonuniform grids are not supported (yet)
 
+f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
+grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
+gp = nodes(grid)   # 100x2 matrix
+values = f(gp[:,0], gp[:,1]).reshape((10,10))
 
-## jitted, non-uniform multilinear interpolation
+# filter values
+from interpolation.splines import filter_cubic
+coeffs = filter_cubic(grid, values) # a 12x12 array
+
+from interpolation.splines import eval_cubic
+points = np.random.random((1000,2))
+eval_cubic(grid, coeffs, points[:,1]) # 2 elements vector
+eval_cubic(grid, coeffs, points)      # 1000x2 matrix      
+out = np.zeros((1000,2))
+eval_cubic(grid, coeffs, points, out) # 1000x2 matrix
+
+```
+
+
+### interp
+
+Simpler interface. Mimmicks default `scipy.interp`: mutlilinear interpolation with constant extrapolation.
 
-There is a simple `interp` function with a flexible API which does multinear on uniform or non uniform cartesian grids.
 
 ```
 ### 1d grid
@@ -66,10 +128,51 @@ interp(x,y,0.5)
 # or at many points:
 u = np.linspace(0,1,1000)   # points
 interp(x,y,u)
+
+```
+
+
+### object interface
+
+This is for compatibility purpose only, until a new jittable model object is found.
+
+```python
+from interpolation.splines import LinearSpline, CubicSpline
+a = np.array([0.0,0.0,0.0])         # lower boundaries
+b = np.array([1.0,1.0,1.0])         # upper boundaries
+orders = np.array([50,50,50])       # 10 points along each dimension
+values = np.random.random(orders)   # values at each node of the grid
+S = np.random.random((10^6,3))    # coordinates at which to evaluate the splines
+
+# multilinear
+lin = LinearSpline(a,b,orders,values)
+V = lin(S)
+# cubic
+spline = CubicSpline(a,b,orders,values) # filter the coefficients
+V = spline(S)                       # interpolates -> (100000,) array
+
 ```
 
+### development notes
 
+Old, unfair timings: (from `misc/speed_comparison.py`)
 
+```
+# interpolate 10^6 points on a 50x50x50 grid.
+Cubic: 0.11488723754882812
+Linear: 0.03426337242126465
+Scipy (linear): 0.6502540111541748
+```
+
+More details are available as an example [notebook](https://github.com/EconForge/interpolation.py/blob/master/examples/cubic_splines_python.ipynb) (outdated)
+
+Missing but available soon:
+- splines at any order
+- derivative
+
+Feasible (some experiments)
+- evaluation on the GPU (with numba.cuda)
+- parallel evaluation (with guvectorize)
 
 
 

interpolation/splines/__init__.py

@@ -2,3 +2,13 @@
 from .multilinear import LinearSpline, LinearSplines
 from .eval_splines import options, eval_linear, eval_cubic
 from .prefilter_cubic import filter_cubic
+from .option_types import options as extrap_options
+
+# dummy functions
+def UCGrid(*args):
+    return tuple(args)
+def CGrid(*args):
+    return tuple(args)
+def nodes(grid):
+    from interpolation.cartesian import mlinspace
+    return mlinspace([g[0] for g in grid],[g[1] for g in grid],[g[2] for g in grid])

interpolation/splines/basis_splines.py

@@ -14,30 +14,30 @@ def B(p, u, i, x):
     else:
         return (((x - u[i]) / (u[i + p] - u[i])) * B(p - 1, u, i, x) + ((u[i + p + 1] - x) / (u[i + p + 1] - u[i + 1])) * B(p - 1, u, i + 1, x))
 
-
-from pylab import *
-
-xvec = linspace(-5, 5, 11)
-
-tvec = linspace(-5, 5, 1000)
-
-yvec = array([B0(xvec, 5, x) for x in tvec])
-y2vec = array([B(2, xvec, 5, x) for x in tvec])
-y3vec = array([B(3, xvec, 5, x) for x in tvec])
-
-
-print(B(3, xvec, 5, 0))
-print(B(3, xvec, 5, 1))
-print(B(3, xvec, 5, 2))
-print(B(3, xvec, 5, 3))
-print(B(3, xvec, 5, 4))
-
-
-plot(xvec, xvec * 0, 'o-')
-
-plot(tvec, yvec)
-plot(tvec, y2vec)
-plot(tvec, y3vec)
-
-
-show()
+#
+# from pylab import *
+#
+# xvec = linspace(-5, 5, 11)
+#
+# tvec = linspace(-5, 5, 1000)
+#
+# yvec = array([B0(xvec, 5, x) for x in tvec])
+# y2vec = array([B(2, xvec, 5, x) for x in tvec])
+# y3vec = array([B(3, xvec, 5, x) for x in tvec])
+#
+#
+# print(B(3, xvec, 5, 0))
+# print(B(3, xvec, 5, 1))
+# print(B(3, xvec, 5, 2))
+# print(B(3, xvec, 5, 3))
+# print(B(3, xvec, 5, 4))
+#
+#
+# plot(xvec, xvec * 0, 'o-')
+#
+# plot(tvec, yvec)
+# plot(tvec, y2vec)
+# plot(tvec, y3vec)
+#
+#
+# show()

interpolation/splines/eval_cubic.py

@@ -1,10 +1,11 @@
 import numpy
 
-
 from .eval_splines import eval_cubic
-## the functions in this file work for any dimension (d<=4)
-## they can optionnally allocate memory for the result
 
+## the functions in this file provide backward compatibility calls
+##
+## they can optionnally allocate memory for the result
+## they work for any dimension, except the functions which compute the gradient
 
 #######################
 # Compatibility calls #
@@ -21,7 +22,7 @@ def get_grid(a,b,n,C):
     print(txt)
     f = source_to_function(txt)
     return f
-    
+
 
 def eval_cubic_spline(a, b, orders, coefs, point):
     """Evaluates a cubic spline at one point

interpolation/splines/multilinear.py

@@ -79,9 +79,11 @@ def set_values(self,values):
 
     def interpolate(self,s):
 
-        from .multilinear_numba import multilinear_interpolation
+        from .eval_splines import eval_linear
+
         s = numpy.ascontiguousarray(s, dtype=self.dtype)
-        a = multilinear_interpolation(self.smin,self.smax,self.orders,self.values,s)
+        grid = tuple((self.smin[i], self.smax[i], self.orders[i]) for i in range(len(self.smin)))
+        a = eval_linear(grid,self.values,s)
         return a
 
     def __call__(self,s):
@@ -162,9 +164,10 @@ def set_values(self,mvalues):
 
     def interpolate(self,s):
 
-        from .multilinear_numba import vec_multilinear_interpolation
+        from .multilinear_numba import eval_linear
+        grid = tuple((self.smin[i], self.smax[i], self.orders[i]) for i in range(len(self.smin)))
         s = numpy.ascontiguousarray(s, dtype=self.dtype)
-        a = vec_multilinear_interpolation(self.smin,self.smax,self.orders,self.mvalues,s)
+        a = eval_linear(grid,self.mvalues,s)
         return a
 
     def __call__(self,s):

interpolation/splines/multilinear_numba.py

@@ -1,3 +1,5 @@
+# these are compatibility calls
+
 import numpy as np
 from numba import njit
 from .eval_splines import eval_linear
@@ -20,4 +22,3 @@ def multilinear_interpolation(smin, smax, orders, values, s, out=None):
         return eval_linear(grid, values, s)
     else:
         eval_linear(grid, values, s, out)
-        

interpolation/splines/splines.py

@@ -86,23 +86,26 @@ def interpolate(self, points, values=None, with_derivatives=False):
 
         import time
 
-        from .eval_cubic import vec_eval_cubic_spline, eval_cubic_spline
+        from .eval_splines import eval_cubic
+
+        grid = tuple((self.a[i], self.b[i], self.orders[i]) for i in range(len(self.a)))
 
         if not np.all( np.isfinite(points)):
             raise Exception('Spline interpolator evaluated at non-finite points.')
 
         if not with_derivatives:
             if points.ndim == 1:
                 # evaluate only on one point
-                return eval_cubic_spline(self.a, self.b, self.orders, self.__coeffs__, points)
+                return eval_cubic(grid, self.__coeffs__, points)
             else:
 
                 N, d = points.shape
                 assert(d==self.d)
                 if values is None:
                     values = np.empty(N, dtype=self.dtype)
-                vec_eval_cubic_spline(self.a, self.b, self.orders, self.__coeffs__, points, values)
+                eval_cubic(grid, self.__coeffs__, points, values)
                 return values
+
         else:
             raise Exception("Not implemented.")
 
@@ -182,10 +185,14 @@ def interpolate(self, points, diff=False):
         N = points.shape[0]
         d = points.shape[1]
 
+        from .eval_splines import eval_cubic
+
+
         if not diff:
-            from .eval_cubic import vec_eval_cubic_splines
+            grid = tuple((self.a[i], self.b[i], self.orders[i]) for i in range(len(self.a)))
+            from .eval_splines import eval_cubic
             values = np.empty((N,n_sp), dtype=float)
-            vec_eval_cubic_splines(self.a, self.b, self.orders, self.__mcoeffs__, points, values)
+            eval_cubic(grid, self.__mcoeffs__, points, values)
             return values
         else:
             from .eval_cubic import vec_eval_cubic_splines_G

interpolation/splines/tests/test_object_api.py

@@ -58,7 +58,7 @@ def ttest_object_api(Obj):
 def ttest_object_vector_api(Obj):
 
     cs = Obj(a,b,orders,mvals)
-    
+
     ii = cs(point)
     iii = cs(points)
     try: