Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_cubic.py: 10%
262 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
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`={} 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, shape (npoints, )
168 A 1-D array of monotonically increasing real values. ``x`` cannot
169 include duplicate values (otherwise f is overspecified)
170 y : ndarray, shape (..., npoints, ...)
171 A N-D array of real values. ``y``'s length along the interpolation
172 axis must be equal to the length of ``x``. Use the ``axis``
173 parameter to select the interpolation axis.
174 axis : int, optional
175 Axis in the ``y`` array corresponding to the x-coordinate values. Defaults
176 to ``axis=0``.
177 extrapolate : bool, optional
178 Whether to extrapolate to out-of-bounds points based on first
179 and last intervals, or to return NaNs.
181 Methods
182 -------
183 __call__
184 derivative
185 antiderivative
186 roots
188 See Also
189 --------
190 CubicHermiteSpline : Piecewise-cubic interpolator.
191 Akima1DInterpolator : Akima 1D interpolator.
192 CubicSpline : Cubic spline data interpolator.
193 PPoly : Piecewise polynomial in terms of coefficients and breakpoints.
195 Notes
196 -----
197 The interpolator preserves monotonicity in the interpolation data and does
198 not overshoot if the data is not smooth.
200 The first derivatives are guaranteed to be continuous, but the second
201 derivatives may jump at :math:`x_k`.
203 Determines the derivatives at the points :math:`x_k`, :math:`f'_k`,
204 by using PCHIP algorithm [1]_.
206 Let :math:`h_k = x_{k+1} - x_k`, and :math:`d_k = (y_{k+1} - y_k) / h_k`
207 are the slopes at internal points :math:`x_k`.
208 If the signs of :math:`d_k` and :math:`d_{k-1}` are different or either of
209 them equals zero, then :math:`f'_k = 0`. Otherwise, it is given by the
210 weighted harmonic mean
212 .. math::
214 \frac{w_1 + w_2}{f'_k} = \frac{w_1}{d_{k-1}} + \frac{w_2}{d_k}
216 where :math:`w_1 = 2 h_k + h_{k-1}` and :math:`w_2 = h_k + 2 h_{k-1}`.
218 The end slopes are set using a one-sided scheme [2]_.
221 References
222 ----------
223 .. [1] F. N. Fritsch and J. Butland,
224 A method for constructing local
225 monotone piecewise cubic interpolants,
226 SIAM J. Sci. Comput., 5(2), 300-304 (1984).
227 :doi:`10.1137/0905021`.
228 .. [2] see, e.g., C. Moler, Numerical Computing with Matlab, 2004.
229 :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 Returns
332 -------
333 y : scalar or array_like
334 The result, of length R or length M or M by R.
336 See Also
337 --------
338 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
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 (npoints, )
379 1-D array of monotonically increasing real values.
380 y : ndarray, shape (..., npoints, ...)
381 N-D array of real values. The length of ``y`` along the interpolation axis
382 must be equal to the length of ``x``. Use the ``axis`` parameter to
383 select the interpolation axis.
384 axis : int, optional
385 Axis in the ``y`` array corresponding to the x-coordinate values. Defaults
386 to ``axis=0``.
388 Methods
389 -------
390 __call__
391 derivative
392 antiderivative
393 roots
395 See Also
396 --------
397 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
398 CubicSpline : Cubic spline data interpolator.
399 PPoly : Piecewise polynomial in terms of coefficients and breakpoints
401 Notes
402 -----
403 .. versionadded:: 0.14
405 Use only for precise data, as the fitted curve passes through the given
406 points exactly. This routine is useful for plotting a pleasingly smooth
407 curve through a few given points for purposes of plotting.
409 References
410 ----------
411 [1] A new method of interpolation and smooth curve fitting based
412 on local procedures. Hiroshi Akima, J. ACM, October 1970, 17(4),
413 589-602.
415 """
417 def __init__(self, x, y, axis=0):
418 # Original implementation in MATLAB by N. Shamsundar (BSD licensed), see
419 # https://www.mathworks.com/matlabcentral/fileexchange/1814-akima-interpolation
420 x, dx, y, axis, _ = prepare_input(x, y, axis)
422 # determine slopes between breakpoints
423 m = np.empty((x.size + 3, ) + y.shape[1:])
424 dx = dx[(slice(None), ) + (None, ) * (y.ndim - 1)]
425 m[2:-2] = np.diff(y, axis=0) / dx
427 # add two additional points on the left ...
428 m[1] = 2. * m[2] - m[3]
429 m[0] = 2. * m[1] - m[2]
430 # ... and on the right
431 m[-2] = 2. * m[-3] - m[-4]
432 m[-1] = 2. * m[-2] - m[-3]
434 # if m1 == m2 != m3 == m4, the slope at the breakpoint is not
435 # defined. This is the fill value:
436 t = .5 * (m[3:] + m[:-3])
437 # get the denominator of the slope t
438 dm = np.abs(np.diff(m, axis=0))
439 f1 = dm[2:]
440 f2 = dm[:-2]
441 f12 = f1 + f2
442 # These are the mask of where the slope at breakpoint is defined:
443 ind = np.nonzero(f12 > 1e-9 * np.max(f12, initial=-np.inf))
444 x_ind, y_ind = ind[0], ind[1:]
445 # Set the slope at breakpoint
446 t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
447 f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
449 super().__init__(x, y, t, axis=0, extrapolate=False)
450 self.axis = axis
452 def extend(self, c, x, right=True):
453 raise NotImplementedError("Extending a 1-D Akima interpolator is not "
454 "yet implemented")
456 # These are inherited from PPoly, but they do not produce an Akima
457 # interpolator. Hence stub them out.
458 @classmethod
459 def from_spline(cls, tck, extrapolate=None):
460 raise NotImplementedError("This method does not make sense for "
461 "an Akima interpolator.")
463 @classmethod
464 def from_bernstein_basis(cls, bp, extrapolate=None):
465 raise NotImplementedError("This method does not make sense for "
466 "an Akima interpolator.")
469class CubicSpline(CubicHermiteSpline):
470 """Cubic spline data interpolator.
472 Interpolate data with a piecewise cubic polynomial which is twice
473 continuously differentiable [1]_. The result is represented as a `PPoly`
474 instance with breakpoints matching the given data.
476 Parameters
477 ----------
478 x : array_like, shape (n,)
479 1-D array containing values of the independent variable.
480 Values must be real, finite and in strictly increasing order.
481 y : array_like
482 Array containing values of the dependent variable. It can have
483 arbitrary number of dimensions, but the length along ``axis``
484 (see below) must match the length of ``x``. Values must be finite.
485 axis : int, optional
486 Axis along which `y` is assumed to be varying. Meaning that for
487 ``x[i]`` the corresponding values are ``np.take(y, i, axis=axis)``.
488 Default is 0.
489 bc_type : string or 2-tuple, optional
490 Boundary condition type. Two additional equations, given by the
491 boundary conditions, are required to determine all coefficients of
492 polynomials on each segment [2]_.
494 If `bc_type` is a string, then the specified condition will be applied
495 at both ends of a spline. Available conditions are:
497 * 'not-a-knot' (default): The first and second segment at a curve end
498 are the same polynomial. It is a good default when there is no
499 information on boundary conditions.
500 * 'periodic': The interpolated functions is assumed to be periodic
501 of period ``x[-1] - x[0]``. The first and last value of `y` must be
502 identical: ``y[0] == y[-1]``. This boundary condition will result in
503 ``y'[0] == y'[-1]`` and ``y''[0] == y''[-1]``.
504 * 'clamped': The first derivative at curves ends are zero. Assuming
505 a 1D `y`, ``bc_type=((1, 0.0), (1, 0.0))`` is the same condition.
506 * 'natural': The second derivative at curve ends are zero. Assuming
507 a 1D `y`, ``bc_type=((2, 0.0), (2, 0.0))`` is the same condition.
509 If `bc_type` is a 2-tuple, the first and the second value will be
510 applied at the curve start and end respectively. The tuple values can
511 be one of the previously mentioned strings (except 'periodic') or a
512 tuple `(order, deriv_values)` allowing to specify arbitrary
513 derivatives at curve ends:
515 * `order`: the derivative order, 1 or 2.
516 * `deriv_value`: array_like containing derivative values, shape must
517 be the same as `y`, excluding ``axis`` dimension. For example, if
518 `y` is 1-D, then `deriv_value` must be a scalar. If `y` is 3-D with
519 the shape (n0, n1, n2) and axis=2, then `deriv_value` must be 2-D
520 and have the shape (n0, n1).
521 extrapolate : {bool, 'periodic', None}, optional
522 If bool, determines whether to extrapolate to out-of-bounds points
523 based on first and last intervals, or to return NaNs. If 'periodic',
524 periodic extrapolation is used. If None (default), ``extrapolate`` is
525 set to 'periodic' for ``bc_type='periodic'`` and to True otherwise.
527 Attributes
528 ----------
529 x : ndarray, shape (n,)
530 Breakpoints. The same ``x`` which was passed to the constructor.
531 c : ndarray, shape (4, n-1, ...)
532 Coefficients of the polynomials on each segment. The trailing
533 dimensions match the dimensions of `y`, excluding ``axis``.
534 For example, if `y` is 1-d, then ``c[k, i]`` is a coefficient for
535 ``(x-x[i])**(3-k)`` on the segment between ``x[i]`` and ``x[i+1]``.
536 axis : int
537 Interpolation axis. The same axis which was passed to the
538 constructor.
540 Methods
541 -------
542 __call__
543 derivative
544 antiderivative
545 integrate
546 roots
548 See Also
549 --------
550 Akima1DInterpolator : Akima 1D interpolator.
551 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator.
552 PPoly : Piecewise polynomial in terms of coefficients and breakpoints.
554 Notes
555 -----
556 Parameters `bc_type` and ``extrapolate`` work independently, i.e. the
557 former controls only construction of a spline, and the latter only
558 evaluation.
560 When a boundary condition is 'not-a-knot' and n = 2, it is replaced by
561 a condition that the first derivative is equal to the linear interpolant
562 slope. When both boundary conditions are 'not-a-knot' and n = 3, the
563 solution is sought as a parabola passing through given points.
565 When 'not-a-knot' boundary conditions is applied to both ends, the
566 resulting spline will be the same as returned by `splrep` (with ``s=0``)
567 and `InterpolatedUnivariateSpline`, but these two methods use a
568 representation in B-spline basis.
570 .. versionadded:: 0.18.0
572 Examples
573 --------
574 In this example the cubic spline is used to interpolate a sampled sinusoid.
575 You can see that the spline continuity property holds for the first and
576 second derivatives and violates only for the third derivative.
578 >>> import numpy as np
579 >>> from scipy.interpolate import CubicSpline
580 >>> import matplotlib.pyplot as plt
581 >>> x = np.arange(10)
582 >>> y = np.sin(x)
583 >>> cs = CubicSpline(x, y)
584 >>> xs = np.arange(-0.5, 9.6, 0.1)
585 >>> fig, ax = plt.subplots(figsize=(6.5, 4))
586 >>> ax.plot(x, y, 'o', label='data')
587 >>> ax.plot(xs, np.sin(xs), label='true')
588 >>> ax.plot(xs, cs(xs), label="S")
589 >>> ax.plot(xs, cs(xs, 1), label="S'")
590 >>> ax.plot(xs, cs(xs, 2), label="S''")
591 >>> ax.plot(xs, cs(xs, 3), label="S'''")
592 >>> ax.set_xlim(-0.5, 9.5)
593 >>> ax.legend(loc='lower left', ncol=2)
594 >>> plt.show()
596 In the second example, the unit circle is interpolated with a spline. A
597 periodic boundary condition is used. You can see that the first derivative
598 values, ds/dx=0, ds/dy=1 at the periodic point (1, 0) are correctly
599 computed. Note that a circle cannot be exactly represented by a cubic
600 spline. To increase precision, more breakpoints would be required.
602 >>> theta = 2 * np.pi * np.linspace(0, 1, 5)
603 >>> y = np.c_[np.cos(theta), np.sin(theta)]
604 >>> cs = CubicSpline(theta, y, bc_type='periodic')
605 >>> print("ds/dx={:.1f} ds/dy={:.1f}".format(cs(0, 1)[0], cs(0, 1)[1]))
606 ds/dx=0.0 ds/dy=1.0
607 >>> xs = 2 * np.pi * np.linspace(0, 1, 100)
608 >>> fig, ax = plt.subplots(figsize=(6.5, 4))
609 >>> ax.plot(y[:, 0], y[:, 1], 'o', label='data')
610 >>> ax.plot(np.cos(xs), np.sin(xs), label='true')
611 >>> ax.plot(cs(xs)[:, 0], cs(xs)[:, 1], label='spline')
612 >>> ax.axes.set_aspect('equal')
613 >>> ax.legend(loc='center')
614 >>> plt.show()
616 The third example is the interpolation of a polynomial y = x**3 on the
617 interval 0 <= x<= 1. A cubic spline can represent this function exactly.
618 To achieve that we need to specify values and first derivatives at
619 endpoints of the interval. Note that y' = 3 * x**2 and thus y'(0) = 0 and
620 y'(1) = 3.
622 >>> cs = CubicSpline([0, 1], [0, 1], bc_type=((1, 0), (1, 3)))
623 >>> x = np.linspace(0, 1)
624 >>> np.allclose(x**3, cs(x))
625 True
627 References
628 ----------
629 .. [1] `Cubic Spline Interpolation
630 <https://en.wikiversity.org/wiki/Cubic_Spline_Interpolation>`_
631 on Wikiversity.
632 .. [2] Carl de Boor, "A Practical Guide to Splines", Springer-Verlag, 1978.
633 """
635 def __init__(self, x, y, axis=0, bc_type='not-a-knot', extrapolate=None):
636 x, dx, y, axis, _ = prepare_input(x, y, axis)
637 n = len(x)
639 bc, y = self._validate_bc(bc_type, y, y.shape[1:], axis)
641 if extrapolate is None:
642 if bc[0] == 'periodic':
643 extrapolate = 'periodic'
644 else:
645 extrapolate = True
647 if y.size == 0:
648 # bail out early for zero-sized arrays
649 s = np.zeros_like(y)
650 else:
651 dxr = dx.reshape([dx.shape[0]] + [1] * (y.ndim - 1))
652 slope = np.diff(y, axis=0) / dxr
654 # If bc is 'not-a-knot' this change is just a convention.
655 # If bc is 'periodic' then we already checked that y[0] == y[-1],
656 # and the spline is just a constant, we handle this case in the
657 # same way by setting the first derivatives to slope, which is 0.
658 if n == 2:
659 if bc[0] in ['not-a-knot', 'periodic']:
660 bc[0] = (1, slope[0])
661 if bc[1] in ['not-a-knot', 'periodic']:
662 bc[1] = (1, slope[0])
664 # This is a special case, when both conditions are 'not-a-knot'
665 # and n == 3. In this case 'not-a-knot' can't be handled regularly
666 # as the both conditions are identical. We handle this case by
667 # constructing a parabola passing through given points.
668 if n == 3 and bc[0] == 'not-a-knot' and bc[1] == 'not-a-knot':
669 A = np.zeros((3, 3)) # This is a standard matrix.
670 b = np.empty((3,) + y.shape[1:], dtype=y.dtype)
672 A[0, 0] = 1
673 A[0, 1] = 1
674 A[1, 0] = dx[1]
675 A[1, 1] = 2 * (dx[0] + dx[1])
676 A[1, 2] = dx[0]
677 A[2, 1] = 1
678 A[2, 2] = 1
680 b[0] = 2 * slope[0]
681 b[1] = 3 * (dxr[0] * slope[1] + dxr[1] * slope[0])
682 b[2] = 2 * slope[1]
684 s = solve(A, b, overwrite_a=True, overwrite_b=True,
685 check_finite=False)
686 elif n == 3 and bc[0] == 'periodic':
687 # In case when number of points is 3 we compute the derivatives
688 # manually
689 s = np.empty((n,) + y.shape[1:], dtype=y.dtype)
690 t = (slope / dxr).sum() / (1. / dxr).sum()
691 s.fill(t)
692 else:
693 # Find derivative values at each x[i] by solving a tridiagonal
694 # system.
695 A = np.zeros((3, n)) # This is a banded matrix representation.
696 b = np.empty((n,) + y.shape[1:], dtype=y.dtype)
698 # Filling the system for i=1..n-2
699 # (x[i-1] - x[i]) * s[i-1] +\
700 # 2 * ((x[i] - x[i-1]) + (x[i+1] - x[i])) * s[i] +\
701 # (x[i] - x[i-1]) * s[i+1] =\
702 # 3 * ((x[i+1] - x[i])*(y[i] - y[i-1])/(x[i] - x[i-1]) +\
703 # (x[i] - x[i-1])*(y[i+1] - y[i])/(x[i+1] - x[i]))
705 A[1, 1:-1] = 2 * (dx[:-1] + dx[1:]) # The diagonal
706 A[0, 2:] = dx[:-1] # The upper diagonal
707 A[-1, :-2] = dx[1:] # The lower diagonal
709 b[1:-1] = 3 * (dxr[1:] * slope[:-1] + dxr[:-1] * slope[1:])
711 bc_start, bc_end = bc
713 if bc_start == 'periodic':
714 # Due to the periodicity, and because y[-1] = y[0], the
715 # linear system has (n-1) unknowns/equations instead of n:
716 A = A[:, 0:-1]
717 A[1, 0] = 2 * (dx[-1] + dx[0])
718 A[0, 1] = dx[-1]
720 b = b[:-1]
722 # Also, due to the periodicity, the system is not tri-diagonal.
723 # We need to compute a "condensed" matrix of shape (n-2, n-2).
724 # See https://web.archive.org/web/20151220180652/http://www.cfm.brown.edu/people/gk/chap6/node14.html
725 # for more explanations.
726 # The condensed matrix is obtained by removing the last column
727 # and last row of the (n-1, n-1) system matrix. The removed
728 # values are saved in scalar variables with the (n-1, n-1)
729 # system matrix indices forming their names:
730 a_m1_0 = dx[-2] # lower left corner value: A[-1, 0]
731 a_m1_m2 = dx[-1]
732 a_m1_m1 = 2 * (dx[-1] + dx[-2])
733 a_m2_m1 = dx[-3]
734 a_0_m1 = dx[0]
736 b[0] = 3 * (dxr[0] * slope[-1] + dxr[-1] * slope[0])
737 b[-1] = 3 * (dxr[-1] * slope[-2] + dxr[-2] * slope[-1])
739 Ac = A[:, :-1]
740 b1 = b[:-1]
741 b2 = np.zeros_like(b1)
742 b2[0] = -a_0_m1
743 b2[-1] = -a_m2_m1
745 # s1 and s2 are the solutions of (n-2, n-2) system
746 s1 = solve_banded((1, 1), Ac, b1, overwrite_ab=False,
747 overwrite_b=False, check_finite=False)
749 s2 = solve_banded((1, 1), Ac, b2, overwrite_ab=False,
750 overwrite_b=False, check_finite=False)
752 # computing the s[n-2] solution:
753 s_m1 = ((b[-1] - a_m1_0 * s1[0] - a_m1_m2 * s1[-1]) /
754 (a_m1_m1 + a_m1_0 * s2[0] + a_m1_m2 * s2[-1]))
756 # s is the solution of the (n, n) system:
757 s = np.empty((n,) + y.shape[1:], dtype=y.dtype)
758 s[:-2] = s1 + s_m1 * s2
759 s[-2] = s_m1
760 s[-1] = s[0]
761 else:
762 if bc_start == 'not-a-knot':
763 A[1, 0] = dx[1]
764 A[0, 1] = x[2] - x[0]
765 d = x[2] - x[0]
766 b[0] = ((dxr[0] + 2*d) * dxr[1] * slope[0] +
767 dxr[0]**2 * slope[1]) / d
768 elif bc_start[0] == 1:
769 A[1, 0] = 1
770 A[0, 1] = 0
771 b[0] = bc_start[1]
772 elif bc_start[0] == 2:
773 A[1, 0] = 2 * dx[0]
774 A[0, 1] = dx[0]
775 b[0] = -0.5 * bc_start[1] * dx[0]**2 + 3 * (y[1] - y[0])
777 if bc_end == 'not-a-knot':
778 A[1, -1] = dx[-2]
779 A[-1, -2] = x[-1] - x[-3]
780 d = x[-1] - x[-3]
781 b[-1] = ((dxr[-1]**2*slope[-2] +
782 (2*d + dxr[-1])*dxr[-2]*slope[-1]) / d)
783 elif bc_end[0] == 1:
784 A[1, -1] = 1
785 A[-1, -2] = 0
786 b[-1] = bc_end[1]
787 elif bc_end[0] == 2:
788 A[1, -1] = 2 * dx[-1]
789 A[-1, -2] = dx[-1]
790 b[-1] = 0.5 * bc_end[1] * dx[-1]**2 + 3 * (y[-1] - y[-2])
792 s = solve_banded((1, 1), A, b, overwrite_ab=True,
793 overwrite_b=True, check_finite=False)
795 super().__init__(x, y, s, axis=0, extrapolate=extrapolate)
796 self.axis = axis
798 @staticmethod
799 def _validate_bc(bc_type, y, expected_deriv_shape, axis):
800 """Validate and prepare boundary conditions.
802 Returns
803 -------
804 validated_bc : 2-tuple
805 Boundary conditions for a curve start and end.
806 y : ndarray
807 y casted to complex dtype if one of the boundary conditions has
808 complex dtype.
809 """
810 if isinstance(bc_type, str):
811 if bc_type == 'periodic':
812 if not np.allclose(y[0], y[-1], rtol=1e-15, atol=1e-15):
813 raise ValueError(
814 "The first and last `y` point along axis {} must "
815 "be identical (within machine precision) when "
816 "bc_type='periodic'.".format(axis))
818 bc_type = (bc_type, bc_type)
820 else:
821 if len(bc_type) != 2:
822 raise ValueError("`bc_type` must contain 2 elements to "
823 "specify start and end conditions.")
825 if 'periodic' in bc_type:
826 raise ValueError("'periodic' `bc_type` is defined for both "
827 "curve ends and cannot be used with other "
828 "boundary conditions.")
830 validated_bc = []
831 for bc in bc_type:
832 if isinstance(bc, str):
833 if bc == 'clamped':
834 validated_bc.append((1, np.zeros(expected_deriv_shape)))
835 elif bc == 'natural':
836 validated_bc.append((2, np.zeros(expected_deriv_shape)))
837 elif bc in ['not-a-knot', 'periodic']:
838 validated_bc.append(bc)
839 else:
840 raise ValueError(f"bc_type={bc} is not allowed.")
841 else:
842 try:
843 deriv_order, deriv_value = bc
844 except Exception as e:
845 raise ValueError(
846 "A specified derivative value must be "
847 "given in the form (order, value)."
848 ) from e
850 if deriv_order not in [1, 2]:
851 raise ValueError("The specified derivative order must "
852 "be 1 or 2.")
854 deriv_value = np.asarray(deriv_value)
855 if deriv_value.shape != expected_deriv_shape:
856 raise ValueError(
857 "`deriv_value` shape {} is not the expected one {}."
858 .format(deriv_value.shape, expected_deriv_shape))
860 if np.issubdtype(deriv_value.dtype, np.complexfloating):
861 y = y.astype(complex, copy=False)
863 validated_bc.append((deriv_order, deriv_value))
865 return validated_bc, y