Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/interpolate/rbf.py : 19%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""rbf - Radial basis functions for interpolation/smoothing scattered N-D data.
3Written by John Travers <jtravs@gmail.com>, February 2007
4Based closely on Matlab code by Alex Chirokov
5Additional, large, improvements by Robert Hetland
6Some additional alterations by Travis Oliphant
7Interpolation with multi-dimensional target domain by Josua Sassen
9Permission to use, modify, and distribute this software is given under the
10terms of the SciPy (BSD style) license. See LICENSE.txt that came with
11this distribution for specifics.
13NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
15Copyright (c) 2006-2007, Robert Hetland <hetland@tamu.edu>
16Copyright (c) 2007, John Travers <jtravs@gmail.com>
18Redistribution and use in source and binary forms, with or without
19modification, are permitted provided that the following conditions are
20met:
22 * Redistributions of source code must retain the above copyright
23 notice, this list of conditions and the following disclaimer.
25 * Redistributions in binary form must reproduce the above
26 copyright notice, this list of conditions and the following
27 disclaimer in the documentation and/or other materials provided
28 with the distribution.
30 * Neither the name of Robert Hetland nor the names of any
31 contributors may be used to endorse or promote products derived
32 from this software without specific prior written permission.
34THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
35"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
36LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
37A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
38OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
39SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
40LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
41DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
42THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
43(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
44OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
45"""
46import numpy as np
48from scipy import linalg
49from scipy.special import xlogy
50from scipy.spatial.distance import cdist, pdist, squareform
52__all__ = ['Rbf']
55class Rbf(object):
56 """
57 Rbf(*args)
59 A class for radial basis function interpolation of functions from
60 N-D scattered data to an M-D domain.
62 Parameters
63 ----------
64 *args : arrays
65 x, y, z, ..., d, where x, y, z, ... are the coordinates of the nodes
66 and d is the array of values at the nodes
67 function : str or callable, optional
68 The radial basis function, based on the radius, r, given by the norm
69 (default is Euclidean distance); the default is 'multiquadric'::
71 'multiquadric': sqrt((r/self.epsilon)**2 + 1)
72 'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1)
73 'gaussian': exp(-(r/self.epsilon)**2)
74 'linear': r
75 'cubic': r**3
76 'quintic': r**5
77 'thin_plate': r**2 * log(r)
79 If callable, then it must take 2 arguments (self, r). The epsilon
80 parameter will be available as self.epsilon. Other keyword
81 arguments passed in will be available as well.
83 epsilon : float, optional
84 Adjustable constant for gaussian or multiquadrics functions
85 - defaults to approximate average distance between nodes (which is
86 a good start).
87 smooth : float, optional
88 Values greater than zero increase the smoothness of the
89 approximation. 0 is for interpolation (default), the function will
90 always go through the nodal points in this case.
91 norm : str, callable, optional
92 A function that returns the 'distance' between two points, with
93 inputs as arrays of positions (x, y, z, ...), and an output as an
94 array of distance. E.g., the default: 'euclidean', such that the result
95 is a matrix of the distances from each point in ``x1`` to each point in
96 ``x2``. For more options, see documentation of
97 `scipy.spatial.distances.cdist`.
98 mode : str, optional
99 Mode of the interpolation, can be '1-D' (default) or 'N-D'. When it is
100 '1-D' the data `d` will be considered as 1-D and flattened
101 internally. When it is 'N-D' the data `d` is assumed to be an array of
102 shape (n_samples, m), where m is the dimension of the target domain.
105 Attributes
106 ----------
107 N : int
108 The number of data points (as determined by the input arrays).
109 di : ndarray
110 The 1-D array of data values at each of the data coordinates `xi`.
111 xi : ndarray
112 The 2-D array of data coordinates.
113 function : str or callable
114 The radial basis function. See description under Parameters.
115 epsilon : float
116 Parameter used by gaussian or multiquadrics functions. See Parameters.
117 smooth : float
118 Smoothing parameter. See description under Parameters.
119 norm : str or callable
120 The distance function. See description under Parameters.
121 mode : str
122 Mode of the interpolation. See description under Parameters.
123 nodes : ndarray
124 A 1-D array of node values for the interpolation.
125 A : internal property, do not use
127 Examples
128 --------
129 >>> from scipy.interpolate import Rbf
130 >>> x, y, z, d = np.random.rand(4, 50)
131 >>> rbfi = Rbf(x, y, z, d) # radial basis function interpolator instance
132 >>> xi = yi = zi = np.linspace(0, 1, 20)
133 >>> di = rbfi(xi, yi, zi) # interpolated values
134 >>> di.shape
135 (20,)
137 """
138 # Available radial basis functions that can be selected as strings;
139 # they all start with _h_ (self._init_function relies on that)
140 def _h_multiquadric(self, r):
141 return np.sqrt((1.0/self.epsilon*r)**2 + 1)
143 def _h_inverse_multiquadric(self, r):
144 return 1.0/np.sqrt((1.0/self.epsilon*r)**2 + 1)
146 def _h_gaussian(self, r):
147 return np.exp(-(1.0/self.epsilon*r)**2)
149 def _h_linear(self, r):
150 return r
152 def _h_cubic(self, r):
153 return r**3
155 def _h_quintic(self, r):
156 return r**5
158 def _h_thin_plate(self, r):
159 return xlogy(r**2, r)
161 # Setup self._function and do smoke test on initial r
162 def _init_function(self, r):
163 if isinstance(self.function, str):
164 self.function = self.function.lower()
165 _mapped = {'inverse': 'inverse_multiquadric',
166 'inverse multiquadric': 'inverse_multiquadric',
167 'thin-plate': 'thin_plate'}
168 if self.function in _mapped:
169 self.function = _mapped[self.function]
171 func_name = "_h_" + self.function
172 if hasattr(self, func_name):
173 self._function = getattr(self, func_name)
174 else:
175 functionlist = [x[3:] for x in dir(self)
176 if x.startswith('_h_')]
177 raise ValueError("function must be a callable or one of " +
178 ", ".join(functionlist))
179 self._function = getattr(self, "_h_"+self.function)
180 elif callable(self.function):
181 allow_one = False
182 if hasattr(self.function, 'func_code') or \
183 hasattr(self.function, '__code__'):
184 val = self.function
185 allow_one = True
186 elif hasattr(self.function, "__call__"):
187 val = self.function.__call__.__func__
188 else:
189 raise ValueError("Cannot determine number of arguments to "
190 "function")
192 argcount = val.__code__.co_argcount
193 if allow_one and argcount == 1:
194 self._function = self.function
195 elif argcount == 2:
196 self._function = self.function.__get__(self, Rbf)
197 else:
198 raise ValueError("Function argument must take 1 or 2 "
199 "arguments.")
201 a0 = self._function(r)
202 if a0.shape != r.shape:
203 raise ValueError("Callable must take array and return array of "
204 "the same shape")
205 return a0
207 def __init__(self, *args, **kwargs):
208 # `args` can be a variable number of arrays; we flatten them and store
209 # them as a single 2-D array `xi` of shape (n_args-1, array_size),
210 # plus a 1-D array `di` for the values.
211 # All arrays must have the same number of elements
212 self.xi = np.asarray([np.asarray(a, dtype=np.float_).flatten()
213 for a in args[:-1]])
214 self.N = self.xi.shape[-1]
216 self.mode = kwargs.pop('mode', '1-D')
218 if self.mode == '1-D':
219 self.di = np.asarray(args[-1]).flatten()
220 self._target_dim = 1
221 elif self.mode == 'N-D':
222 self.di = np.asarray(args[-1])
223 self._target_dim = self.di.shape[-1]
224 else:
225 raise ValueError("Mode has to be 1-D or N-D.")
227 if not all([x.size == self.di.shape[0] for x in self.xi]):
228 raise ValueError("All arrays must be equal length.")
230 self.norm = kwargs.pop('norm', 'euclidean')
231 self.epsilon = kwargs.pop('epsilon', None)
232 if self.epsilon is None:
233 # default epsilon is the "the average distance between nodes" based
234 # on a bounding hypercube
235 ximax = np.amax(self.xi, axis=1)
236 ximin = np.amin(self.xi, axis=1)
237 edges = ximax - ximin
238 edges = edges[np.nonzero(edges)]
239 self.epsilon = np.power(np.prod(edges)/self.N, 1.0/edges.size)
241 self.smooth = kwargs.pop('smooth', 0.0)
242 self.function = kwargs.pop('function', 'multiquadric')
244 # attach anything left in kwargs to self for use by any user-callable
245 # function or to save on the object returned.
246 for item, value in kwargs.items():
247 setattr(self, item, value)
249 # Compute weights
250 if self._target_dim > 1: # If we have more than one target dimension,
251 # we first factorize the matrix
252 self.nodes = np.zeros((self.N, self._target_dim), dtype=self.di.dtype)
253 lu, piv = linalg.lu_factor(self.A)
254 for i in range(self._target_dim):
255 self.nodes[:, i] = linalg.lu_solve((lu, piv), self.di[:, i])
256 else:
257 self.nodes = linalg.solve(self.A, self.di)
259 @property
260 def A(self):
261 # this only exists for backwards compatibility: self.A was available
262 # and, at least technically, public.
263 r = squareform(pdist(self.xi.T, self.norm)) # Pairwise norm
264 return self._init_function(r) - np.eye(self.N)*self.smooth
266 def _call_norm(self, x1, x2):
267 return cdist(x1.T, x2.T, self.norm)
269 def __call__(self, *args):
270 args = [np.asarray(x) for x in args]
271 if not all([x.shape == y.shape for x in args for y in args]):
272 raise ValueError("Array lengths must be equal")
273 if self._target_dim > 1:
274 shp = args[0].shape + (self._target_dim,)
275 else:
276 shp = args[0].shape
277 xa = np.asarray([a.flatten() for a in args], dtype=np.float_)
278 r = self._call_norm(xa, self.xi)
279 return np.dot(self._function(r), self.nodes).reshape(shp)