Merge pull request #38 from EconForge/albop/code_generation

Albop/code generation

Author: albop

README.md

@@ -2,56 +2,119 @@
 
 [![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])       # 50 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
 
-Other available features are:
-- linear extrapolation
-- computation of derivatives
-- interpolate many functions at once (multi-splines or vector valued splines)
+# compute values on grid points
+values = f(gp[:,0], gp[:,1]).reshape((10,10))
 
-Experimental
-- evaluation on the GPU (with numba.cuda)
-- parallel evaluation (with guvectorize)
+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
 
-In the near future:
-- JIT classes for all interpolation objects
+# 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
 
 ## jitted, non-uniform multilinear interpolation
 
-There is a simple `interp` function with a flexible API which does multinear on uniform or non uniform cartesian grids.
+# 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
+
+
+
+# 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))
+
+# 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.
+
 
 ```
 ### 1d grid
@@ -66,10 +129,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])       # 50 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/multilinear/fungen.py

@@ -185,7 +185,8 @@ def gen_tensor_reduction(X, symbs, inds=[]):
 
 @generated_jit(nopython=True)
 def tensor_reduction(C,l):
-    ex = gen_tensor_reduction('C', ['l[{}]'.format(i) for i in range(len(l.types))])
+    d = len(l.types)
+    ex = gen_tensor_reduction('C', ['l[{}]'.format(i) for i in range(d)])
     dd = dict()
     s = "def tensor_reduction(C,l): return {}".format(ex)
     eval(compile(ast.parse(s),'<string>','exec'), dd)

interpolation/splines/__init__.py

@@ -1,2 +1,14 @@
 from .splines import CubicSpline, CubicSplines
 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()