Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_rbfinterp.py: 13%
136 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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 : (P, N) array_like
140 Data point coordinates.
141 d : (P, ...) array_like
142 Data values at `y`.
143 neighbors : int, optional
144 If specified, the value of the interpolant at each evaluation point
145 will be computed using only this many nearest data points. All the data
146 points are used by default.
147 smoothing : float or (P,) array_like, optional
148 Smoothing parameter. The interpolant perfectly fits the data when this
149 is set to 0. For large values, the interpolant approaches a least
150 squares fit of a polynomial with the specified degree. Default is 0.
151 kernel : str, optional
152 Type of RBF. This should be one of
154 - 'linear' : ``-r``
155 - 'thin_plate_spline' : ``r**2 * log(r)``
156 - 'cubic' : ``r**3``
157 - 'quintic' : ``-r**5``
158 - 'multiquadric' : ``-sqrt(1 + r**2)``
159 - 'inverse_multiquadric' : ``1/sqrt(1 + r**2)``
160 - 'inverse_quadratic' : ``1/(1 + r**2)``
161 - 'gaussian' : ``exp(-r**2)``
163 Default is 'thin_plate_spline'.
164 epsilon : float, optional
165 Shape parameter that scales the input to the RBF. If `kernel` is
166 'linear', 'thin_plate_spline', 'cubic', or 'quintic', this defaults to
167 1 and can be ignored because it has the same effect as scaling the
168 smoothing parameter. Otherwise, this must be specified.
169 degree : int, optional
170 Degree of the added polynomial. For some RBFs the interpolant may not
171 be well-posed if the polynomial degree is too small. Those RBFs and
172 their corresponding minimum degrees are
174 - 'multiquadric' : 0
175 - 'linear' : 0
176 - 'thin_plate_spline' : 1
177 - 'cubic' : 1
178 - 'quintic' : 2
180 The default value is the minimum degree for `kernel` or 0 if there is
181 no minimum degree. Set this to -1 for no added polynomial.
183 Notes
184 -----
185 An RBF is a scalar valued function in N-dimensional space whose value at
186 :math:`x` can be expressed in terms of :math:`r=||x - c||`, where :math:`c`
187 is the center of the RBF.
189 An RBF interpolant for the vector of data values :math:`d`, which are from
190 locations :math:`y`, is a linear combination of RBFs centered at :math:`y`
191 plus a polynomial with a specified degree. The RBF interpolant is written
192 as
194 .. math::
195 f(x) = K(x, y) a + P(x) b,
197 where :math:`K(x, y)` is a matrix of RBFs with centers at :math:`y`
198 evaluated at the points :math:`x`, and :math:`P(x)` is a matrix of
199 monomials, which span polynomials with the specified degree, evaluated at
200 :math:`x`. The coefficients :math:`a` and :math:`b` are the solution to the
201 linear equations
203 .. math::
204 (K(y, y) + \\lambda I) a + P(y) b = d
206 and
208 .. math::
209 P(y)^T a = 0,
211 where :math:`\\lambda` is a non-negative smoothing parameter that controls
212 how well we want to fit the data. The data are fit exactly when the
213 smoothing parameter is 0.
215 The above system is uniquely solvable if the following requirements are
216 met:
218 - :math:`P(y)` must have full column rank. :math:`P(y)` always has full
219 column rank when `degree` is -1 or 0. When `degree` is 1,
220 :math:`P(y)` has full column rank if the data point locations are not
221 all collinear (N=2), coplanar (N=3), etc.
222 - If `kernel` is 'multiquadric', 'linear', 'thin_plate_spline',
223 'cubic', or 'quintic', then `degree` must not be lower than the
224 minimum value listed above.
225 - If `smoothing` is 0, then each data point location must be distinct.
227 When using an RBF that is not scale invariant ('multiquadric',
228 'inverse_multiquadric', 'inverse_quadratic', or 'gaussian'), an appropriate
229 shape parameter must be chosen (e.g., through cross validation). Smaller
230 values for the shape parameter correspond to wider RBFs. The problem can
231 become ill-conditioned or singular when the shape parameter is too small.
233 The memory required to solve for the RBF interpolation coefficients
234 increases quadratically with the number of data points, which can become
235 impractical when interpolating more than about a thousand data points.
236 To overcome memory limitations for large interpolation problems, the
237 `neighbors` argument can be specified to compute an RBF interpolant for
238 each evaluation point using only the nearest data points.
240 .. versionadded:: 1.7.0
242 See Also
243 --------
244 NearestNDInterpolator
245 LinearNDInterpolator
246 CloughTocher2DInterpolator
248 References
249 ----------
250 .. [1] Fasshauer, G., 2007. Meshfree Approximation Methods with Matlab.
251 World Scientific Publishing Co.
253 .. [2] http://amadeus.math.iit.edu/~fass/603_ch3.pdf
255 .. [3] Wahba, G., 1990. Spline Models for Observational Data. SIAM.
257 .. [4] http://pages.stat.wisc.edu/~wahba/stat860public/lect/lect8/lect8.pdf
259 Examples
260 --------
261 Demonstrate interpolating scattered data to a grid in 2-D.
263 >>> import numpy as np
264 >>> import matplotlib.pyplot as plt
265 >>> from scipy.interpolate import RBFInterpolator
266 >>> from scipy.stats.qmc import Halton
268 >>> rng = np.random.default_rng()
269 >>> xobs = 2*Halton(2, seed=rng).random(100) - 1
270 >>> yobs = np.sum(xobs, axis=1)*np.exp(-6*np.sum(xobs**2, axis=1))
272 >>> xgrid = np.mgrid[-1:1:50j, -1:1:50j]
273 >>> xflat = xgrid.reshape(2, -1).T
274 >>> yflat = RBFInterpolator(xobs, yobs)(xflat)
275 >>> ygrid = yflat.reshape(50, 50)
277 >>> fig, ax = plt.subplots()
278 >>> ax.pcolormesh(*xgrid, ygrid, vmin=-0.25, vmax=0.25, shading='gouraud')
279 >>> p = ax.scatter(*xobs.T, c=yobs, s=50, ec='k', vmin=-0.25, vmax=0.25)
280 >>> fig.colorbar(p)
281 >>> plt.show()
283 """
285 def __init__(self, y, d,
286 neighbors=None,
287 smoothing=0.0,
288 kernel="thin_plate_spline",
289 epsilon=None,
290 degree=None):
291 y = np.asarray(y, dtype=float, order="C")
292 if y.ndim != 2:
293 raise ValueError("`y` must be a 2-dimensional array.")
295 ny, ndim = y.shape
297 d_dtype = complex if np.iscomplexobj(d) else float
298 d = np.asarray(d, dtype=d_dtype, order="C")
299 if d.shape[0] != ny:
300 raise ValueError(
301 f"Expected the first axis of `d` to have length {ny}."
302 )
304 d_shape = d.shape[1:]
305 d = d.reshape((ny, -1))
306 # If `d` is complex, convert it to a float array with twice as many
307 # columns. Otherwise, the LHS matrix would need to be converted to
308 # complex and take up 2x more memory than necessary.
309 d = d.view(float)
311 if np.isscalar(smoothing):
312 smoothing = np.full(ny, smoothing, dtype=float)
313 else:
314 smoothing = np.asarray(smoothing, dtype=float, order="C")
315 if smoothing.shape != (ny,):
316 raise ValueError(
317 "Expected `smoothing` to be a scalar or have shape "
318 f"({ny},)."
319 )
321 kernel = kernel.lower()
322 if kernel not in _AVAILABLE:
323 raise ValueError(f"`kernel` must be one of {_AVAILABLE}.")
325 if epsilon is None:
326 if kernel in _SCALE_INVARIANT:
327 epsilon = 1.0
328 else:
329 raise ValueError(
330 "`epsilon` must be specified if `kernel` is not one of "
331 f"{_SCALE_INVARIANT}."
332 )
333 else:
334 epsilon = float(epsilon)
336 min_degree = _NAME_TO_MIN_DEGREE.get(kernel, -1)
337 if degree is None:
338 degree = max(min_degree, 0)
339 else:
340 degree = int(degree)
341 if degree < -1:
342 raise ValueError("`degree` must be at least -1.")
343 elif degree < min_degree:
344 warnings.warn(
345 f"`degree` should not be below {min_degree} when `kernel` "
346 f"is '{kernel}'. The interpolant may not be uniquely "
347 "solvable, and the smoothing parameter may have an "
348 "unintuitive effect.",
349 UserWarning
350 )
352 if neighbors is None:
353 nobs = ny
354 else:
355 # Make sure the number of nearest neighbors used for interpolation
356 # does not exceed the number of observations.
357 neighbors = int(min(neighbors, ny))
358 nobs = neighbors
360 powers = _monomial_powers(ndim, degree)
361 # The polynomial matrix must have full column rank in order for the
362 # interpolant to be well-posed, which is not possible if there are
363 # fewer observations than monomials.
364 if powers.shape[0] > nobs:
365 raise ValueError(
366 f"At least {powers.shape[0]} data points are required when "
367 f"`degree` is {degree} and the number of dimensions is {ndim}."
368 )
370 if neighbors is None:
371 shift, scale, coeffs = _build_and_solve_system(
372 y, d, smoothing, kernel, epsilon, powers
373 )
375 # Make these attributes private since they do not always exist.
376 self._shift = shift
377 self._scale = scale
378 self._coeffs = coeffs
380 else:
381 self._tree = KDTree(y)
383 self.y = y
384 self.d = d
385 self.d_shape = d_shape
386 self.d_dtype = d_dtype
387 self.neighbors = neighbors
388 self.smoothing = smoothing
389 self.kernel = kernel
390 self.epsilon = epsilon
391 self.powers = powers
393 def _chunk_evaluator(
394 self,
395 x,
396 y,
397 shift,
398 scale,
399 coeffs,
400 memory_budget=1000000
401 ):
402 """
403 Evaluate the interpolation while controlling memory consumption.
404 We chunk the input if we need more memory than specified.
406 Parameters
407 ----------
408 x : (Q, N) float ndarray
409 array of points on which to evaluate
410 y: (P, N) float ndarray
411 array of points on which we know function values
412 shift: (N, ) ndarray
413 Domain shift used to create the polynomial matrix.
414 scale : (N,) float ndarray
415 Domain scaling used to create the polynomial matrix.
416 coeffs: (P+R, S) float ndarray
417 Coefficients in front of basis functions
418 memory_budget: int
419 Total amount of memory (in units of sizeof(float)) we wish
420 to devote for storing the array of coefficients for
421 interpolated points. If we need more memory than that, we
422 chunk the input.
424 Returns
425 -------
426 (Q, S) float ndarray
427 Interpolated array
428 """
429 nx, ndim = x.shape
430 if self.neighbors is None:
431 nnei = len(y)
432 else:
433 nnei = self.neighbors
434 # in each chunk we consume the same space we already occupy
435 chunksize = memory_budget // ((self.powers.shape[0] + nnei)) + 1
436 if chunksize <= nx:
437 out = np.empty((nx, self.d.shape[1]), dtype=float)
438 for i in range(0, nx, chunksize):
439 vec = _build_evaluation_coefficients(
440 x[i:i + chunksize, :],
441 y,
442 self.kernel,
443 self.epsilon,
444 self.powers,
445 shift,
446 scale)
447 out[i:i + chunksize, :] = np.dot(vec, coeffs)
448 else:
449 vec = _build_evaluation_coefficients(
450 x,
451 y,
452 self.kernel,
453 self.epsilon,
454 self.powers,
455 shift,
456 scale)
457 out = np.dot(vec, coeffs)
458 return out
460 def __call__(self, x):
461 """Evaluate the interpolant at `x`.
463 Parameters
464 ----------
465 x : (Q, N) array_like
466 Evaluation point coordinates.
468 Returns
469 -------
470 (Q, ...) ndarray
471 Values of the interpolant at `x`.
473 """
474 x = np.asarray(x, dtype=float, order="C")
475 if x.ndim != 2:
476 raise ValueError("`x` must be a 2-dimensional array.")
478 nx, ndim = x.shape
479 if ndim != self.y.shape[1]:
480 raise ValueError("Expected the second axis of `x` to have length "
481 f"{self.y.shape[1]}.")
483 # Our memory budget for storing RBF coefficients is
484 # based on how many floats in memory we already occupy
485 # If this number is below 1e6 we just use 1e6
486 # This memory budget is used to decide how we chunk
487 # the inputs
488 memory_budget = max(x.size + self.y.size + self.d.size, 1000000)
490 if self.neighbors is None:
491 out = self._chunk_evaluator(
492 x,
493 self.y,
494 self._shift,
495 self._scale,
496 self._coeffs,
497 memory_budget=memory_budget)
498 else:
499 # Get the indices of the k nearest observation points to each
500 # evaluation point.
501 _, yindices = self._tree.query(x, self.neighbors)
502 if self.neighbors == 1:
503 # `KDTree` squeezes the output when neighbors=1.
504 yindices = yindices[:, None]
506 # Multiple evaluation points may have the same neighborhood of
507 # observation points. Make the neighborhoods unique so that we only
508 # compute the interpolation coefficients once for each
509 # neighborhood.
510 yindices = np.sort(yindices, axis=1)
511 yindices, inv = np.unique(yindices, return_inverse=True, axis=0)
512 # `inv` tells us which neighborhood will be used by each evaluation
513 # point. Now we find which evaluation points will be using each
514 # neighborhood.
515 xindices = [[] for _ in range(len(yindices))]
516 for i, j in enumerate(inv):
517 xindices[j].append(i)
519 out = np.empty((nx, self.d.shape[1]), dtype=float)
520 for xidx, yidx in zip(xindices, yindices):
521 # `yidx` are the indices of the observations in this
522 # neighborhood. `xidx` are the indices of the evaluation points
523 # that are using this neighborhood.
524 xnbr = x[xidx]
525 ynbr = self.y[yidx]
526 dnbr = self.d[yidx]
527 snbr = self.smoothing[yidx]
528 shift, scale, coeffs = _build_and_solve_system(
529 ynbr,
530 dnbr,
531 snbr,
532 self.kernel,
533 self.epsilon,
534 self.powers,
535 )
536 out[xidx] = self._chunk_evaluator(
537 xnbr,
538 ynbr,
539 shift,
540 scale,
541 coeffs,
542 memory_budget=memory_budget)
544 out = out.view(self.d_dtype)
545 out = out.reshape((nx, ) + self.d_shape)
546 return out