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