1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
from __future__ import division
import numpy as np
from geotiepoints.multilinear_cython import multilinear_interpolation
def mlinspace(smin, smax, orders):
if len(orders) == 1:
res = np.atleast_2d(
np.linspace(np.array(smin), np.array(smax), np.array(orders)))
return res.copy() # workaround for strange bug
else:
meshes = np.meshgrid(
*[np.linspace(smin[i], smax[i], orders[i]) for i in range(len(orders))], indexing='ij')
return np.vstack([l.flatten() for l in meshes])
class MultilinearInterpolator:
"""Multilinear interpolation.
Methods
-------
smin, smax, orders : iterable objects
Specifies the boundaries of a cartesian grid, with the number of points along each dimension.
values : array_like (2D), optional
Each line enumerates values taken by a function on the cartesian grid, with last index varying faster.
Attributes
----------
smin, smax, orders :
Boundaries and number of points along each dimension.
d : number of dimensions
grid : array_like (2D)
Enumerate the approximation grid. Each column contains coordinates of a point in R^d.
values :
Values on the grid to be interpolated.
Example
-------
smin = [-1,-1]
smax = [1,1]
orders = [5,5]
f = lambda x: np.vstack([np.sqrt( x[0,:]**2 + x[1,:]**2 ),
np.power( x[0,:]**3 + x[1,:]**3, 1.0/3.0 )])
interp = MultilinearInterpolator(smin,smax,orders)
interp.set_values( f(interp.grid) )
random_points = np.random.random( (2, 1000) )
interpolated_values = interp(random_points)
exact_values = f(random_points)
"""
__grid__ = None
def __init__(self, smin, smax, orders, values=None, dtype=np.float64):
self.smin = np.asarray(smin, dtype=dtype)
self.smax = np.asarray(smax, dtype=dtype)
self.orders = np.asarray(orders, dtype="long")
self.d = len(orders)
self.dtype = dtype
if values is not None:
self.set_values(values)
@property
def grid(self):
if self.__grid__ is None:
self.__grid__ = mlinspace(self.smin, self.smax, self.orders)
return self.__grid__
def set_values(self, values):
self.values = np.ascontiguousarray(values, dtype=self.dtype)
def interpolate(self, s):
s = np.ascontiguousarray(s, dtype=self.dtype)
a = multilinear_interpolation(
self.smin, self.smax, self.orders, self.values, s)
return a
def __call__(self, s):
return self.interpolate(s)
|
