Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_cubic.py: 10%
262 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Interpolation algorithms using piecewise cubic polynomials."""
3import numpy as np
5from . import PPoly
6from ._polyint import _isscalar
7from scipy.linalg import solve_banded, solve
10__all__ = ["CubicHermiteSpline", "PchipInterpolator", "pchip_interpolate",
11 "Akima1DInterpolator", "CubicSpline"]
14def prepare_input(x, y, axis, dydx=None):
15 """Prepare input for cubic spline interpolators.
17 All data are converted to numpy arrays and checked for correctness.
18 Axes equal to `axis` of arrays `y` and `dydx` are moved to be the 0th
19 axis. The value of `axis` is converted to lie in
20 [0, number of dimensions of `y`).
21 """
23 x, y = map(np.asarray, (x, y))
24 if np.issubdtype(x.dtype, np.complexfloating):
25 raise ValueError("`x` must contain real values.")
26 x = x.astype(float)
28 if np.issubdtype(y.dtype, np.complexfloating):
29 dtype = complex
30 else:
31 dtype = float
33 if dydx is not None:
34 dydx = np.asarray(dydx)
35 if y.shape != dydx.shape:
36 raise ValueError("The shapes of `y` and `dydx` must be identical.")
37 if np.issubdtype(dydx.dtype, np.complexfloating):
38 dtype = complex
39 dydx = dydx.astype(dtype, copy=False)
41 y = y.astype(dtype, copy=False)
42 axis = axis % y.ndim
43 if x.ndim != 1:
44 raise ValueError("`x` must be 1-dimensional.")
45 if x.shape[0] < 2:
46 raise ValueError("`x` must contain at least 2 elements.")
47 if x.shape[0] != y.shape[axis]:
48 raise ValueError("The length of `y` along `axis`={0} doesn't "
49 "match the length of `x`".format(axis))
51 if not np.all(np.isfinite(x)):
52 raise ValueError("`x` must contain only finite values.")
53 if not np.all(np.isfinite(y)):
54 raise ValueError("`y` must contain only finite values.")
56 if dydx is not None and not np.all(np.isfinite(dydx)):
57 raise ValueError("`dydx` must contain only finite values.")
59 dx = np.diff(x)
60 if np.any(dx <= 0):
61 raise ValueError("`x` must be strictly increasing sequence.")
63 y = np.moveaxis(y, axis, 0)
64 if dydx is not None:
65 dydx = np.moveaxis(dydx, axis, 0)
67 return x, dx, y, axis, dydx
70class CubicHermiteSpline(PPoly):
71 """Piecewise-cubic interpolator matching values and first derivatives.
73 The result is represented as a `PPoly` instance.
75 Parameters
76 ----------
77 x : array_like, shape (n,)
78 1-D array containing values of the independent variable.
79 Values must be real, finite and in strictly increasing order.
80 y : array_like
81 Array containing values of the dependent variable. It can have
82 arbitrary number of dimensions, but the length along ``axis``
83 (see below) must match the length of ``x``. Values must be finite.
84 dydx : array_like
85 Array containing derivatives of the dependent variable. It can have
86 arbitrary number of dimensions, but the length along ``axis``
87 (see below) must match the length of ``x``. Values must be finite.
88 axis : int, optional
89 Axis along which `y` is assumed to be varying. Meaning that for
90 ``x[i]`` the corresponding values are ``np.take(y, i, axis=axis)``.
91 Default is 0.
92 extrapolate : {bool, 'periodic', None}, optional
93 If bool, determines whether to extrapolate to out-of-bounds points
94 based on first and last intervals, or to return NaNs. If 'periodic',
95 periodic extrapolation is used. If None (default), it is set to True.
97 Attributes
98 ----------
99 x : ndarray, shape (n,)
100 Breakpoints. The same ``x`` which was passed to the constructor.
101 c : ndarray, shape (4, n-1, ...)
102 Coefficients of the polynomials on each segment. The trailing
103 dimensions match the dimensions of `y`, excluding ``axis``.
104 For example, if `y` is 1-D, then ``c[k, i]`` is a coefficient for
105 ``(x-x[i])**(3-k)`` on the segment between ``x[i]`` and ``x[i+1]``.
106 axis : int
107 Interpolation axis. The same axis which was passed to the
108 constructor.
110 Methods
111 -------
112 __call__
113 derivative
114 antiderivative
115 integrate
116 roots
118 See Also
119 --------
120 Akima1DInterpolator : Akima 1D interpolator.
121 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
122 CubicSpline : Cubic spline data interpolator.
123 PPoly : Piecewise polynomial in terms of coefficients and breakpoints
125 Notes
126 -----
127 If you want to create a higher-order spline matching higher-order
128 derivatives, use `BPoly.from_derivatives`.
130 References
131 ----------
132 .. [1] `Cubic Hermite spline
133 <https://en.wikipedia.org/wiki/Cubic_Hermite_spline>`_
134 on Wikipedia.
135 """
137 def __init__(self, x, y, dydx, axis=0, extrapolate=None):
138 if extrapolate is None:
139 extrapolate = True
141 x, dx, y, axis, dydx = prepare_input(x, y, axis, dydx)
143 dxr = dx.reshape([dx.shape[0]] + [1] * (y.ndim - 1))
144 slope = np.diff(y, axis=0) / dxr
145 t = (dydx[:-1] + dydx[1:] - 2 * slope) / dxr
147 c = np.empty((4, len(x) - 1) + y.shape[1:], dtype=t.dtype)
148 c[0] = t / dxr
149 c[1] = (slope - dydx[:-1]) / dxr - t
150 c[2] = dydx[:-1]
151 c[3] = y[:-1]
153 super().__init__(c, x, extrapolate=extrapolate)
154 self.axis = axis
157class PchipInterpolator(CubicHermiteSpline):
158 r"""PCHIP 1-D monotonic cubic interpolation.
160 ``x`` and ``y`` are arrays of values used to approximate some function f,
161 with ``y = f(x)``. The interpolant uses monotonic cubic splines
162 to find the value of new points. (PCHIP stands for Piecewise Cubic
163 Hermite Interpolating Polynomial).
165 Parameters
166 ----------
167 x : ndarray
168 A 1-D array of monotonically increasing real values. ``x`` cannot
169 include duplicate values (otherwise f is overspecified)
170 y : ndarray
171 A 1-D array of real values. ``y``'s length along the interpolation
172 axis must be equal to the length of ``x``. If N-D array, use ``axis``
173 parameter to select correct axis.
174 axis : int, optional
175 Axis in the y array corresponding to the x-coordinate values.
176 extrapolate : bool, optional
177 Whether to extrapolate to out-of-bounds points based on first
178 and last intervals, or to return NaNs.
180 Methods
181 -------
182 __call__
183 derivative
184 antiderivative
185 roots
187 See Also
188 --------
189 CubicHermiteSpline : Piecewise-cubic interpolator.
190 Akima1DInterpolator : Akima 1D interpolator.
191 CubicSpline : Cubic spline data interpolator.
192 PPoly : Piecewise polynomial in terms of coefficients and breakpoints.
194 Notes
195 -----
196 The interpolator preserves monotonicity in the interpolation data and does
197 not overshoot if the data is not smooth.
199 The first derivatives are guaranteed to be continuous, but the second
200 derivatives may jump at :math:`x_k`.
202 Determines the derivatives at the points :math:`x_k`, :math:`f'_k`,
203 by using PCHIP algorithm [1]_.
205 Let :math:`h_k = x_{k+1} - x_k`, and :math:`d_k = (y_{k+1} - y_k) / h_k`
206 are the slopes at internal points :math:`x_k`.
207 If the signs of :math:`d_k` and :math:`d_{k-1}` are different or either of
208 them equals zero, then :math:`f'_k = 0`. Otherwise, it is given by the
209 weighted harmonic mean
211 .. math::
213 \frac{w_1 + w_2}{f'_k} = \frac{w_1}{d_{k-1}} + \frac{w_2}{d_k}
215 where :math:`w_1 = 2 h_k + h_{k-1}` and :math:`w_2 = h_k + 2 h_{k-1}`.
217 The end slopes are set using a one-sided scheme [2]_.
220 References
221 ----------
222 .. [1] F. N. Fritsch and J. Butland,
223 A method for constructing local
224 monotone piecewise cubic interpolants,
225 SIAM J. Sci. Comput., 5(2), 300-304 (1984).
226 :doi:`10.1137/0905021`.
227 .. [2] see, e.g., C. Moler, Numerical Computing with Matlab, 2004.
228 :doi:`10.1137/1.9780898717952`
231 """
233 def __init__(self, x, y, axis=0, extrapolate=None):
234 x, _, y, axis, _ = prepare_input(x, y, axis)
235 xp = x.reshape((x.shape[0],) + (1,)*(y.ndim-1))
236 dk = self._find_derivatives(xp, y)
237 super().__init__(x, y, dk, axis=0, extrapolate=extrapolate)
238 self.axis = axis
240 @staticmethod
241 def _edge_case(h0, h1, m0, m1):
242 # one-sided three-point estimate for the derivative
243 d = ((2*h0 + h1)*m0 - h0*m1) / (h0 + h1)
245 # try to preserve shape
246 mask = np.sign(d) != np.sign(m0)
247 mask2 = (np.sign(m0) != np.sign(m1)) & (np.abs(d) > 3.*np.abs(m0))
248 mmm = (~mask) & mask2
250 d[mask] = 0.
251 d[mmm] = 3.*m0[mmm]
253 return d
255 @staticmethod
256 def _find_derivatives(x, y):
257 # Determine the derivatives at the points y_k, d_k, by using
258 # PCHIP algorithm is:
259 # We choose the derivatives at the point x_k by
260 # Let m_k be the slope of the kth segment (between k and k+1)
261 # If m_k=0 or m_{k-1}=0 or sgn(m_k) != sgn(m_{k-1}) then d_k == 0
262 # else use weighted harmonic mean:
263 # w_1 = 2h_k + h_{k-1}, w_2 = h_k + 2h_{k-1}
264 # 1/d_k = 1/(w_1 + w_2)*(w_1 / m_k + w_2 / m_{k-1})
265 # where h_k is the spacing between x_k and x_{k+1}
266 y_shape = y.shape
267 if y.ndim == 1:
268 # So that _edge_case doesn't end up assigning to scalars
269 x = x[:, None]
270 y = y[:, None]
272 hk = x[1:] - x[:-1]
273 mk = (y[1:] - y[:-1]) / hk
275 if y.shape[0] == 2:
276 # edge case: only have two points, use linear interpolation
277 dk = np.zeros_like(y)
278 dk[0] = mk
279 dk[1] = mk
280 return dk.reshape(y_shape)
282 smk = np.sign(mk)
283 condition = (smk[1:] != smk[:-1]) | (mk[1:] == 0) | (mk[:-1] == 0)
285 w1 = 2*hk[1:] + hk[:-1]
286 w2 = hk[1:] + 2*hk[:-1]
288 # values where division by zero occurs will be excluded
289 # by 'condition' afterwards
290 with np.errstate(divide='ignore', invalid='ignore'):
291 whmean = (w1/mk[:-1] + w2/mk[1:]) / (w1 + w2)
293 dk = np.zeros_like(y)
294 dk[1:-1][condition] = 0.0
295 dk[1:-1][~condition] = 1.0 / whmean[~condition]
297 # special case endpoints, as suggested in
298 # Cleve Moler, Numerical Computing with MATLAB, Chap 3.6 (pchiptx.m)
299 dk[0] = PchipInterpolator._edge_case(hk[0], hk[1], mk[0], mk[1])
300 dk[-1] = PchipInterpolator._edge_case(hk[-1], hk[-2], mk[-1], mk[-2])
302 return dk.reshape(y_shape)
305def pchip_interpolate(xi, yi, x, der=0, axis=0):
306 """
307 Convenience function for pchip interpolation.
309 xi and yi are arrays of values used to approximate some function f,
310 with ``yi = f(xi)``. The interpolant uses monotonic cubic splines
311 to find the value of new points x and the derivatives there.
313 See `scipy.interpolate.PchipInterpolator` for details.
315 Parameters
316 ----------
317 xi : array_like
318 A sorted list of x-coordinates, of length N.
319 yi : array_like
320 A 1-D array of real values. `yi`'s length along the interpolation
321 axis must be equal to the length of `xi`. If N-D array, use axis
322 parameter to select correct axis.
323 x : scalar or array_like
324 Of length M.
325 der : int or list, optional
326 Derivatives to extract. The 0th derivative can be included to
327 return the function value.
328 axis : int, optional
329 Axis in the yi array corresponding to the x-coordinate values.
331 See Also
332 --------
333 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
335 Returns
336 -------
337 y : scalar or array_like
338 The result, of length R or length M or M by R,
340 Examples
341 --------
342 We can interpolate 2D observed data using pchip interpolation:
344 >>> import numpy as np
345 >>> import matplotlib.pyplot as plt
346 >>> from scipy.interpolate import pchip_interpolate
347 >>> x_observed = np.linspace(0.0, 10.0, 11)
348 >>> y_observed = np.sin(x_observed)
349 >>> x = np.linspace(min(x_observed), max(x_observed), num=100)
350 >>> y = pchip_interpolate(x_observed, y_observed, x)
351 >>> plt.plot(x_observed, y_observed, "o", label="observation")
352 >>> plt.plot(x, y, label="pchip interpolation")
353 >>> plt.legend()
354 >>> plt.show()
356 """
357 P = PchipInterpolator(xi, yi, axis=axis)
359 if der == 0:
360 return P(x)
361 elif _isscalar(der):
362 return P.derivative(der)(x)
363 else:
364 return [P.derivative(nu)(x) for nu in der]
367class Akima1DInterpolator(CubicHermiteSpline):
368 """
369 Akima interpolator
371 Fit piecewise cubic polynomials, given vectors x and y. The interpolation
372 method by Akima uses a continuously differentiable sub-spline built from
373 piecewise cubic polynomials. The resultant curve passes through the given
374 data points and will appear smooth and natural.
376 Parameters
377 ----------
378 x : ndarray, shape (m, )
379 1-D array of monotonically increasing real values.
380 y : ndarray, shape (m, ...)
381 N-D array of real values. The length of ``y`` along the first axis
382 must be equal to the length of ``x``.
383 axis : int, optional
384 Specifies the axis of ``y`` along which to interpolate. Interpolation
385 defaults to the first axis of ``y``.
387 Methods
388 -------
389 __call__
390 derivative
391 antiderivative
392 roots
394 See Also
395 --------
396 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
397 CubicSpline : Cubic spline data interpolator.
398 PPoly : Piecewise polynomial in terms of coefficients and breakpoints
400 Notes
401 -----
402 .. versionadded:: 0.14
404 Use only for precise data, as the fitted curve passes through the given
405 points exactly. This routine is useful for plotting a pleasingly smooth
406 curve through a few given points for purposes of plotting.
408 References
409 ----------
410 [1] A new method of interpolation and smooth curve fitting based
411 on local procedures. Hiroshi Akima, J. ACM, October 1970, 17(4),
412 589-602.
414 """
416 def __init__(self, x, y, axis=0):
417 # Original implementation in MATLAB by N. Shamsundar (BSD licensed), see
418 # https://www.mathworks.com/matlabcentral/fileexchange/1814-akima-interpolation
419 x, dx, y, axis, _ = prepare_input(x, y, axis)
421 # determine slopes between breakpoints
422 m = np.empty((x.size + 3, ) + y.shape[1:])
423 dx = dx[(slice(None), ) + (None, ) * (y.ndim - 1)]
424 m[2:-2] = np.diff(y, axis=0) / dx
426 # add two additional points on the left ...
427 m[1] = 2. * m[2] - m[3]
428 m[0] = 2. * m[1] - m[2]
429 # ... and on the right
430 m[-2] = 2. * m[-3] - m[-4]
431 m[-1] = 2. * m[-2] - m[-3]
433 # if m1 == m2 != m3 == m4, the slope at the breakpoint is not
434 # defined. This is the fill value:
435 t = .5 * (m[3:] + m[:-3])
436 # get the denominator of the slope t
437 dm = np.abs(np.diff(m, axis=0))
438 f1 = dm[2:]
439 f2 = dm[:-2]
440 f12 = f1 + f2
441 # These are the mask of where the slope at breakpoint is defined:
442 ind = np.nonzero(f12 > 1e-9 * np.max(f12, initial=-np.inf))
443 x_ind, y_ind = ind[0], ind[1:]
444 # Set the slope at breakpoint
445 t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
446 f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
448 super().__init__(x, y, t, axis=0, extrapolate=False)
449 self.axis = axis
451 def extend(self, c, x, right=True):
452 raise NotImplementedError("Extending a 1-D Akima interpolator is not "
453 "yet implemented")
455 # These are inherited from PPoly, but they do not produce an Akima
456 # interpolator. Hence stub them out.
457 @classmethod
458 def from_spline(cls, tck, extrapolate=None):
459 raise NotImplementedError("This method does not make sense for "
460 "an Akima interpolator.")
462 @classmethod
463 def from_bernstein_basis(cls, bp, extrapolate=None):
464 raise NotImplementedError("This method does not make sense for "
465 "an Akima interpolator.")
468class CubicSpline(CubicHermiteSpline):
469 """Cubic spline data interpolator.
471 Interpolate data with a piecewise cubic polynomial which is twice
472 continuously differentiable [1]_. The result is represented as a `PPoly`
473 instance with breakpoints matching the given data.
475 Parameters
476 ----------
477 x : array_like, shape (n,)
478 1-D array containing values of the independent variable.
479 Values must be real, finite and in strictly increasing order.
480 y : array_like
481 Array containing values of the dependent variable. It can have
482 arbitrary number of dimensions, but the length along ``axis``
483 (see below) must match the length of ``x``. Values must be finite.
484 axis : int, optional
485 Axis along which `y` is assumed to be varying. Meaning that for
486 ``x[i]`` the corresponding values are ``np.take(y, i, axis=axis)``.
487 Default is 0.
488 bc_type : string or 2-tuple, optional
489 Boundary condition type. Two additional equations, given by the
490 boundary conditions, are required to determine all coefficients of
491 polynomials on each segment [2]_.
493 If `bc_type` is a string, then the specified condition will be applied
494 at both ends of a spline. Available conditions are:
496 * 'not-a-knot' (default): The first and second segment at a curve end
497 are the same polynomial. It is a good default when there is no
498 information on boundary conditions.
499 * 'periodic': The interpolated functions is assumed to be periodic
500 of period ``x[-1] - x[0]``. The first and last value of `y` must be
501 identical: ``y[0] == y[-1]``. This boundary condition will result in
502 ``y'[0] == y'[-1]`` and ``y''[0] == y''[-1]``.
503 * 'clamped': The first derivative at curves ends are zero. Assuming
504 a 1D `y`, ``bc_type=((1, 0.0), (1, 0.0))`` is the same condition.
505 * 'natural': The second derivative at curve ends are zero. Assuming
506 a 1D `y`, ``bc_type=((2, 0.0), (2, 0.0))`` is the same condition.
508 If `bc_type` is a 2-tuple, the first and the second value will be
509 applied at the curve start and end respectively. The tuple values can
510 be one of the previously mentioned strings (except 'periodic') or a
511 tuple `(order, deriv_values)` allowing to specify arbitrary
512 derivatives at curve ends:
514 * `order`: the derivative order, 1 or 2.
515 * `deriv_value`: array_like containing derivative values, shape must
516 be the same as `y`, excluding ``axis`` dimension. For example, if
517 `y` is 1-D, then `deriv_value` must be a scalar. If `y` is 3-D with
518 the shape (n0, n1, n2) and axis=2, then `deriv_value` must be 2-D
519 and have the shape (n0, n1).
520 extrapolate : {bool, 'periodic', None}, optional
521 If bool, determines whether to extrapolate to out-of-bounds points
522 based on first and last intervals, or to return NaNs. If 'periodic',
523 periodic extrapolation is used. If None (default), ``extrapolate`` is
524 set to 'periodic' for ``bc_type='periodic'`` and to True otherwise.
526 Attributes
527 ----------
528 x : ndarray, shape (n,)
529 Breakpoints. The same ``x`` which was passed to the constructor.
530 c : ndarray, shape (4, n-1, ...)
531 Coefficients of the polynomials on each segment. The trailing
532 dimensions match the dimensions of `y`, excluding ``axis``.
533 For example, if `y` is 1-d, then ``c[k, i]`` is a coefficient for
534 ``(x-x[i])**(3-k)`` on the segment between ``x[i]`` and ``x[i+1]``.
535 axis : int
536 Interpolation axis. The same axis which was passed to the
537 constructor.
539 Methods
540 -------
541 __call__
542 derivative
543 antiderivative
544 integrate
545 roots
547 See Also
548 --------
549 Akima1DInterpolator : Akima 1D interpolator.
550 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
551 PPoly : Piecewise polynomial in terms of coefficients and breakpoints.
553 Notes
554 -----
555 Parameters `bc_type` and ``extrapolate`` work independently, i.e. the
556 former controls only construction of a spline, and the latter only
557 evaluation.
559 When a boundary condition is 'not-a-knot' and n = 2, it is replaced by
560 a condition that the first derivative is equal to the linear interpolant
561 slope. When both boundary conditions are 'not-a-knot' and n = 3, the
562 solution is sought as a parabola passing through given points.
564 When 'not-a-knot' boundary conditions is applied to both ends, the
565 resulting spline will be the same as returned by `splrep` (with ``s=0``)
566 and `InterpolatedUnivariateSpline`, but these two methods use a
567 representation in B-spline basis.
569 .. versionadded:: 0.18.0
571 Examples
572 --------
573 In this example the cubic spline is used to interpolate a sampled sinusoid.
574 You can see that the spline continuity property holds for the first and
575 second derivatives and violates only for the third derivative.
577 >>> import numpy as np
578 >>> from scipy.interpolate import CubicSpline
579 >>> import matplotlib.pyplot as plt
580 >>> x = np.arange(10)
581 >>> y = np.sin(x)
582 >>> cs = CubicSpline(x, y)
583 >>> xs = np.arange(-0.5, 9.6, 0.1)
584 >>> fig, ax = plt.subplots(figsize=(6.5, 4))
585 >>> ax.plot(x, y, 'o', label='data')
586 >>> ax.plot(xs, np.sin(xs), label='true')
587 >>> ax.plot(xs, cs(xs), label="S")
588 >>> ax.plot(xs, cs(xs, 1), label="S'")
589 >>> ax.plot(xs, cs(xs, 2), label="S''")
590 >>> ax.plot(xs, cs(xs, 3), label="S'''")
591 >>> ax.set_xlim(-0.5, 9.5)
592 >>> ax.legend(loc='lower left', ncol=2)
593 >>> plt.show()
595 In the second example, the unit circle is interpolated with a spline. A
596 periodic boundary condition is used. You can see that the first derivative
597 values, ds/dx=0, ds/dy=1 at the periodic point (1, 0) are correctly
598 computed. Note that a circle cannot be exactly represented by a cubic
599 spline. To increase precision, more breakpoints would be required.
601 >>> theta = 2 * np.pi * np.linspace(0, 1, 5)
602 >>> y = np.c_[np.cos(theta), np.sin(theta)]
603 >>> cs = CubicSpline(theta, y, bc_type='periodic')
604 >>> print("ds/dx={:.1f} ds/dy={:.1f}".format(cs(0, 1)[0], cs(0, 1)[1]))
605 ds/dx=0.0 ds/dy=1.0
606 >>> xs = 2 * np.pi * np.linspace(0, 1, 100)
607 >>> fig, ax = plt.subplots(figsize=(6.5, 4))
608 >>> ax.plot(y[:, 0], y[:, 1], 'o', label='data')
609 >>> ax.plot(np.cos(xs), np.sin(xs), label='true')
610 >>> ax.plot(cs(xs)[:, 0], cs(xs)[:, 1], label='spline')
611 >>> ax.axes.set_aspect('equal')
612 >>> ax.legend(loc='center')
613 >>> plt.show()
615 The third example is the interpolation of a polynomial y = x**3 on the
616 interval 0 <= x<= 1. A cubic spline can represent this function exactly.
617 To achieve that we need to specify values and first derivatives at
618 endpoints of the interval. Note that y' = 3 * x**2 and thus y'(0) = 0 and
619 y'(1) = 3.
621 >>> cs = CubicSpline([0, 1], [0, 1], bc_type=((1, 0), (1, 3)))
622 >>> x = np.linspace(0, 1)
623 >>> np.allclose(x**3, cs(x))
624 True
626 References
627 ----------
628 .. [1] `Cubic Spline Interpolation
629 <https://en.wikiversity.org/wiki/Cubic_Spline_Interpolation>`_
630 on Wikiversity.
631 .. [2] Carl de Boor, "A Practical Guide to Splines", Springer-Verlag, 1978.
632 """
634 def __init__(self, x, y, axis=0, bc_type='not-a-knot', extrapolate=None):
635 x, dx, y, axis, _ = prepare_input(x, y, axis)
636 n = len(x)
638 bc, y = self._validate_bc(bc_type, y, y.shape[1:], axis)
640 if extrapolate is None:
641 if bc[0] == 'periodic':
642 extrapolate = 'periodic'
643 else:
644 extrapolate = True
646 if y.size == 0:
647 # bail out early for zero-sized arrays
648 s = np.zeros_like(y)
649 else:
650 dxr = dx.reshape([dx.shape[0]] + [1] * (y.ndim - 1))
651 slope = np.diff(y, axis=0) / dxr
653 # If bc is 'not-a-knot' this change is just a convention.
654 # If bc is 'periodic' then we already checked that y[0] == y[-1],
655 # and the spline is just a constant, we handle this case in the
656 # same way by setting the first derivatives to slope, which is 0.
657 if n == 2:
658 if bc[0] in ['not-a-knot', 'periodic']:
659 bc[0] = (1, slope[0])
660 if bc[1] in ['not-a-knot', 'periodic']:
661 bc[1] = (1, slope[0])
663 # This is a special case, when both conditions are 'not-a-knot'
664 # and n == 3. In this case 'not-a-knot' can't be handled regularly
665 # as the both conditions are identical. We handle this case by
666 # constructing a parabola passing through given points.
667 if n == 3 and bc[0] == 'not-a-knot' and bc[1] == 'not-a-knot':
668 A = np.zeros((3, 3)) # This is a standard matrix.
669 b = np.empty((3,) + y.shape[1:], dtype=y.dtype)
671 A[0, 0] = 1
672 A[0, 1] = 1
673 A[1, 0] = dx[1]
674 A[1, 1] = 2 * (dx[0] + dx[1])
675 A[1, 2] = dx[0]
676 A[2, 1] = 1
677 A[2, 2] = 1
679 b[0] = 2 * slope[0]
680 b[1] = 3 * (dxr[0] * slope[1] + dxr[1] * slope[0])
681 b[2] = 2 * slope[1]
683 s = solve(A, b, overwrite_a=True, overwrite_b=True,
684 check_finite=False)
685 elif n == 3 and bc[0] == 'periodic':
686 # In case when number of points is 3 we compute the derivatives
687 # manually
688 s = np.empty((n,) + y.shape[1:], dtype=y.dtype)
689 t = (slope / dxr).sum() / (1. / dxr).sum()
690 s.fill(t)
691 else:
692 # Find derivative values at each x[i] by solving a tridiagonal
693 # system.
694 A = np.zeros((3, n)) # This is a banded matrix representation.
695 b = np.empty((n,) + y.shape[1:], dtype=y.dtype)
697 # Filling the system for i=1..n-2
698 # (x[i-1] - x[i]) * s[i-1] +\
699 # 2 * ((x[i] - x[i-1]) + (x[i+1] - x[i])) * s[i] +\
700 # (x[i] - x[i-1]) * s[i+1] =\
701 # 3 * ((x[i+1] - x[i])*(y[i] - y[i-1])/(x[i] - x[i-1]) +\
702 # (x[i] - x[i-1])*(y[i+1] - y[i])/(x[i+1] - x[i]))
704 A[1, 1:-1] = 2 * (dx[:-1] + dx[1:]) # The diagonal
705 A[0, 2:] = dx[:-1] # The upper diagonal
706 A[-1, :-2] = dx[1:] # The lower diagonal
708 b[1:-1] = 3 * (dxr[1:] * slope[:-1] + dxr[:-1] * slope[1:])
710 bc_start, bc_end = bc
712 if bc_start == 'periodic':
713 # Due to the periodicity, and because y[-1] = y[0], the
714 # linear system has (n-1) unknowns/equations instead of n:
715 A = A[:, 0:-1]
716 A[1, 0] = 2 * (dx[-1] + dx[0])
717 A[0, 1] = dx[-1]
719 b = b[:-1]
721 # Also, due to the periodicity, the system is not tri-diagonal.
722 # We need to compute a "condensed" matrix of shape (n-2, n-2).
723 # See https://web.archive.org/web/20151220180652/http://www.cfm.brown.edu/people/gk/chap6/node14.html
724 # for more explanations.
725 # The condensed matrix is obtained by removing the last column
726 # and last row of the (n-1, n-1) system matrix. The removed
727 # values are saved in scalar variables with the (n-1, n-1)
728 # system matrix indices forming their names:
729 a_m1_0 = dx[-2] # lower left corner value: A[-1, 0]
730 a_m1_m2 = dx[-1]
731 a_m1_m1 = 2 * (dx[-1] + dx[-2])
732 a_m2_m1 = dx[-3]
733 a_0_m1 = dx[0]
735 b[0] = 3 * (dxr[0] * slope[-1] + dxr[-1] * slope[0])
736 b[-1] = 3 * (dxr[-1] * slope[-2] + dxr[-2] * slope[-1])
738 Ac = A[:, :-1]
739 b1 = b[:-1]
740 b2 = np.zeros_like(b1)
741 b2[0] = -a_0_m1
742 b2[-1] = -a_m2_m1
744 # s1 and s2 are the solutions of (n-2, n-2) system
745 s1 = solve_banded((1, 1), Ac, b1, overwrite_ab=False,
746 overwrite_b=False, check_finite=False)
748 s2 = solve_banded((1, 1), Ac, b2, overwrite_ab=False,
749 overwrite_b=False, check_finite=False)
751 # computing the s[n-2] solution:
752 s_m1 = ((b[-1] - a_m1_0 * s1[0] - a_m1_m2 * s1[-1]) /
753 (a_m1_m1 + a_m1_0 * s2[0] + a_m1_m2 * s2[-1]))
755 # s is the solution of the (n, n) system:
756 s = np.empty((n,) + y.shape[1:], dtype=y.dtype)
757 s[:-2] = s1 + s_m1 * s2
758 s[-2] = s_m1
759 s[-1] = s[0]
760 else:
761 if bc_start == 'not-a-knot':
762 A[1, 0] = dx[1]
763 A[0, 1] = x[2] - x[0]
764 d = x[2] - x[0]
765 b[0] = ((dxr[0] + 2*d) * dxr[1] * slope[0] +
766 dxr[0]**2 * slope[1]) / d
767 elif bc_start[0] == 1:
768 A[1, 0] = 1
769 A[0, 1] = 0
770 b[0] = bc_start[1]
771 elif bc_start[0] == 2:
772 A[1, 0] = 2 * dx[0]
773 A[0, 1] = dx[0]
774 b[0] = -0.5 * bc_start[1] * dx[0]**2 + 3 * (y[1] - y[0])
776 if bc_end == 'not-a-knot':
777 A[1, -1] = dx[-2]
778 A[-1, -2] = x[-1] - x[-3]
779 d = x[-1] - x[-3]
780 b[-1] = ((dxr[-1]**2*slope[-2] +
781 (2*d + dxr[-1])*dxr[-2]*slope[-1]) / d)
782 elif bc_end[0] == 1:
783 A[1, -1] = 1
784 A[-1, -2] = 0
785 b[-1] = bc_end[1]
786 elif bc_end[0] == 2:
787 A[1, -1] = 2 * dx[-1]
788 A[-1, -2] = dx[-1]
789 b[-1] = 0.5 * bc_end[1] * dx[-1]**2 + 3 * (y[-1] - y[-2])
791 s = solve_banded((1, 1), A, b, overwrite_ab=True,
792 overwrite_b=True, check_finite=False)
794 super().__init__(x, y, s, axis=0, extrapolate=extrapolate)
795 self.axis = axis
797 @staticmethod
798 def _validate_bc(bc_type, y, expected_deriv_shape, axis):
799 """Validate and prepare boundary conditions.
801 Returns
802 -------
803 validated_bc : 2-tuple
804 Boundary conditions for a curve start and end.
805 y : ndarray
806 y casted to complex dtype if one of the boundary conditions has
807 complex dtype.
808 """
809 if isinstance(bc_type, str):
810 if bc_type == 'periodic':
811 if not np.allclose(y[0], y[-1], rtol=1e-15, atol=1e-15):
812 raise ValueError(
813 "The first and last `y` point along axis {} must "
814 "be identical (within machine precision) when "
815 "bc_type='periodic'.".format(axis))
817 bc_type = (bc_type, bc_type)
819 else:
820 if len(bc_type) != 2:
821 raise ValueError("`bc_type` must contain 2 elements to "
822 "specify start and end conditions.")
824 if 'periodic' in bc_type:
825 raise ValueError("'periodic' `bc_type` is defined for both "
826 "curve ends and cannot be used with other "
827 "boundary conditions.")
829 validated_bc = []
830 for bc in bc_type:
831 if isinstance(bc, str):
832 if bc == 'clamped':
833 validated_bc.append((1, np.zeros(expected_deriv_shape)))
834 elif bc == 'natural':
835 validated_bc.append((2, np.zeros(expected_deriv_shape)))
836 elif bc in ['not-a-knot', 'periodic']:
837 validated_bc.append(bc)
838 else:
839 raise ValueError("bc_type={} is not allowed.".format(bc))
840 else:
841 try:
842 deriv_order, deriv_value = bc
843 except Exception as e:
844 raise ValueError(
845 "A specified derivative value must be "
846 "given in the form (order, value)."
847 ) from e
849 if deriv_order not in [1, 2]:
850 raise ValueError("The specified derivative order must "
851 "be 1 or 2.")
853 deriv_value = np.asarray(deriv_value)
854 if deriv_value.shape != expected_deriv_shape:
855 raise ValueError(
856 "`deriv_value` shape {} is not the expected one {}."
857 .format(deriv_value.shape, expected_deriv_shape))
859 if np.issubdtype(deriv_value.dtype, np.complexfloating):
860 y = y.astype(complex, copy=False)
862 validated_bc.append((deriv_order, deriv_value))
864 return validated_bc, y