Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_rbfinterp.py: 13%
136 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Module for RBF interpolation."""
2import warnings
3from itertools import combinations_with_replacement
5import numpy as np
6from numpy.linalg import LinAlgError
7from scipy.spatial import KDTree
8from scipy.special import comb
9from scipy.linalg.lapack import dgesv # type: ignore[attr-defined]
11from ._rbfinterp_pythran import (_build_system,
12 _build_evaluation_coefficients,
13 _polynomial_matrix)
16__all__ = ["RBFInterpolator"]
19# These RBFs are implemented.
20_AVAILABLE = {
21 "linear",
22 "thin_plate_spline",
23 "cubic",
24 "quintic",
25 "multiquadric",
26 "inverse_multiquadric",
27 "inverse_quadratic",
28 "gaussian"
29 }
32# The shape parameter does not need to be specified when using these RBFs.
33_SCALE_INVARIANT = {"linear", "thin_plate_spline", "cubic", "quintic"}
36# For RBFs that are conditionally positive definite of order m, the interpolant
37# should include polynomial terms with degree >= m - 1. Define the minimum
38# degrees here. These values are from Chapter 8 of Fasshauer's "Meshfree
39# Approximation Methods with MATLAB". The RBFs that are not in this dictionary
40# are positive definite and do not need polynomial terms.
41_NAME_TO_MIN_DEGREE = {
42 "multiquadric": 0,
43 "linear": 0,
44 "thin_plate_spline": 1,
45 "cubic": 1,
46 "quintic": 2
47 }
50def _monomial_powers(ndim, degree):
51 """Return the powers for each monomial in a polynomial.
53 Parameters
54 ----------
55 ndim : int
56 Number of variables in the polynomial.
57 degree : int
58 Degree of the polynomial.
60 Returns
61 -------
62 (nmonos, ndim) int ndarray
63 Array where each row contains the powers for each variable in a
64 monomial.
66 """
67 nmonos = comb(degree + ndim, ndim, exact=True)
68 out = np.zeros((nmonos, ndim), dtype=int)
69 count = 0
70 for deg in range(degree + 1):
71 for mono in combinations_with_replacement(range(ndim), deg):
72 # `mono` is a tuple of variables in the current monomial with
73 # multiplicity indicating power (e.g., (0, 1, 1) represents x*y**2)
74 for var in mono:
75 out[count, var] += 1
77 count += 1
79 return out
82def _build_and_solve_system(y, d, smoothing, kernel, epsilon, powers):
83 """Build and solve the RBF interpolation system of equations.
85 Parameters
86 ----------
87 y : (P, N) float ndarray
88 Data point coordinates.
89 d : (P, S) float ndarray
90 Data values at `y`.
91 smoothing : (P,) float ndarray
92 Smoothing parameter for each data point.
93 kernel : str
94 Name of the RBF.
95 epsilon : float
96 Shape parameter.
97 powers : (R, N) int ndarray
98 The exponents for each monomial in the polynomial.
100 Returns
101 -------
102 coeffs : (P + R, S) float ndarray
103 Coefficients for each RBF and monomial.
104 shift : (N,) float ndarray
105 Domain shift used to create the polynomial matrix.
106 scale : (N,) float ndarray
107 Domain scaling used to create the polynomial matrix.
109 """
110 lhs, rhs, shift, scale = _build_system(
111 y, d, smoothing, kernel, epsilon, powers
112 )
113 _, _, coeffs, info = dgesv(lhs, rhs, overwrite_a=True, overwrite_b=True)
114 if info < 0:
115 raise ValueError(f"The {-info}-th argument had an illegal value.")
116 elif info > 0:
117 msg = "Singular matrix."
118 nmonos = powers.shape[0]
119 if nmonos > 0:
120 pmat = _polynomial_matrix((y - shift)/scale, powers)
121 rank = np.linalg.matrix_rank(pmat)
122 if rank < nmonos:
123 msg = (
124 "Singular matrix. The matrix of monomials evaluated at "
125 "the data point coordinates does not have full column "
126 f"rank ({rank}/{nmonos})."
127 )
129 raise LinAlgError(msg)
131 return shift, scale, coeffs
134class RBFInterpolator:
135 """Radial basis function (RBF) interpolation in N dimensions.
137 Parameters
138 ----------
139 y : (npoints, ndims) array_like
140 2-D array of data point coordinates.
141 d : (npoints, ...) array_like
142 N-D array of data values at `y`. The length of `d` along the first
143 axis must be equal to the length of `y`. Unlike some interpolators, the
144 interpolation axis cannot be changed.
145 neighbors : int, optional
146 If specified, the value of the interpolant at each evaluation point
147 will be computed using only this many nearest data points. All the data
148 points are used by default.
149 smoothing : float or (npoints, ) array_like, optional
150 Smoothing parameter. The interpolant perfectly fits the data when this
151 is set to 0. For large values, the interpolant approaches a least
152 squares fit of a polynomial with the specified degree. Default is 0.
153 kernel : str, optional
154 Type of RBF. This should be one of
156 - 'linear' : ``-r``
157 - 'thin_plate_spline' : ``r**2 * log(r)``
158 - 'cubic' : ``r**3``
159 - 'quintic' : ``-r**5``
160 - 'multiquadric' : ``-sqrt(1 + r**2)``
161 - 'inverse_multiquadric' : ``1/sqrt(1 + r**2)``
162 - 'inverse_quadratic' : ``1/(1 + r**2)``
163 - 'gaussian' : ``exp(-r**2)``
165 Default is 'thin_plate_spline'.
166 epsilon : float, optional
167 Shape parameter that scales the input to the RBF. If `kernel` is
168 'linear', 'thin_plate_spline', 'cubic', or 'quintic', this defaults to
169 1 and can be ignored because it has the same effect as scaling the
170 smoothing parameter. Otherwise, this must be specified.
171 degree : int, optional
172 Degree of the added polynomial. For some RBFs the interpolant may not
173 be well-posed if the polynomial degree is too small. Those RBFs and
174 their corresponding minimum degrees are
176 - 'multiquadric' : 0
177 - 'linear' : 0
178 - 'thin_plate_spline' : 1
179 - 'cubic' : 1
180 - 'quintic' : 2
182 The default value is the minimum degree for `kernel` or 0 if there is
183 no minimum degree. Set this to -1 for no added polynomial.
185 Notes
186 -----
187 An RBF is a scalar valued function in N-dimensional space whose value at
188 :math:`x` can be expressed in terms of :math:`r=||x - c||`, where :math:`c`
189 is the center of the RBF.
191 An RBF interpolant for the vector of data values :math:`d`, which are from
192 locations :math:`y`, is a linear combination of RBFs centered at :math:`y`
193 plus a polynomial with a specified degree. The RBF interpolant is written
194 as
196 .. math::
197 f(x) = K(x, y) a + P(x) b,
199 where :math:`K(x, y)` is a matrix of RBFs with centers at :math:`y`
200 evaluated at the points :math:`x`, and :math:`P(x)` is a matrix of
201 monomials, which span polynomials with the specified degree, evaluated at
202 :math:`x`. The coefficients :math:`a` and :math:`b` are the solution to the
203 linear equations
205 .. math::
206 (K(y, y) + \\lambda I) a + P(y) b = d
208 and
210 .. math::
211 P(y)^T a = 0,
213 where :math:`\\lambda` is a non-negative smoothing parameter that controls
214 how well we want to fit the data. The data are fit exactly when the
215 smoothing parameter is 0.
217 The above system is uniquely solvable if the following requirements are
218 met:
220 - :math:`P(y)` must have full column rank. :math:`P(y)` always has full
221 column rank when `degree` is -1 or 0. When `degree` is 1,
222 :math:`P(y)` has full column rank if the data point locations are not
223 all collinear (N=2), coplanar (N=3), etc.
224 - If `kernel` is 'multiquadric', 'linear', 'thin_plate_spline',
225 'cubic', or 'quintic', then `degree` must not be lower than the
226 minimum value listed above.
227 - If `smoothing` is 0, then each data point location must be distinct.
229 When using an RBF that is not scale invariant ('multiquadric',
230 'inverse_multiquadric', 'inverse_quadratic', or 'gaussian'), an appropriate
231 shape parameter must be chosen (e.g., through cross validation). Smaller
232 values for the shape parameter correspond to wider RBFs. The problem can
233 become ill-conditioned or singular when the shape parameter is too small.
235 The memory required to solve for the RBF interpolation coefficients
236 increases quadratically with the number of data points, which can become
237 impractical when interpolating more than about a thousand data points.
238 To overcome memory limitations for large interpolation problems, the
239 `neighbors` argument can be specified to compute an RBF interpolant for
240 each evaluation point using only the nearest data points.
242 .. versionadded:: 1.7.0
244 See Also
245 --------
246 NearestNDInterpolator
247 LinearNDInterpolator
248 CloughTocher2DInterpolator
250 References
251 ----------
252 .. [1] Fasshauer, G., 2007. Meshfree Approximation Methods with Matlab.
253 World Scientific Publishing Co.
255 .. [2] http://amadeus.math.iit.edu/~fass/603_ch3.pdf
257 .. [3] Wahba, G., 1990. Spline Models for Observational Data. SIAM.
259 .. [4] http://pages.stat.wisc.edu/~wahba/stat860public/lect/lect8/lect8.pdf
261 Examples
262 --------
263 Demonstrate interpolating scattered data to a grid in 2-D.
265 >>> import numpy as np
266 >>> import matplotlib.pyplot as plt
267 >>> from scipy.interpolate import RBFInterpolator
268 >>> from scipy.stats.qmc import Halton
270 >>> rng = np.random.default_rng()
271 >>> xobs = 2*Halton(2, seed=rng).random(100) - 1
272 >>> yobs = np.sum(xobs, axis=1)*np.exp(-6*np.sum(xobs**2, axis=1))
274 >>> xgrid = np.mgrid[-1:1:50j, -1:1:50j]
275 >>> xflat = xgrid.reshape(2, -1).T
276 >>> yflat = RBFInterpolator(xobs, yobs)(xflat)
277 >>> ygrid = yflat.reshape(50, 50)
279 >>> fig, ax = plt.subplots()
280 >>> ax.pcolormesh(*xgrid, ygrid, vmin=-0.25, vmax=0.25, shading='gouraud')
281 >>> p = ax.scatter(*xobs.T, c=yobs, s=50, ec='k', vmin=-0.25, vmax=0.25)
282 >>> fig.colorbar(p)
283 >>> plt.show()
285 """
287 def __init__(self, y, d,
288 neighbors=None,
289 smoothing=0.0,
290 kernel="thin_plate_spline",
291 epsilon=None,
292 degree=None):
293 y = np.asarray(y, dtype=float, order="C")
294 if y.ndim != 2:
295 raise ValueError("`y` must be a 2-dimensional array.")
297 ny, ndim = y.shape
299 d_dtype = complex if np.iscomplexobj(d) else float
300 d = np.asarray(d, dtype=d_dtype, order="C")
301 if d.shape[0] != ny:
302 raise ValueError(
303 f"Expected the first axis of `d` to have length {ny}."
304 )
306 d_shape = d.shape[1:]
307 d = d.reshape((ny, -1))
308 # If `d` is complex, convert it to a float array with twice as many
309 # columns. Otherwise, the LHS matrix would need to be converted to
310 # complex and take up 2x more memory than necessary.
311 d = d.view(float)
313 if np.isscalar(smoothing):
314 smoothing = np.full(ny, smoothing, dtype=float)
315 else:
316 smoothing = np.asarray(smoothing, dtype=float, order="C")
317 if smoothing.shape != (ny,):
318 raise ValueError(
319 "Expected `smoothing` to be a scalar or have shape "
320 f"({ny},)."
321 )
323 kernel = kernel.lower()
324 if kernel not in _AVAILABLE:
325 raise ValueError(f"`kernel` must be one of {_AVAILABLE}.")
327 if epsilon is None:
328 if kernel in _SCALE_INVARIANT:
329 epsilon = 1.0
330 else:
331 raise ValueError(
332 "`epsilon` must be specified if `kernel` is not one of "
333 f"{_SCALE_INVARIANT}."
334 )
335 else:
336 epsilon = float(epsilon)
338 min_degree = _NAME_TO_MIN_DEGREE.get(kernel, -1)
339 if degree is None:
340 degree = max(min_degree, 0)
341 else:
342 degree = int(degree)
343 if degree < -1:
344 raise ValueError("`degree` must be at least -1.")
345 elif degree < min_degree:
346 warnings.warn(
347 f"`degree` should not be below {min_degree} when `kernel` "
348 f"is '{kernel}'. The interpolant may not be uniquely "
349 "solvable, and the smoothing parameter may have an "
350 "unintuitive effect.",
351 UserWarning
352 )
354 if neighbors is None:
355 nobs = ny
356 else:
357 # Make sure the number of nearest neighbors used for interpolation
358 # does not exceed the number of observations.
359 neighbors = int(min(neighbors, ny))
360 nobs = neighbors
362 powers = _monomial_powers(ndim, degree)
363 # The polynomial matrix must have full column rank in order for the
364 # interpolant to be well-posed, which is not possible if there are
365 # fewer observations than monomials.
366 if powers.shape[0] > nobs:
367 raise ValueError(
368 f"At least {powers.shape[0]} data points are required when "
369 f"`degree` is {degree} and the number of dimensions is {ndim}."
370 )
372 if neighbors is None:
373 shift, scale, coeffs = _build_and_solve_system(
374 y, d, smoothing, kernel, epsilon, powers
375 )
377 # Make these attributes private since they do not always exist.
378 self._shift = shift
379 self._scale = scale
380 self._coeffs = coeffs
382 else:
383 self._tree = KDTree(y)
385 self.y = y
386 self.d = d
387 self.d_shape = d_shape
388 self.d_dtype = d_dtype
389 self.neighbors = neighbors
390 self.smoothing = smoothing
391 self.kernel = kernel
392 self.epsilon = epsilon
393 self.powers = powers
395 def _chunk_evaluator(
396 self,
397 x,
398 y,
399 shift,
400 scale,
401 coeffs,
402 memory_budget=1000000
403 ):
404 """
405 Evaluate the interpolation while controlling memory consumption.
406 We chunk the input if we need more memory than specified.
408 Parameters
409 ----------
410 x : (Q, N) float ndarray
411 array of points on which to evaluate
412 y: (P, N) float ndarray
413 array of points on which we know function values
414 shift: (N, ) ndarray
415 Domain shift used to create the polynomial matrix.
416 scale : (N,) float ndarray
417 Domain scaling used to create the polynomial matrix.
418 coeffs: (P+R, S) float ndarray
419 Coefficients in front of basis functions
420 memory_budget: int
421 Total amount of memory (in units of sizeof(float)) we wish
422 to devote for storing the array of coefficients for
423 interpolated points. If we need more memory than that, we
424 chunk the input.
426 Returns
427 -------
428 (Q, S) float ndarray
429 Interpolated array
430 """
431 nx, ndim = x.shape
432 if self.neighbors is None:
433 nnei = len(y)
434 else:
435 nnei = self.neighbors
436 # in each chunk we consume the same space we already occupy
437 chunksize = memory_budget // (self.powers.shape[0] + nnei) + 1
438 if chunksize <= nx:
439 out = np.empty((nx, self.d.shape[1]), dtype=float)
440 for i in range(0, nx, chunksize):
441 vec = _build_evaluation_coefficients(
442 x[i:i + chunksize, :],
443 y,
444 self.kernel,
445 self.epsilon,
446 self.powers,
447 shift,
448 scale)
449 out[i:i + chunksize, :] = np.dot(vec, coeffs)
450 else:
451 vec = _build_evaluation_coefficients(
452 x,
453 y,
454 self.kernel,
455 self.epsilon,
456 self.powers,
457 shift,
458 scale)
459 out = np.dot(vec, coeffs)
460 return out
462 def __call__(self, x):
463 """Evaluate the interpolant at `x`.
465 Parameters
466 ----------
467 x : (Q, N) array_like
468 Evaluation point coordinates.
470 Returns
471 -------
472 (Q, ...) ndarray
473 Values of the interpolant at `x`.
475 """
476 x = np.asarray(x, dtype=float, order="C")
477 if x.ndim != 2:
478 raise ValueError("`x` must be a 2-dimensional array.")
480 nx, ndim = x.shape
481 if ndim != self.y.shape[1]:
482 raise ValueError("Expected the second axis of `x` to have length "
483 f"{self.y.shape[1]}.")
485 # Our memory budget for storing RBF coefficients is
486 # based on how many floats in memory we already occupy
487 # If this number is below 1e6 we just use 1e6
488 # This memory budget is used to decide how we chunk
489 # the inputs
490 memory_budget = max(x.size + self.y.size + self.d.size, 1000000)
492 if self.neighbors is None:
493 out = self._chunk_evaluator(
494 x,
495 self.y,
496 self._shift,
497 self._scale,
498 self._coeffs,
499 memory_budget=memory_budget)
500 else:
501 # Get the indices of the k nearest observation points to each
502 # evaluation point.
503 _, yindices = self._tree.query(x, self.neighbors)
504 if self.neighbors == 1:
505 # `KDTree` squeezes the output when neighbors=1.
506 yindices = yindices[:, None]
508 # Multiple evaluation points may have the same neighborhood of
509 # observation points. Make the neighborhoods unique so that we only
510 # compute the interpolation coefficients once for each
511 # neighborhood.
512 yindices = np.sort(yindices, axis=1)
513 yindices, inv = np.unique(yindices, return_inverse=True, axis=0)
514 # `inv` tells us which neighborhood will be used by each evaluation
515 # point. Now we find which evaluation points will be using each
516 # neighborhood.
517 xindices = [[] for _ in range(len(yindices))]
518 for i, j in enumerate(inv):
519 xindices[j].append(i)
521 out = np.empty((nx, self.d.shape[1]), dtype=float)
522 for xidx, yidx in zip(xindices, yindices):
523 # `yidx` are the indices of the observations in this
524 # neighborhood. `xidx` are the indices of the evaluation points
525 # that are using this neighborhood.
526 xnbr = x[xidx]
527 ynbr = self.y[yidx]
528 dnbr = self.d[yidx]
529 snbr = self.smoothing[yidx]
530 shift, scale, coeffs = _build_and_solve_system(
531 ynbr,
532 dnbr,
533 snbr,
534 self.kernel,
535 self.epsilon,
536 self.powers,
537 )
538 out[xidx] = self._chunk_evaluator(
539 xnbr,
540 ynbr,
541 shift,
542 scale,
543 coeffs,
544 memory_budget=memory_budget)
546 out = out.view(self.d_dtype)
547 out = out.reshape((nx, ) + self.d_shape)
548 return out