|
1 |
| -# Optimized interpolation routines in Python / numba |
2 |
| - |
3 | 1 |  
|
4 | 2 |
|
5 | 3 |
|
6 |
| -The library contains: |
7 |
| -- `splines.*`: fast numba-compatible multilinear and cubic interpolation |
8 |
| -- `multilinear.*`: fast numba-compatible multilinear interpolation (alternative implementation) |
9 |
| -- `smolyak.*`: smolyak polynomials |
10 |
| -- complete polynomials |
11 |
| - |
12 |
| -## install |
13 |
| - |
14 |
| -Install latest version: |
15 |
| - |
16 |
| -- from conda: `conda install -c conda-forge interpolation` |
17 |
| -- from PyPI: `pip install interpolation` |
18 |
| - |
19 |
| -Latest development version from git: |
20 |
| - |
21 |
| -``` |
22 |
| -pip install poetry # if not already installed |
23 |
| -pip install git+https://github.com/econforge/interpolation.py.git/ |
24 |
| -``` |
25 |
| - |
26 |
| -The project uses a `pyproject.toml` file instead of `setup.py` and other legacy configuration files. For those used to development installation, this is feasible using `dephell`: |
27 |
| - |
28 |
| -``` |
29 |
| -pip install dephell # if not already installed |
30 |
| -dephell --from=pyproject.toml --to=setup.py # only once |
31 |
| -pip install -e . # like old times |
32 |
| -``` |
33 |
| - |
34 |
| -## multilinear and cubic interpolation |
35 |
| - |
36 |
| -Fast numba-accelerated interpolation routines |
37 |
| -for multilinear and cubic interpolation, with any number of dimensions. |
38 |
| -Several interfaces are provided. |
39 |
| - |
40 |
| -### eval_linear |
41 |
| - |
42 |
| -Preferred interface for multilinear interpolation. It can interpolate on uniform |
43 |
| -and nonuniform cartesian grids. Several extrapolation options are available. |
44 |
| - |
45 |
| - |
46 |
| -```python |
47 |
| -import numpy as np |
48 |
| - |
49 |
| -from interpolation.splines import UCGrid, CGrid, nodes |
50 |
| - |
51 |
| -# we interpolate function |
52 |
| -f = lambda x,y: np.sin(np.sqrt(x**2+y**2+0.00001))/np.sqrt(x**2+y**2+0.00001) |
53 |
| - |
54 |
| -# uniform cartesian grid |
55 |
| -grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10)) |
56 |
| - |
57 |
| -# get grid points |
58 |
| -gp = nodes(grid) # 100x2 matrix |
59 |
| - |
60 |
| -# compute values on grid points |
61 |
| -values = f(gp[:,0], gp[:,1]).reshape((10,10)) |
62 |
| - |
63 |
| -from interpolation.splines import eval_linear |
64 |
| -# interpolate at one point |
65 |
| -point = np.array([0.1,0.45]) # 1d array |
66 |
| -val = eval_linear(grid, values, point) # float |
67 |
| - |
68 |
| -# interpolate at many points: |
69 |
| -points = np.random.random((10000,2)) |
70 |
| -eval_linear(grid, values, points) # 10000 vector |
71 |
| - |
72 |
| -# output can be preallocated |
73 |
| -out = np.zeros(10000) |
74 |
| -eval_linear(grid, values, points, out) # 10000 vector |
75 |
| - |
76 |
| -## jitted, non-uniform multilinear interpolation |
77 |
| - |
78 |
| -# default calls extrapolate data by using the nearest value inside the grid |
79 |
| -# other extrapolation options can be chosen among NEAREST, LINEAR, CONSTANT |
80 |
| - |
81 |
| -from interpolation.splines import extrap_options as xto |
82 |
| -points = np.random.random((100,2))*3-1 |
83 |
| -eval_linear(grid, values, points, xto.NEAREST) # 100 |
84 |
| -eval_linear(grid, values, points, xto.LINEAR) # 10000 vector |
85 |
| -eval_linear(grid, values, points, xto.CONSTANT) # 10000 vector |
86 |
| - |
87 |
| - |
88 |
| -# one can also approximate on nonuniform cartesian grids |
89 |
| - |
90 |
| -grid = CGrid(np.linspace(-1,1,100)**3, (-1.0, 1.0, 10)) |
91 |
| - |
92 |
| -points = np.random.random((10000,2)) |
93 |
| -eval_linear(grid, values, points) # 10000 vector |
94 |
| - |
95 |
| - |
96 |
| -# it is also possible to interpolate vector-valued functions in the following way |
97 |
| - |
98 |
| -f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001) |
99 |
| -g = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001) |
100 |
| -grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10)) |
101 |
| -gp = nodes(grid) # 100x2 matrix |
102 |
| -mvalues = np.concatenate([ |
103 |
| - f(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None], |
104 |
| - g(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None] |
105 |
| -],axis=2) # 10x10x2 array |
106 |
| -points = np.random.random((1000,2)) |
107 |
| -eval_linear(grid, mvalues, points[:,1]) # 2 elements vector |
108 |
| -eval_linear(grid, mvalues, points) # 1000x2 matrix |
109 |
| -out = np.zeros((1000,2)) |
110 |
| -eval_linear(grid, mvalues, points, out) # 1000x2 matrix |
111 |
| - |
112 |
| - |
113 |
| - |
114 |
| -# finally, the same syntax can be used to interpolate using cubic splines |
115 |
| -# one just needs to prefilter the coefficients first |
116 |
| -# the same set of options apply but nonuniform grids are not supported (yet) |
117 |
| - |
118 |
| -f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001) |
119 |
| -grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10)) |
120 |
| -gp = nodes(grid) # 100x2 matrix |
121 |
| -values = f(gp[:,0], gp[:,1]).reshape((10,10)) |
122 |
| - |
123 |
| -# filter values |
124 |
| -from interpolation.splines import filter_cubic |
125 |
| -coeffs = filter_cubic(grid, values) # a 12x12 array |
126 |
| - |
127 |
| -from interpolation.splines import eval_cubic |
128 |
| -points = np.random.random((1000,2)) |
129 |
| -eval_cubic(grid, coeffs, points[:,1]) # 2 elements vector |
130 |
| -eval_cubic(grid, coeffs, points) # 1000x2 matrix |
131 |
| -out = np.zeros((1000,2)) |
132 |
| -eval_cubic(grid, coeffs, points, out) # 1000x2 matrix |
133 |
| - |
134 |
| -``` |
135 |
| - |
136 |
| -*Remark*: the arguably strange syntax for the extapolation option comes from the fact the actualy method called must be determined by type inference. So `eval_linear(..., extrap_method='linear')` would not work because the type of the last argument is always a string. Instead, we use opts.CONSTANT and opts.LINEAR for instance which have different numba types. |
137 |
| - |
138 |
| -Despite what it looks `UCGrid` and `CGRid` are not objects but functions which return very simple python structures that is a tuple of its arguments. For instance, `((0.0,1.0, 10), (0.0,1.0,20))` represents a 2d square discretized with 10 points along the first dimension and 20 along the second dimension. Similarly `(np.array([0.0, 0.1, 0.3, 1.0]), (0.0, 1.0, 20))` represents a square nonuniformly discretized along the first dimension (with 3 points) but uniformly along the second one. Now type dispatch is very sensitive to the exact types (floats vs ints), (tuple vs lists) which is potentially error-prone. Eventually, the functions `UCGrid` and `CGrid` will provide some type check and sensible conversions where it applies. This may change when if a parameterized structure-like object appear in numba. |
139 |
| - |
140 |
| -### interp |
141 |
| - |
142 |
| -Simpler interface. Mimmicks default `scipy.interp`: mutlilinear interpolation with constant extrapolation. |
143 |
| - |
144 |
| - |
145 |
| -``` |
146 |
| -### 1d grid |
147 |
| -from interpolation import interp |
148 |
| -
|
149 |
| -x = np.linspace(0,1,100)**2 # non-uniform points |
150 |
| -y = np.linspace(0,1,100) # values |
151 |
| -
|
152 |
| -# interpolate at one point: |
153 |
| -interp(x,y,0.5) |
154 |
| -
|
155 |
| -# or at many points: |
156 |
| -u = np.linspace(0,1,1000) # points |
157 |
| -interp(x,y,u) |
158 |
| -
|
159 |
| -``` |
160 |
| - |
161 |
| - |
162 |
| -### object interface |
163 |
| - |
164 |
| -This is for compatibility purpose only, until a new jittable model object is found. |
165 |
| - |
166 |
| -```python |
167 |
| -from interpolation.splines import LinearSpline, CubicSpline |
168 |
| -a = np.array([0.0,0.0,0.0]) # lower boundaries |
169 |
| -b = np.array([1.0,1.0,1.0]) # upper boundaries |
170 |
| -orders = np.array([50,50,50]) # 50 points along each dimension |
171 |
| -values = np.random.random(orders) # values at each node of the grid |
172 |
| -S = np.random.random((10**6,3)) # coordinates at which to evaluate the splines |
173 |
| - |
174 |
| -# multilinear |
175 |
| -lin = LinearSpline(a,b,orders,values) |
176 |
| -V = lin(S) |
177 |
| -# cubic |
178 |
| -spline = CubicSpline(a,b,orders,values) # filter the coefficients |
179 |
| -V = spline(S) # interpolates -> (100000,) array |
180 |
| - |
181 |
| -``` |
182 |
| - |
183 |
| -### development notes |
184 |
| - |
185 |
| -Old, unfair timings: (from `misc/speed_comparison.py`) |
186 |
| - |
187 |
| -``` |
188 |
| -# interpolate 10^6 points on a 50x50x50 grid. |
189 |
| -Cubic: 0.11488723754882812 |
190 |
| -Linear: 0.03426337242126465 |
191 |
| -Scipy (linear): 0.6502540111541748 |
192 |
| -``` |
193 |
| - |
194 |
| -More details are available as an example [notebook](https://github.com/EconForge/interpolation.py/blob/master/examples/cubic_splines_python.ipynb) (outdated) |
195 |
| - |
196 |
| -Missing but available soon: |
197 |
| -- splines at any order |
198 |
| - |
199 |
| -Feasible (some experiments) |
200 |
| -- evaluation on the GPU (with numba.cuda) |
201 |
| -- parallel evaluation (with guvectorize) |
202 |
| - |
203 |
| - |
204 |
| - |
205 |
| -## smolyak interpolation |
206 |
| - |
207 |
| -See [testfile](https://github.com/EconForge/interpolation.py/blob/master/interpolation/smolyak/tests/test_interp.py) for examples. |
| 4 | +See the [doc](https://www.econforge.org/interpolation.py). |
0 commit comments