Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_rbf.py: 20%
97 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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:
56 """
57 Rbf(*args, **kwargs)
59 A class for radial basis function interpolation of functions from
60 N-D scattered data to an M-D domain.
62 .. note::
63 `Rbf` is legacy code, for new usage please use `RBFInterpolator`
64 instead.
66 Parameters
67 ----------
68 *args : arrays
69 x, y, z, ..., d, where x, y, z, ... are the coordinates of the nodes
70 and d is the array of values at the nodes
71 function : str or callable, optional
72 The radial basis function, based on the radius, r, given by the norm
73 (default is Euclidean distance); the default is 'multiquadric'::
75 'multiquadric': sqrt((r/self.epsilon)**2 + 1)
76 'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1)
77 'gaussian': exp(-(r/self.epsilon)**2)
78 'linear': r
79 'cubic': r**3
80 'quintic': r**5
81 'thin_plate': r**2 * log(r)
83 If callable, then it must take 2 arguments (self, r). The epsilon
84 parameter will be available as self.epsilon. Other keyword
85 arguments passed in will be available as well.
87 epsilon : float, optional
88 Adjustable constant for gaussian or multiquadrics functions
89 - defaults to approximate average distance between nodes (which is
90 a good start).
91 smooth : float, optional
92 Values greater than zero increase the smoothness of the
93 approximation. 0 is for interpolation (default), the function will
94 always go through the nodal points in this case.
95 norm : str, callable, optional
96 A function that returns the 'distance' between two points, with
97 inputs as arrays of positions (x, y, z, ...), and an output as an
98 array of distance. E.g., the default: 'euclidean', such that the result
99 is a matrix of the distances from each point in ``x1`` to each point in
100 ``x2``. For more options, see documentation of
101 `scipy.spatial.distances.cdist`.
102 mode : str, optional
103 Mode of the interpolation, can be '1-D' (default) or 'N-D'. When it is
104 '1-D' the data `d` will be considered as 1-D and flattened
105 internally. When it is 'N-D' the data `d` is assumed to be an array of
106 shape (n_samples, m), where m is the dimension of the target domain.
109 Attributes
110 ----------
111 N : int
112 The number of data points (as determined by the input arrays).
113 di : ndarray
114 The 1-D array of data values at each of the data coordinates `xi`.
115 xi : ndarray
116 The 2-D array of data coordinates.
117 function : str or callable
118 The radial basis function. See description under Parameters.
119 epsilon : float
120 Parameter used by gaussian or multiquadrics functions. See Parameters.
121 smooth : float
122 Smoothing parameter. See description under Parameters.
123 norm : str or callable
124 The distance function. See description under Parameters.
125 mode : str
126 Mode of the interpolation. See description under Parameters.
127 nodes : ndarray
128 A 1-D array of node values for the interpolation.
129 A : internal property, do not use
131 See Also
132 --------
133 RBFInterpolator
135 Examples
136 --------
137 >>> import numpy as np
138 >>> from scipy.interpolate import Rbf
139 >>> rng = np.random.default_rng()
140 >>> x, y, z, d = rng.random((4, 50))
141 >>> rbfi = Rbf(x, y, z, d) # radial basis function interpolator instance
142 >>> xi = yi = zi = np.linspace(0, 1, 20)
143 >>> di = rbfi(xi, yi, zi) # interpolated values
144 >>> di.shape
145 (20,)
147 """
148 # Available radial basis functions that can be selected as strings;
149 # they all start with _h_ (self._init_function relies on that)
150 def _h_multiquadric(self, r):
151 return np.sqrt((1.0/self.epsilon*r)**2 + 1)
153 def _h_inverse_multiquadric(self, r):
154 return 1.0/np.sqrt((1.0/self.epsilon*r)**2 + 1)
156 def _h_gaussian(self, r):
157 return np.exp(-(1.0/self.epsilon*r)**2)
159 def _h_linear(self, r):
160 return r
162 def _h_cubic(self, r):
163 return r**3
165 def _h_quintic(self, r):
166 return r**5
168 def _h_thin_plate(self, r):
169 return xlogy(r**2, r)
171 # Setup self._function and do smoke test on initial r
172 def _init_function(self, r):
173 if isinstance(self.function, str):
174 self.function = self.function.lower()
175 _mapped = {'inverse': 'inverse_multiquadric',
176 'inverse multiquadric': 'inverse_multiquadric',
177 'thin-plate': 'thin_plate'}
178 if self.function in _mapped:
179 self.function = _mapped[self.function]
181 func_name = "_h_" + self.function
182 if hasattr(self, func_name):
183 self._function = getattr(self, func_name)
184 else:
185 functionlist = [x[3:] for x in dir(self)
186 if x.startswith('_h_')]
187 raise ValueError("function must be a callable or one of " +
188 ", ".join(functionlist))
189 self._function = getattr(self, "_h_"+self.function)
190 elif callable(self.function):
191 allow_one = False
192 if hasattr(self.function, 'func_code') or \
193 hasattr(self.function, '__code__'):
194 val = self.function
195 allow_one = True
196 elif hasattr(self.function, "__call__"):
197 val = self.function.__call__.__func__
198 else:
199 raise ValueError("Cannot determine number of arguments to "
200 "function")
202 argcount = val.__code__.co_argcount
203 if allow_one and argcount == 1:
204 self._function = self.function
205 elif argcount == 2:
206 self._function = self.function.__get__(self, Rbf)
207 else:
208 raise ValueError("Function argument must take 1 or 2 "
209 "arguments.")
211 a0 = self._function(r)
212 if a0.shape != r.shape:
213 raise ValueError("Callable must take array and return array of "
214 "the same shape")
215 return a0
217 def __init__(self, *args, **kwargs):
218 # `args` can be a variable number of arrays; we flatten them and store
219 # them as a single 2-D array `xi` of shape (n_args-1, array_size),
220 # plus a 1-D array `di` for the values.
221 # All arrays must have the same number of elements
222 self.xi = np.asarray([np.asarray(a, dtype=np.float_).flatten()
223 for a in args[:-1]])
224 self.N = self.xi.shape[-1]
226 self.mode = kwargs.pop('mode', '1-D')
228 if self.mode == '1-D':
229 self.di = np.asarray(args[-1]).flatten()
230 self._target_dim = 1
231 elif self.mode == 'N-D':
232 self.di = np.asarray(args[-1])
233 self._target_dim = self.di.shape[-1]
234 else:
235 raise ValueError("Mode has to be 1-D or N-D.")
237 if not all([x.size == self.di.shape[0] for x in self.xi]):
238 raise ValueError("All arrays must be equal length.")
240 self.norm = kwargs.pop('norm', 'euclidean')
241 self.epsilon = kwargs.pop('epsilon', None)
242 if self.epsilon is None:
243 # default epsilon is the "the average distance between nodes" based
244 # on a bounding hypercube
245 ximax = np.amax(self.xi, axis=1)
246 ximin = np.amin(self.xi, axis=1)
247 edges = ximax - ximin
248 edges = edges[np.nonzero(edges)]
249 self.epsilon = np.power(np.prod(edges)/self.N, 1.0/edges.size)
251 self.smooth = kwargs.pop('smooth', 0.0)
252 self.function = kwargs.pop('function', 'multiquadric')
254 # attach anything left in kwargs to self for use by any user-callable
255 # function or to save on the object returned.
256 for item, value in kwargs.items():
257 setattr(self, item, value)
259 # Compute weights
260 if self._target_dim > 1: # If we have more than one target dimension,
261 # we first factorize the matrix
262 self.nodes = np.zeros((self.N, self._target_dim), dtype=self.di.dtype)
263 lu, piv = linalg.lu_factor(self.A)
264 for i in range(self._target_dim):
265 self.nodes[:, i] = linalg.lu_solve((lu, piv), self.di[:, i])
266 else:
267 self.nodes = linalg.solve(self.A, self.di)
269 @property
270 def A(self):
271 # this only exists for backwards compatibility: self.A was available
272 # and, at least technically, public.
273 r = squareform(pdist(self.xi.T, self.norm)) # Pairwise norm
274 return self._init_function(r) - np.eye(self.N)*self.smooth
276 def _call_norm(self, x1, x2):
277 return cdist(x1.T, x2.T, self.norm)
279 def __call__(self, *args):
280 args = [np.asarray(x) for x in args]
281 if not all([x.shape == y.shape for x in args for y in args]):
282 raise ValueError("Array lengths must be equal")
283 if self._target_dim > 1:
284 shp = args[0].shape + (self._target_dim,)
285 else:
286 shp = args[0].shape
287 xa = np.asarray([a.flatten() for a in args], dtype=np.float_)
288 r = self._call_norm(xa, self.xi)
289 return np.dot(self._function(r), self.nodes).reshape(shp)