Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_bsplines.py: 8%
579 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import operator
2from math import prod
4import numpy as np
5from numpy.core.multiarray import normalize_axis_index
6from scipy.linalg import (get_lapack_funcs, LinAlgError,
7 cholesky_banded, cho_solve_banded,
8 solve, solve_banded)
9from scipy.optimize import minimize_scalar
10from . import _bspl
11from . import _fitpack_impl
12from scipy.sparse import csr_array
13from scipy.special import poch
14from itertools import combinations
16__all__ = ["BSpline", "make_interp_spline", "make_lsq_spline",
17 "make_smoothing_spline"]
20def _get_dtype(dtype):
21 """Return np.complex128 for complex dtypes, np.float64 otherwise."""
22 if np.issubdtype(dtype, np.complexfloating):
23 return np.complex_
24 else:
25 return np.float_
28def _as_float_array(x, check_finite=False):
29 """Convert the input into a C contiguous float array.
31 NB: Upcasts half- and single-precision floats to double precision.
32 """
33 x = np.ascontiguousarray(x)
34 dtyp = _get_dtype(x.dtype)
35 x = x.astype(dtyp, copy=False)
36 if check_finite and not np.isfinite(x).all():
37 raise ValueError("Array must not contain infs or nans.")
38 return x
41def _dual_poly(j, k, t, y):
42 """
43 Dual polynomial of the B-spline B_{j,k,t} -
44 polynomial which is associated with B_{j,k,t}:
45 $p_{j,k}(y) = (y - t_{j+1})(y - t_{j+2})...(y - t_{j+k})$
46 """
47 if k == 0:
48 return 1
49 return np.prod([(y - t[j + i]) for i in range(1, k + 1)])
52def _diff_dual_poly(j, k, y, d, t):
53 """
54 d-th derivative of the dual polynomial $p_{j,k}(y)$
55 """
56 if d == 0:
57 return _dual_poly(j, k, t, y)
58 if d == k:
59 return poch(1, k)
60 comb = list(combinations(range(j + 1, j + k + 1), d))
61 res = 0
62 for i in range(len(comb) * len(comb[0])):
63 res += np.prod([(y - t[j + p]) for p in range(1, k + 1)
64 if (j + p) not in comb[i//d]])
65 return res
68class BSpline:
69 r"""Univariate spline in the B-spline basis.
71 .. math::
73 S(x) = \sum_{j=0}^{n-1} c_j B_{j, k; t}(x)
75 where :math:`B_{j, k; t}` are B-spline basis functions of degree `k`
76 and knots `t`.
78 Parameters
79 ----------
80 t : ndarray, shape (n+k+1,)
81 knots
82 c : ndarray, shape (>=n, ...)
83 spline coefficients
84 k : int
85 B-spline degree
86 extrapolate : bool or 'periodic', optional
87 whether to extrapolate beyond the base interval, ``t[k] .. t[n]``,
88 or to return nans.
89 If True, extrapolates the first and last polynomial pieces of b-spline
90 functions active on the base interval.
91 If 'periodic', periodic extrapolation is used.
92 Default is True.
93 axis : int, optional
94 Interpolation axis. Default is zero.
96 Attributes
97 ----------
98 t : ndarray
99 knot vector
100 c : ndarray
101 spline coefficients
102 k : int
103 spline degree
104 extrapolate : bool
105 If True, extrapolates the first and last polynomial pieces of b-spline
106 functions active on the base interval.
107 axis : int
108 Interpolation axis.
109 tck : tuple
110 A read-only equivalent of ``(self.t, self.c, self.k)``
112 Methods
113 -------
114 __call__
115 basis_element
116 derivative
117 antiderivative
118 integrate
119 construct_fast
120 design_matrix
121 from_power_basis
123 Notes
124 -----
125 B-spline basis elements are defined via
127 .. math::
129 B_{i, 0}(x) = 1, \textrm{if $t_i \le x < t_{i+1}$, otherwise $0$,}
131 B_{i, k}(x) = \frac{x - t_i}{t_{i+k} - t_i} B_{i, k-1}(x)
132 + \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} B_{i+1, k-1}(x)
134 **Implementation details**
136 - At least ``k+1`` coefficients are required for a spline of degree `k`,
137 so that ``n >= k+1``. Additional coefficients, ``c[j]`` with
138 ``j > n``, are ignored.
140 - B-spline basis elements of degree `k` form a partition of unity on the
141 *base interval*, ``t[k] <= x <= t[n]``.
144 Examples
145 --------
147 Translating the recursive definition of B-splines into Python code, we have:
149 >>> def B(x, k, i, t):
150 ... if k == 0:
151 ... return 1.0 if t[i] <= x < t[i+1] else 0.0
152 ... if t[i+k] == t[i]:
153 ... c1 = 0.0
154 ... else:
155 ... c1 = (x - t[i])/(t[i+k] - t[i]) * B(x, k-1, i, t)
156 ... if t[i+k+1] == t[i+1]:
157 ... c2 = 0.0
158 ... else:
159 ... c2 = (t[i+k+1] - x)/(t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t)
160 ... return c1 + c2
162 >>> def bspline(x, t, c, k):
163 ... n = len(t) - k - 1
164 ... assert (n >= k+1) and (len(c) >= n)
165 ... return sum(c[i] * B(x, k, i, t) for i in range(n))
167 Note that this is an inefficient (if straightforward) way to
168 evaluate B-splines --- this spline class does it in an equivalent,
169 but much more efficient way.
171 Here we construct a quadratic spline function on the base interval
172 ``2 <= x <= 4`` and compare with the naive way of evaluating the spline:
174 >>> from scipy.interpolate import BSpline
175 >>> k = 2
176 >>> t = [0, 1, 2, 3, 4, 5, 6]
177 >>> c = [-1, 2, 0, -1]
178 >>> spl = BSpline(t, c, k)
179 >>> spl(2.5)
180 array(1.375)
181 >>> bspline(2.5, t, c, k)
182 1.375
184 Note that outside of the base interval results differ. This is because
185 `BSpline` extrapolates the first and last polynomial pieces of B-spline
186 functions active on the base interval.
188 >>> import matplotlib.pyplot as plt
189 >>> import numpy as np
190 >>> fig, ax = plt.subplots()
191 >>> xx = np.linspace(1.5, 4.5, 50)
192 >>> ax.plot(xx, [bspline(x, t, c ,k) for x in xx], 'r-', lw=3, label='naive')
193 >>> ax.plot(xx, spl(xx), 'b-', lw=4, alpha=0.7, label='BSpline')
194 >>> ax.grid(True)
195 >>> ax.legend(loc='best')
196 >>> plt.show()
199 References
200 ----------
201 .. [1] Tom Lyche and Knut Morken, Spline methods,
202 http://www.uio.no/studier/emner/matnat/ifi/INF-MAT5340/v05/undervisningsmateriale/
203 .. [2] Carl de Boor, A practical guide to splines, Springer, 2001.
205 """
207 def __init__(self, t, c, k, extrapolate=True, axis=0):
208 super().__init__()
210 self.k = operator.index(k)
211 self.c = np.asarray(c)
212 self.t = np.ascontiguousarray(t, dtype=np.float64)
214 if extrapolate == 'periodic':
215 self.extrapolate = extrapolate
216 else:
217 self.extrapolate = bool(extrapolate)
219 n = self.t.shape[0] - self.k - 1
221 axis = normalize_axis_index(axis, self.c.ndim)
223 # Note that the normalized axis is stored in the object.
224 self.axis = axis
225 if axis != 0:
226 # roll the interpolation axis to be the first one in self.c
227 # More specifically, the target shape for self.c is (n, ...),
228 # and axis !=0 means that we have c.shape (..., n, ...)
229 # ^
230 # axis
231 self.c = np.moveaxis(self.c, axis, 0)
233 if k < 0:
234 raise ValueError("Spline order cannot be negative.")
235 if self.t.ndim != 1:
236 raise ValueError("Knot vector must be one-dimensional.")
237 if n < self.k + 1:
238 raise ValueError("Need at least %d knots for degree %d" %
239 (2*k + 2, k))
240 if (np.diff(self.t) < 0).any():
241 raise ValueError("Knots must be in a non-decreasing order.")
242 if len(np.unique(self.t[k:n+1])) < 2:
243 raise ValueError("Need at least two internal knots.")
244 if not np.isfinite(self.t).all():
245 raise ValueError("Knots should not have nans or infs.")
246 if self.c.ndim < 1:
247 raise ValueError("Coefficients must be at least 1-dimensional.")
248 if self.c.shape[0] < n:
249 raise ValueError("Knots, coefficients and degree are inconsistent.")
251 dt = _get_dtype(self.c.dtype)
252 self.c = np.ascontiguousarray(self.c, dtype=dt)
254 @classmethod
255 def construct_fast(cls, t, c, k, extrapolate=True, axis=0):
256 """Construct a spline without making checks.
258 Accepts same parameters as the regular constructor. Input arrays
259 `t` and `c` must of correct shape and dtype.
260 """
261 self = object.__new__(cls)
262 self.t, self.c, self.k = t, c, k
263 self.extrapolate = extrapolate
264 self.axis = axis
265 return self
267 @property
268 def tck(self):
269 """Equivalent to ``(self.t, self.c, self.k)`` (read-only).
270 """
271 return self.t, self.c, self.k
273 @classmethod
274 def basis_element(cls, t, extrapolate=True):
275 """Return a B-spline basis element ``B(x | t[0], ..., t[k+1])``.
277 Parameters
278 ----------
279 t : ndarray, shape (k+2,)
280 internal knots
281 extrapolate : bool or 'periodic', optional
282 whether to extrapolate beyond the base interval, ``t[0] .. t[k+1]``,
283 or to return nans.
284 If 'periodic', periodic extrapolation is used.
285 Default is True.
287 Returns
288 -------
289 basis_element : callable
290 A callable representing a B-spline basis element for the knot
291 vector `t`.
293 Notes
294 -----
295 The degree of the B-spline, `k`, is inferred from the length of `t` as
296 ``len(t)-2``. The knot vector is constructed by appending and prepending
297 ``k+1`` elements to internal knots `t`.
299 Examples
300 --------
302 Construct a cubic B-spline:
304 >>> import numpy as np
305 >>> from scipy.interpolate import BSpline
306 >>> b = BSpline.basis_element([0, 1, 2, 3, 4])
307 >>> k = b.k
308 >>> b.t[k:-k]
309 array([ 0., 1., 2., 3., 4.])
310 >>> k
311 3
313 Construct a quadratic B-spline on ``[0, 1, 1, 2]``, and compare
314 to its explicit form:
316 >>> t = [0, 1, 1, 2]
317 >>> b = BSpline.basis_element(t)
318 >>> def f(x):
319 ... return np.where(x < 1, x*x, (2. - x)**2)
321 >>> import matplotlib.pyplot as plt
322 >>> fig, ax = plt.subplots()
323 >>> x = np.linspace(0, 2, 51)
324 >>> ax.plot(x, b(x), 'g', lw=3)
325 >>> ax.plot(x, f(x), 'r', lw=8, alpha=0.4)
326 >>> ax.grid(True)
327 >>> plt.show()
329 """
330 k = len(t) - 2
331 t = _as_float_array(t)
332 t = np.r_[(t[0]-1,) * k, t, (t[-1]+1,) * k]
333 c = np.zeros_like(t)
334 c[k] = 1.
335 return cls.construct_fast(t, c, k, extrapolate)
337 @classmethod
338 def design_matrix(cls, x, t, k, extrapolate=False):
339 """
340 Returns a design matrix as a CSR format sparse array.
342 Parameters
343 ----------
344 x : array_like, shape (n,)
345 Points to evaluate the spline at.
346 t : array_like, shape (nt,)
347 Sorted 1D array of knots.
348 k : int
349 B-spline degree.
350 extrapolate : bool or 'periodic', optional
351 Whether to extrapolate based on the first and last intervals
352 or raise an error. If 'periodic', periodic extrapolation is used.
353 Default is False.
355 .. versionadded:: 1.10.0
357 Returns
358 -------
359 design_matrix : `csr_array` object
360 Sparse matrix in CSR format where each row contains all the basis
361 elements of the input row (first row = basis elements of x[0],
362 ..., last row = basis elements x[-1]).
364 Examples
365 --------
366 Construct a design matrix for a B-spline
368 >>> from scipy.interpolate import make_interp_spline, BSpline
369 >>> import numpy as np
370 >>> x = np.linspace(0, np.pi * 2, 4)
371 >>> y = np.sin(x)
372 >>> k = 3
373 >>> bspl = make_interp_spline(x, y, k=k)
374 >>> design_matrix = bspl.design_matrix(x, bspl.t, k)
375 >>> design_matrix.toarray()
376 [[1. , 0. , 0. , 0. ],
377 [0.2962963 , 0.44444444, 0.22222222, 0.03703704],
378 [0.03703704, 0.22222222, 0.44444444, 0.2962963 ],
379 [0. , 0. , 0. , 1. ]]
381 Construct a design matrix for some vector of knots
383 >>> k = 2
384 >>> t = [-1, 0, 1, 2, 3, 4, 5, 6]
385 >>> x = [1, 2, 3, 4]
386 >>> design_matrix = BSpline.design_matrix(x, t, k).toarray()
387 >>> design_matrix
388 [[0.5, 0.5, 0. , 0. , 0. ],
389 [0. , 0.5, 0.5, 0. , 0. ],
390 [0. , 0. , 0.5, 0.5, 0. ],
391 [0. , 0. , 0. , 0.5, 0.5]]
393 This result is equivalent to the one created in the sparse format
395 >>> c = np.eye(len(t) - k - 1)
396 >>> design_matrix_gh = BSpline(t, c, k)(x)
397 >>> np.allclose(design_matrix, design_matrix_gh, atol=1e-14)
398 True
400 Notes
401 -----
402 .. versionadded:: 1.8.0
404 In each row of the design matrix all the basis elements are evaluated
405 at the certain point (first row - x[0], ..., last row - x[-1]).
407 `nt` is a length of the vector of knots: as far as there are
408 `nt - k - 1` basis elements, `nt` should be not less than `2 * k + 2`
409 to have at least `k + 1` basis element.
411 Out of bounds `x` raises a ValueError.
412 """
413 x = _as_float_array(x, True)
414 t = _as_float_array(t, True)
416 if extrapolate != 'periodic':
417 extrapolate = bool(extrapolate)
419 if k < 0:
420 raise ValueError("Spline order cannot be negative.")
421 if t.ndim != 1 or np.any(t[1:] < t[:-1]):
422 raise ValueError(f"Expect t to be a 1-D sorted array_like, but "
423 f"got t={t}.")
424 # There are `nt - k - 1` basis elements in a BSpline built on the
425 # vector of knots with length `nt`, so to have at least `k + 1` basis
426 # elements we need to have at least `2 * k + 2` elements in the vector
427 # of knots.
428 if len(t) < 2 * k + 2:
429 raise ValueError(f"Length t is not enough for k={k}.")
431 if extrapolate == 'periodic':
432 # With periodic extrapolation we map x to the segment
433 # [t[k], t[n]].
434 n = t.size - k - 1
435 x = t[k] + (x - t[k]) % (t[n] - t[k])
436 extrapolate = False
437 elif not extrapolate and (
438 (min(x) < t[k]) or (max(x) > t[t.shape[0] - k - 1])
439 ):
440 # Checks from `find_interval` function
441 raise ValueError(f'Out of bounds w/ x = {x}.')
443 # Compute number of non-zeros of final CSR array in order to determine
444 # the dtype of indices and indptr of the CSR array.
445 n = x.shape[0]
446 nnz = n * (k + 1)
447 if nnz < np.iinfo(np.int32).max:
448 int_dtype = np.int32
449 else:
450 int_dtype = np.int64
451 # Preallocate indptr and indices
452 indices = np.empty(n * (k + 1), dtype=int_dtype)
453 indptr = np.arange(0, (n + 1) * (k + 1), k + 1, dtype=int_dtype)
455 # indptr is not passed to Cython as it is already fully computed
456 data, indices = _bspl._make_design_matrix(
457 x, t, k, extrapolate, indices
458 )
459 return csr_array(
460 (data, indices, indptr),
461 shape=(x.shape[0], t.shape[0] - k - 1)
462 )
464 def __call__(self, x, nu=0, extrapolate=None):
465 """
466 Evaluate a spline function.
468 Parameters
469 ----------
470 x : array_like
471 points to evaluate the spline at.
472 nu : int, optional
473 derivative to evaluate (default is 0).
474 extrapolate : bool or 'periodic', optional
475 whether to extrapolate based on the first and last intervals
476 or return nans. If 'periodic', periodic extrapolation is used.
477 Default is `self.extrapolate`.
479 Returns
480 -------
481 y : array_like
482 Shape is determined by replacing the interpolation axis
483 in the coefficient array with the shape of `x`.
485 """
486 if extrapolate is None:
487 extrapolate = self.extrapolate
488 x = np.asarray(x)
489 x_shape, x_ndim = x.shape, x.ndim
490 x = np.ascontiguousarray(x.ravel(), dtype=np.float_)
492 # With periodic extrapolation we map x to the segment
493 # [self.t[k], self.t[n]].
494 if extrapolate == 'periodic':
495 n = self.t.size - self.k - 1
496 x = self.t[self.k] + (x - self.t[self.k]) % (self.t[n] -
497 self.t[self.k])
498 extrapolate = False
500 out = np.empty((len(x), prod(self.c.shape[1:])), dtype=self.c.dtype)
501 self._ensure_c_contiguous()
502 self._evaluate(x, nu, extrapolate, out)
503 out = out.reshape(x_shape + self.c.shape[1:])
504 if self.axis != 0:
505 # transpose to move the calculated values to the interpolation axis
506 l = list(range(out.ndim))
507 l = l[x_ndim:x_ndim+self.axis] + l[:x_ndim] + l[x_ndim+self.axis:]
508 out = out.transpose(l)
509 return out
511 def _evaluate(self, xp, nu, extrapolate, out):
512 _bspl.evaluate_spline(self.t, self.c.reshape(self.c.shape[0], -1),
513 self.k, xp, nu, extrapolate, out)
515 def _ensure_c_contiguous(self):
516 """
517 c and t may be modified by the user. The Cython code expects
518 that they are C contiguous.
520 """
521 if not self.t.flags.c_contiguous:
522 self.t = self.t.copy()
523 if not self.c.flags.c_contiguous:
524 self.c = self.c.copy()
526 def derivative(self, nu=1):
527 """Return a B-spline representing the derivative.
529 Parameters
530 ----------
531 nu : int, optional
532 Derivative order.
533 Default is 1.
535 Returns
536 -------
537 b : BSpline object
538 A new instance representing the derivative.
540 See Also
541 --------
542 splder, splantider
544 """
545 c = self.c
546 # pad the c array if needed
547 ct = len(self.t) - len(c)
548 if ct > 0:
549 c = np.r_[c, np.zeros((ct,) + c.shape[1:])]
550 tck = _fitpack_impl.splder((self.t, c, self.k), nu)
551 return self.construct_fast(*tck, extrapolate=self.extrapolate,
552 axis=self.axis)
554 def antiderivative(self, nu=1):
555 """Return a B-spline representing the antiderivative.
557 Parameters
558 ----------
559 nu : int, optional
560 Antiderivative order. Default is 1.
562 Returns
563 -------
564 b : BSpline object
565 A new instance representing the antiderivative.
567 Notes
568 -----
569 If antiderivative is computed and ``self.extrapolate='periodic'``,
570 it will be set to False for the returned instance. This is done because
571 the antiderivative is no longer periodic and its correct evaluation
572 outside of the initially given x interval is difficult.
574 See Also
575 --------
576 splder, splantider
578 """
579 c = self.c
580 # pad the c array if needed
581 ct = len(self.t) - len(c)
582 if ct > 0:
583 c = np.r_[c, np.zeros((ct,) + c.shape[1:])]
584 tck = _fitpack_impl.splantider((self.t, c, self.k), nu)
586 if self.extrapolate == 'periodic':
587 extrapolate = False
588 else:
589 extrapolate = self.extrapolate
591 return self.construct_fast(*tck, extrapolate=extrapolate,
592 axis=self.axis)
594 def integrate(self, a, b, extrapolate=None):
595 """Compute a definite integral of the spline.
597 Parameters
598 ----------
599 a : float
600 Lower limit of integration.
601 b : float
602 Upper limit of integration.
603 extrapolate : bool or 'periodic', optional
604 whether to extrapolate beyond the base interval,
605 ``t[k] .. t[-k-1]``, or take the spline to be zero outside of the
606 base interval. If 'periodic', periodic extrapolation is used.
607 If None (default), use `self.extrapolate`.
609 Returns
610 -------
611 I : array_like
612 Definite integral of the spline over the interval ``[a, b]``.
614 Examples
615 --------
616 Construct the linear spline ``x if x < 1 else 2 - x`` on the base
617 interval :math:`[0, 2]`, and integrate it
619 >>> from scipy.interpolate import BSpline
620 >>> b = BSpline.basis_element([0, 1, 2])
621 >>> b.integrate(0, 1)
622 array(0.5)
624 If the integration limits are outside of the base interval, the result
625 is controlled by the `extrapolate` parameter
627 >>> b.integrate(-1, 1)
628 array(0.0)
629 >>> b.integrate(-1, 1, extrapolate=False)
630 array(0.5)
632 >>> import matplotlib.pyplot as plt
633 >>> fig, ax = plt.subplots()
634 >>> ax.grid(True)
635 >>> ax.axvline(0, c='r', lw=5, alpha=0.5) # base interval
636 >>> ax.axvline(2, c='r', lw=5, alpha=0.5)
637 >>> xx = [-1, 1, 2]
638 >>> ax.plot(xx, b(xx))
639 >>> plt.show()
641 """
642 if extrapolate is None:
643 extrapolate = self.extrapolate
645 # Prepare self.t and self.c.
646 self._ensure_c_contiguous()
648 # Swap integration bounds if needed.
649 sign = 1
650 if b < a:
651 a, b = b, a
652 sign = -1
653 n = self.t.size - self.k - 1
655 if extrapolate != "periodic" and not extrapolate:
656 # Shrink the integration interval, if needed.
657 a = max(a, self.t[self.k])
658 b = min(b, self.t[n])
660 if self.c.ndim == 1:
661 # Fast path: use FITPACK's routine
662 # (cf _fitpack_impl.splint).
663 integral = _fitpack_impl.splint(a, b, self.tck)
664 return integral * sign
666 out = np.empty((2, prod(self.c.shape[1:])), dtype=self.c.dtype)
668 # Compute the antiderivative.
669 c = self.c
670 ct = len(self.t) - len(c)
671 if ct > 0:
672 c = np.r_[c, np.zeros((ct,) + c.shape[1:])]
673 ta, ca, ka = _fitpack_impl.splantider((self.t, c, self.k), 1)
675 if extrapolate == 'periodic':
676 # Split the integral into the part over period (can be several
677 # of them) and the remaining part.
679 ts, te = self.t[self.k], self.t[n]
680 period = te - ts
681 interval = b - a
682 n_periods, left = divmod(interval, period)
684 if n_periods > 0:
685 # Evaluate the difference of antiderivatives.
686 x = np.asarray([ts, te], dtype=np.float_)
687 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1),
688 ka, x, 0, False, out)
689 integral = out[1] - out[0]
690 integral *= n_periods
691 else:
692 integral = np.zeros((1, prod(self.c.shape[1:])),
693 dtype=self.c.dtype)
695 # Map a to [ts, te], b is always a + left.
696 a = ts + (a - ts) % period
697 b = a + left
699 # If b <= te then we need to integrate over [a, b], otherwise
700 # over [a, te] and from xs to what is remained.
701 if b <= te:
702 x = np.asarray([a, b], dtype=np.float_)
703 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1),
704 ka, x, 0, False, out)
705 integral += out[1] - out[0]
706 else:
707 x = np.asarray([a, te], dtype=np.float_)
708 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1),
709 ka, x, 0, False, out)
710 integral += out[1] - out[0]
712 x = np.asarray([ts, ts + b - te], dtype=np.float_)
713 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1),
714 ka, x, 0, False, out)
715 integral += out[1] - out[0]
716 else:
717 # Evaluate the difference of antiderivatives.
718 x = np.asarray([a, b], dtype=np.float_)
719 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1),
720 ka, x, 0, extrapolate, out)
721 integral = out[1] - out[0]
723 integral *= sign
724 return integral.reshape(ca.shape[1:])
726 @classmethod
727 def from_power_basis(cls, pp, bc_type='not-a-knot'):
728 r"""
729 Construct a polynomial in the B-spline basis
730 from a piecewise polynomial in the power basis.
732 For now, accepts ``CubicSpline`` instances only.
734 Parameters
735 ----------
736 pp : CubicSpline
737 A piecewise polynomial in the power basis, as created
738 by ``CubicSpline``
739 bc_type : string, optional
740 Boundary condition type as in ``CubicSpline``: one of the
741 ``not-a-knot``, ``natural``, ``clamped``, or ``periodic``.
742 Necessary for construction an instance of ``BSpline`` class.
743 Default is ``not-a-knot``.
745 Returns
746 -------
747 b : BSpline object
748 A new instance representing the initial polynomial
749 in the B-spline basis.
751 Notes
752 -----
753 .. versionadded:: 1.8.0
755 Accepts only ``CubicSpline`` instances for now.
757 The algorithm follows from differentiation
758 the Marsden's identity [1]: each of coefficients of spline
759 interpolation function in the B-spline basis is computed as follows:
761 .. math::
763 c_j = \sum_{m=0}^{k} \frac{(k-m)!}{k!}
764 c_{m,i} (-1)^{k-m} D^m p_{j,k}(x_i)
766 :math:`c_{m, i}` - a coefficient of CubicSpline,
767 :math:`D^m p_{j, k}(x_i)` - an m-th defivative of a dual polynomial
768 in :math:`x_i`.
770 ``k`` always equals 3 for now.
772 First ``n - 2`` coefficients are computed in :math:`x_i = x_j`, e.g.
774 .. math::
776 c_1 = \sum_{m=0}^{k} \frac{(k-1)!}{k!} c_{m,1} D^m p_{j,3}(x_1)
778 Last ``nod + 2`` coefficients are computed in ``x[-2]``,
779 ``nod`` - number of derivatives at the ends.
781 For example, consider :math:`x = [0, 1, 2, 3, 4]`,
782 :math:`y = [1, 1, 1, 1, 1]` and bc_type = ``natural``
784 The coefficients of CubicSpline in the power basis:
786 :math:`[[0, 0, 0, 0, 0], [0, 0, 0, 0, 0],
787 [0, 0, 0, 0, 0], [1, 1, 1, 1, 1]]`
789 The knot vector: :math:`t = [0, 0, 0, 0, 1, 2, 3, 4, 4, 4, 4]`
791 In this case
793 .. math::
795 c_j = \frac{0!}{k!} c_{3, i} k! = c_{3, i} = 1,~j = 0, ..., 6
797 References
798 ----------
799 .. [1] Tom Lyche and Knut Morken, Spline Methods, 2005, Section 3.1.2
801 """
802 from ._cubic import CubicSpline
803 if not isinstance(pp, CubicSpline):
804 raise NotImplementedError("Only CubicSpline objects are accepted"
805 "for now. Got %s instead." % type(pp))
806 x = pp.x
807 coef = pp.c
808 k = pp.c.shape[0] - 1
809 n = x.shape[0]
811 if bc_type == 'not-a-knot':
812 t = _not_a_knot(x, k)
813 elif bc_type == 'natural' or bc_type == 'clamped':
814 t = _augknt(x, k)
815 elif bc_type == 'periodic':
816 t = _periodic_knots(x, k)
817 else:
818 raise TypeError('Unknown boundary condition: %s' % bc_type)
820 nod = t.shape[0] - (n + k + 1) # number of derivatives at the ends
821 c = np.zeros(n + nod, dtype=pp.c.dtype)
822 for m in range(k + 1):
823 for i in range(n - 2):
824 c[i] += poch(k + 1, -m) * coef[m, i]\
825 * np.power(-1, k - m)\
826 * _diff_dual_poly(i, k, x[i], m, t)
827 for j in range(n - 2, n + nod):
828 c[j] += poch(k + 1, -m) * coef[m, n - 2]\
829 * np.power(-1, k - m)\
830 * _diff_dual_poly(j, k, x[n - 2], m, t)
831 return cls.construct_fast(t, c, k, pp.extrapolate, pp.axis)
834#################################
835# Interpolating spline helpers #
836#################################
838def _not_a_knot(x, k):
839 """Given data x, construct the knot vector w/ not-a-knot BC.
840 cf de Boor, XIII(12)."""
841 x = np.asarray(x)
842 if k % 2 != 1:
843 raise ValueError("Odd degree for now only. Got %s." % k)
845 m = (k - 1) // 2
846 t = x[m+1:-m-1]
847 t = np.r_[(x[0],)*(k+1), t, (x[-1],)*(k+1)]
848 return t
851def _augknt(x, k):
852 """Construct a knot vector appropriate for the order-k interpolation."""
853 return np.r_[(x[0],)*k, x, (x[-1],)*k]
856def _convert_string_aliases(deriv, target_shape):
857 if isinstance(deriv, str):
858 if deriv == "clamped":
859 deriv = [(1, np.zeros(target_shape))]
860 elif deriv == "natural":
861 deriv = [(2, np.zeros(target_shape))]
862 else:
863 raise ValueError("Unknown boundary condition : %s" % deriv)
864 return deriv
867def _process_deriv_spec(deriv):
868 if deriv is not None:
869 try:
870 ords, vals = zip(*deriv)
871 except TypeError as e:
872 msg = ("Derivatives, `bc_type`, should be specified as a pair of "
873 "iterables of pairs of (order, value).")
874 raise ValueError(msg) from e
875 else:
876 ords, vals = [], []
877 return np.atleast_1d(ords, vals)
880def _woodbury_algorithm(A, ur, ll, b, k):
881 '''
882 Solve a cyclic banded linear system with upper right
883 and lower blocks of size ``(k-1) / 2`` using
884 the Woodbury formula
886 Parameters
887 ----------
888 A : 2-D array, shape(k, n)
889 Matrix of diagonals of original matrix (see
890 ``solve_banded`` documentation).
891 ur : 2-D array, shape(bs, bs)
892 Upper right block matrix.
893 ll : 2-D array, shape(bs, bs)
894 Lower left block matrix.
895 b : 1-D array, shape(n,)
896 Vector of constant terms of the system of linear equations.
897 k : int
898 B-spline degree.
900 Returns
901 -------
902 c : 1-D array, shape(n,)
903 Solution of the original system of linear equations.
905 Notes
906 -----
907 This algorithm works only for systems with banded matrix A plus
908 a correction term U @ V.T, where the matrix U @ V.T gives upper right
909 and lower left block of A
910 The system is solved with the following steps:
911 1. New systems of linear equations are constructed:
912 A @ z_i = u_i,
913 u_i - columnn vector of U,
914 i = 1, ..., k - 1
915 2. Matrix Z is formed from vectors z_i:
916 Z = [ z_1 | z_2 | ... | z_{k - 1} ]
917 3. Matrix H = (1 + V.T @ Z)^{-1}
918 4. The system A' @ y = b is solved
919 5. x = y - Z @ (H @ V.T @ y)
920 Also, ``n`` should be greater than ``k``, otherwise corner block
921 elements will intersect with diagonals.
923 Examples
924 --------
925 Consider the case of n = 8, k = 5 (size of blocks - 2 x 2).
926 The matrix of a system: U: V:
927 x x x * * a b a b 0 0 0 0 1 0
928 x x x x * * c 0 c 0 0 0 0 0 1
929 x x x x x * * 0 0 0 0 0 0 0 0
930 * x x x x x * 0 0 0 0 0 0 0 0
931 * * x x x x x 0 0 0 0 0 0 0 0
932 d * * x x x x 0 0 d 0 1 0 0 0
933 e f * * x x x 0 0 e f 0 1 0 0
935 References
936 ----------
937 .. [1] William H. Press, Saul A. Teukolsky, William T. Vetterling
938 and Brian P. Flannery, Numerical Recipes, 2007, Section 2.7.3
940 '''
941 k_mod = k - k % 2
942 bs = int((k - 1) / 2) + (k + 1) % 2
944 n = A.shape[1] + 1
945 U = np.zeros((n - 1, k_mod))
946 VT = np.zeros((k_mod, n - 1)) # V transpose
948 # upper right block
949 U[:bs, :bs] = ur
950 VT[np.arange(bs), np.arange(bs) - bs] = 1
952 # lower left block
953 U[-bs:, -bs:] = ll
954 VT[np.arange(bs) - bs, np.arange(bs)] = 1
956 Z = solve_banded((bs, bs), A, U)
958 H = solve(np.identity(k_mod) + VT @ Z, np.identity(k_mod))
960 y = solve_banded((bs, bs), A, b)
961 c = y - Z @ (H @ (VT @ y))
963 return c
966def _periodic_knots(x, k):
967 '''
968 returns vector of nodes on circle
969 '''
970 xc = np.copy(x)
971 n = len(xc)
972 if k % 2 == 0:
973 dx = np.diff(xc)
974 xc[1: -1] -= dx[:-1] / 2
975 dx = np.diff(xc)
976 t = np.zeros(n + 2 * k)
977 t[k: -k] = xc
978 for i in range(0, k):
979 # filling first `k` elements in descending order
980 t[k - i - 1] = t[k - i] - dx[-(i % (n - 1)) - 1]
981 # filling last `k` elements in ascending order
982 t[-k + i] = t[-k + i - 1] + dx[i % (n - 1)]
983 return t
986def _make_interp_per_full_matr(x, y, t, k):
987 '''
988 Returns a solution of a system for B-spline interpolation with periodic
989 boundary conditions. First ``k - 1`` rows of matrix are condtions of
990 periodicity (continuity of ``k - 1`` derivatives at the boundary points).
991 Last ``n`` rows are interpolation conditions.
992 RHS is ``k - 1`` zeros and ``n`` ordinates in this case.
994 Parameters
995 ----------
996 x : 1-D array, shape (n,)
997 Values of x - coordinate of a given set of points.
998 y : 1-D array, shape (n,)
999 Values of y - coordinate of a given set of points.
1000 t : 1-D array, shape(n+2*k,)
1001 Vector of knots.
1002 k : int
1003 The maximum degree of spline
1005 Returns
1006 -------
1007 c : 1-D array, shape (n+k-1,)
1008 B-spline coefficients
1010 Notes
1011 -----
1012 ``t`` is supposed to be taken on circle.
1014 '''
1016 x, y, t = map(np.asarray, (x, y, t))
1018 n = x.size
1019 # LHS: the collocation matrix + derivatives at edges
1020 matr = np.zeros((n + k - 1, n + k - 1))
1022 # derivatives at x[0] and x[-1]:
1023 for i in range(k - 1):
1024 bb = _bspl.evaluate_all_bspl(t, k, x[0], k, nu=i + 1)
1025 matr[i, : k + 1] += bb
1026 bb = _bspl.evaluate_all_bspl(t, k, x[-1], n + k - 1, nu=i + 1)[:-1]
1027 matr[i, -k:] -= bb
1029 # collocation matrix
1030 for i in range(n):
1031 xval = x[i]
1032 # find interval
1033 if xval == t[k]:
1034 left = k
1035 else:
1036 left = np.searchsorted(t, xval) - 1
1038 # fill a row
1039 bb = _bspl.evaluate_all_bspl(t, k, xval, left)
1040 matr[i + k - 1, left-k:left+1] = bb
1042 # RHS
1043 b = np.r_[[0] * (k - 1), y]
1045 c = solve(matr, b)
1046 return c
1049def _make_periodic_spline(x, y, t, k, axis):
1050 '''
1051 Compute the (coefficients of) interpolating B-spline with periodic
1052 boundary conditions.
1054 Parameters
1055 ----------
1056 x : array_like, shape (n,)
1057 Abscissas.
1058 y : array_like, shape (n,)
1059 Ordinates.
1060 k : int
1061 B-spline degree.
1062 t : array_like, shape (n + 2 * k,).
1063 Knots taken on a circle, ``k`` on the left and ``k`` on the right
1064 of the vector ``x``.
1066 Returns
1067 -------
1068 b : a BSpline object of the degree ``k`` and with knots ``t``.
1070 Notes
1071 -----
1072 The original system is formed by ``n + k - 1`` equations where the first
1073 ``k - 1`` of them stand for the ``k - 1`` derivatives continuity on the
1074 edges while the other equations correspond to an interpolating case
1075 (matching all the input points). Due to a special form of knot vector, it
1076 can be proved that in the original system the first and last ``k``
1077 coefficients of a spline function are the same, respectively. It follows
1078 from the fact that all ``k - 1`` derivatives are equal term by term at ends
1079 and that the matrix of the original system of linear equations is
1080 non-degenerate. So, we can reduce the number of equations to ``n - 1``
1081 (first ``k - 1`` equations could be reduced). Another trick of this
1082 implementation is cyclic shift of values of B-splines due to equality of
1083 ``k`` unknown coefficients. With this we can receive matrix of the system
1084 with upper right and lower left blocks, and ``k`` diagonals. It allows
1085 to use Woodbury formula to optimize the computations.
1087 '''
1088 n = y.shape[0]
1090 extradim = prod(y.shape[1:])
1091 y_new = y.reshape(n, extradim)
1092 c = np.zeros((n + k - 1, extradim))
1094 # n <= k case is solved with full matrix
1095 if n <= k:
1096 for i in range(extradim):
1097 c[:, i] = _make_interp_per_full_matr(x, y_new[:, i], t, k)
1098 c = np.ascontiguousarray(c.reshape((n + k - 1,) + y.shape[1:]))
1099 return BSpline.construct_fast(t, c, k, extrapolate='periodic', axis=axis)
1101 nt = len(t) - k - 1
1103 # size of block elements
1104 kul = int(k / 2)
1106 # kl = ku = k
1107 ab = np.zeros((3 * k + 1, nt), dtype=np.float_, order='F')
1109 # upper right and lower left blocks
1110 ur = np.zeros((kul, kul))
1111 ll = np.zeros_like(ur)
1113 # `offset` is made to shift all the non-zero elements to the end of the
1114 # matrix
1115 _bspl._colloc(x, t, k, ab, offset=k)
1117 # remove zeros before the matrix
1118 ab = ab[-k - (k + 1) % 2:, :]
1120 # The least elements in rows (except repetitions) are diagonals
1121 # of block matrices. Upper right matrix is an upper triangular
1122 # matrix while lower left is a lower triangular one.
1123 for i in range(kul):
1124 ur += np.diag(ab[-i - 1, i: kul], k=i)
1125 ll += np.diag(ab[i, -kul - (k % 2): n - 1 + 2 * kul - i], k=-i)
1127 # remove elements that occur in the last point
1128 # (first and last points are equivalent)
1129 A = ab[:, kul: -k + kul]
1131 for i in range(extradim):
1132 cc = _woodbury_algorithm(A, ur, ll, y_new[:, i][:-1], k)
1133 c[:, i] = np.concatenate((cc[-kul:], cc, cc[:kul + k % 2]))
1134 c = np.ascontiguousarray(c.reshape((n + k - 1,) + y.shape[1:]))
1135 return BSpline.construct_fast(t, c, k, extrapolate='periodic', axis=axis)
1138def make_interp_spline(x, y, k=3, t=None, bc_type=None, axis=0,
1139 check_finite=True):
1140 """Compute the (coefficients of) interpolating B-spline.
1142 Parameters
1143 ----------
1144 x : array_like, shape (n,)
1145 Abscissas.
1146 y : array_like, shape (n, ...)
1147 Ordinates.
1148 k : int, optional
1149 B-spline degree. Default is cubic, ``k = 3``.
1150 t : array_like, shape (nt + k + 1,), optional.
1151 Knots.
1152 The number of knots needs to agree with the number of data points and
1153 the number of derivatives at the edges. Specifically, ``nt - n`` must
1154 equal ``len(deriv_l) + len(deriv_r)``.
1155 bc_type : 2-tuple or None
1156 Boundary conditions.
1157 Default is None, which means choosing the boundary conditions
1158 automatically. Otherwise, it must be a length-two tuple where the first
1159 element (``deriv_l``) sets the boundary conditions at ``x[0]`` and
1160 the second element (``deriv_r``) sets the boundary conditions at
1161 ``x[-1]``. Each of these must be an iterable of pairs
1162 ``(order, value)`` which gives the values of derivatives of specified
1163 orders at the given edge of the interpolation interval.
1164 Alternatively, the following string aliases are recognized:
1166 * ``"clamped"``: The first derivatives at the ends are zero. This is
1167 equivalent to ``bc_type=([(1, 0.0)], [(1, 0.0)])``.
1168 * ``"natural"``: The second derivatives at ends are zero. This is
1169 equivalent to ``bc_type=([(2, 0.0)], [(2, 0.0)])``.
1170 * ``"not-a-knot"`` (default): The first and second segments are the
1171 same polynomial. This is equivalent to having ``bc_type=None``.
1172 * ``"periodic"``: The values and the first ``k-1`` derivatives at the
1173 ends are equivalent.
1175 axis : int, optional
1176 Interpolation axis. Default is 0.
1177 check_finite : bool, optional
1178 Whether to check that the input arrays contain only finite numbers.
1179 Disabling may give a performance gain, but may result in problems
1180 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1181 Default is True.
1183 Returns
1184 -------
1185 b : a BSpline object of the degree ``k`` and with knots ``t``.
1187 See Also
1188 --------
1189 BSpline : base class representing the B-spline objects
1190 CubicSpline : a cubic spline in the polynomial basis
1191 make_lsq_spline : a similar factory function for spline fitting
1192 UnivariateSpline : a wrapper over FITPACK spline fitting routines
1193 splrep : a wrapper over FITPACK spline fitting routines
1195 Examples
1196 --------
1198 Use cubic interpolation on Chebyshev nodes:
1200 >>> import numpy as np
1201 >>> import matplotlib.pyplot as plt
1202 >>> def cheb_nodes(N):
1203 ... jj = 2.*np.arange(N) + 1
1204 ... x = np.cos(np.pi * jj / 2 / N)[::-1]
1205 ... return x
1207 >>> x = cheb_nodes(20)
1208 >>> y = np.sqrt(1 - x**2)
1210 >>> from scipy.interpolate import BSpline, make_interp_spline
1211 >>> b = make_interp_spline(x, y)
1212 >>> np.allclose(b(x), y)
1213 True
1215 Note that the default is a cubic spline with a not-a-knot boundary condition
1217 >>> b.k
1218 3
1220 Here we use a 'natural' spline, with zero 2nd derivatives at edges:
1222 >>> l, r = [(2, 0.0)], [(2, 0.0)]
1223 >>> b_n = make_interp_spline(x, y, bc_type=(l, r)) # or, bc_type="natural"
1224 >>> np.allclose(b_n(x), y)
1225 True
1226 >>> x0, x1 = x[0], x[-1]
1227 >>> np.allclose([b_n(x0, 2), b_n(x1, 2)], [0, 0])
1228 True
1230 Interpolation of parametric curves is also supported. As an example, we
1231 compute a discretization of a snail curve in polar coordinates
1233 >>> phi = np.linspace(0, 2.*np.pi, 40)
1234 >>> r = 0.3 + np.cos(phi)
1235 >>> x, y = r*np.cos(phi), r*np.sin(phi) # convert to Cartesian coordinates
1237 Build an interpolating curve, parameterizing it by the angle
1239 >>> spl = make_interp_spline(phi, np.c_[x, y])
1241 Evaluate the interpolant on a finer grid (note that we transpose the result
1242 to unpack it into a pair of x- and y-arrays)
1244 >>> phi_new = np.linspace(0, 2.*np.pi, 100)
1245 >>> x_new, y_new = spl(phi_new).T
1247 Plot the result
1249 >>> plt.plot(x, y, 'o')
1250 >>> plt.plot(x_new, y_new, '-')
1251 >>> plt.show()
1253 Build a B-spline curve with 2 dimensional y
1255 >>> x = np.linspace(0, 2*np.pi, 10)
1256 >>> y = np.array([np.sin(x), np.cos(x)])
1258 Periodic condition is satisfied because y coordinates of points on the ends
1259 are equivalent
1261 >>> ax = plt.axes(projection='3d')
1262 >>> xx = np.linspace(0, 2*np.pi, 100)
1263 >>> bspl = make_interp_spline(x, y, k=5, bc_type='periodic', axis=1)
1264 >>> ax.plot3D(xx, *bspl(xx))
1265 >>> ax.scatter3D(x, *y, color='red')
1266 >>> plt.show()
1268 """
1269 # convert string aliases for the boundary conditions
1270 if bc_type is None or bc_type == 'not-a-knot' or bc_type == 'periodic':
1271 deriv_l, deriv_r = None, None
1272 elif isinstance(bc_type, str):
1273 deriv_l, deriv_r = bc_type, bc_type
1274 else:
1275 try:
1276 deriv_l, deriv_r = bc_type
1277 except TypeError as e:
1278 raise ValueError("Unknown boundary condition: %s" % bc_type) from e
1280 y = np.asarray(y)
1282 axis = normalize_axis_index(axis, y.ndim)
1284 x = _as_float_array(x, check_finite)
1285 y = _as_float_array(y, check_finite)
1287 y = np.moveaxis(y, axis, 0) # now internally interp axis is zero
1289 # sanity check the input
1290 if bc_type == 'periodic' and not np.allclose(y[0], y[-1], atol=1e-15):
1291 raise ValueError("First and last points does not match while "
1292 "periodic case expected")
1293 if x.size != y.shape[0]:
1294 raise ValueError('Shapes of x {} and y {} are incompatible'
1295 .format(x.shape, y.shape))
1296 if np.any(x[1:] == x[:-1]):
1297 raise ValueError("Expect x to not have duplicates")
1298 if x.ndim != 1 or np.any(x[1:] < x[:-1]):
1299 raise ValueError("Expect x to be a 1D strictly increasing sequence.")
1301 # special-case k=0 right away
1302 if k == 0:
1303 if any(_ is not None for _ in (t, deriv_l, deriv_r)):
1304 raise ValueError("Too much info for k=0: t and bc_type can only "
1305 "be None.")
1306 t = np.r_[x, x[-1]]
1307 c = np.asarray(y)
1308 c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
1309 return BSpline.construct_fast(t, c, k, axis=axis)
1311 # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
1312 if k == 1 and t is None:
1313 if not (deriv_l is None and deriv_r is None):
1314 raise ValueError("Too much info for k=1: bc_type can only be None.")
1315 t = np.r_[x[0], x, x[-1]]
1316 c = np.asarray(y)
1317 c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
1318 return BSpline.construct_fast(t, c, k, axis=axis)
1320 k = operator.index(k)
1322 if bc_type == 'periodic' and t is not None:
1323 raise NotImplementedError("For periodic case t is constructed "
1324 "automatically and can not be passed "
1325 "manually")
1327 # come up with a sensible knot vector, if needed
1328 if t is None:
1329 if deriv_l is None and deriv_r is None:
1330 if bc_type == 'periodic':
1331 t = _periodic_knots(x, k)
1332 elif k == 2:
1333 # OK, it's a bit ad hoc: Greville sites + omit
1334 # 2nd and 2nd-to-last points, a la not-a-knot
1335 t = (x[1:] + x[:-1]) / 2.
1336 t = np.r_[(x[0],)*(k+1),
1337 t[1:-1],
1338 (x[-1],)*(k+1)]
1339 else:
1340 t = _not_a_knot(x, k)
1341 else:
1342 t = _augknt(x, k)
1344 t = _as_float_array(t, check_finite)
1346 if k < 0:
1347 raise ValueError("Expect non-negative k.")
1348 if t.ndim != 1 or np.any(t[1:] < t[:-1]):
1349 raise ValueError("Expect t to be a 1-D sorted array_like.")
1350 if t.size < x.size + k + 1:
1351 raise ValueError('Got %d knots, need at least %d.' %
1352 (t.size, x.size + k + 1))
1353 if (x[0] < t[k]) or (x[-1] > t[-k]):
1354 raise ValueError('Out of bounds w/ x = %s.' % x)
1356 if bc_type == 'periodic':
1357 return _make_periodic_spline(x, y, t, k, axis)
1359 # Here : deriv_l, r = [(nu, value), ...]
1360 deriv_l = _convert_string_aliases(deriv_l, y.shape[1:])
1361 deriv_l_ords, deriv_l_vals = _process_deriv_spec(deriv_l)
1362 nleft = deriv_l_ords.shape[0]
1364 deriv_r = _convert_string_aliases(deriv_r, y.shape[1:])
1365 deriv_r_ords, deriv_r_vals = _process_deriv_spec(deriv_r)
1366 nright = deriv_r_ords.shape[0]
1368 # have `n` conditions for `nt` coefficients; need nt-n derivatives
1369 n = x.size
1370 nt = t.size - k - 1
1372 if nt - n != nleft + nright:
1373 raise ValueError("The number of derivatives at boundaries does not "
1374 "match: expected {}, got {}+{}".format(nt-n, nleft, nright))
1376 # bail out if the `y` array is zero-sized
1377 if y.size == 0:
1378 c = np.zeros((nt,) + y.shape[1:], dtype=float)
1379 return BSpline.construct_fast(t, c, k, axis=axis)
1381 # set up the LHS: the collocation matrix + derivatives at boundaries
1382 kl = ku = k
1383 ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
1384 _bspl._colloc(x, t, k, ab, offset=nleft)
1385 if nleft > 0:
1386 _bspl._handle_lhs_derivatives(t, k, x[0], ab, kl, ku, deriv_l_ords)
1387 if nright > 0:
1388 _bspl._handle_lhs_derivatives(t, k, x[-1], ab, kl, ku, deriv_r_ords,
1389 offset=nt-nright)
1391 # set up the RHS: values to interpolate (+ derivative values, if any)
1392 extradim = prod(y.shape[1:])
1393 rhs = np.empty((nt, extradim), dtype=y.dtype)
1394 if nleft > 0:
1395 rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
1396 rhs[nleft:nt - nright] = y.reshape(-1, extradim)
1397 if nright > 0:
1398 rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
1400 # solve Ab @ x = rhs; this is the relevant part of linalg.solve_banded
1401 if check_finite:
1402 ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
1403 gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
1404 lu, piv, c, info = gbsv(kl, ku, ab, rhs,
1405 overwrite_ab=True, overwrite_b=True)
1407 if info > 0:
1408 raise LinAlgError("Collocation matrix is singular.")
1409 elif info < 0:
1410 raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
1412 c = np.ascontiguousarray(c.reshape((nt,) + y.shape[1:]))
1413 return BSpline.construct_fast(t, c, k, axis=axis)
1416def make_lsq_spline(x, y, t, k=3, w=None, axis=0, check_finite=True):
1417 r"""Compute the (coefficients of) an LSQ (Least SQuared) based
1418 fitting B-spline.
1420 The result is a linear combination
1422 .. math::
1424 S(x) = \sum_j c_j B_j(x; t)
1426 of the B-spline basis elements, :math:`B_j(x; t)`, which minimizes
1428 .. math::
1430 \sum_{j} \left( w_j \times (S(x_j) - y_j) \right)^2
1432 Parameters
1433 ----------
1434 x : array_like, shape (m,)
1435 Abscissas.
1436 y : array_like, shape (m, ...)
1437 Ordinates.
1438 t : array_like, shape (n + k + 1,).
1439 Knots.
1440 Knots and data points must satisfy Schoenberg-Whitney conditions.
1441 k : int, optional
1442 B-spline degree. Default is cubic, ``k = 3``.
1443 w : array_like, shape (m,), optional
1444 Weights for spline fitting. Must be positive. If ``None``,
1445 then weights are all equal.
1446 Default is ``None``.
1447 axis : int, optional
1448 Interpolation axis. Default is zero.
1449 check_finite : bool, optional
1450 Whether to check that the input arrays contain only finite numbers.
1451 Disabling may give a performance gain, but may result in problems
1452 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1453 Default is True.
1455 Returns
1456 -------
1457 b : a BSpline object of the degree ``k`` with knots ``t``.
1459 See Also
1460 --------
1461 BSpline : base class representing the B-spline objects
1462 make_interp_spline : a similar factory function for interpolating splines
1463 LSQUnivariateSpline : a FITPACK-based spline fitting routine
1464 splrep : a FITPACK-based fitting routine
1466 Notes
1467 -----
1468 The number of data points must be larger than the spline degree ``k``.
1470 Knots ``t`` must satisfy the Schoenberg-Whitney conditions,
1471 i.e., there must be a subset of data points ``x[j]`` such that
1472 ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``.
1474 Examples
1475 --------
1476 Generate some noisy data:
1478 >>> import numpy as np
1479 >>> import matplotlib.pyplot as plt
1480 >>> rng = np.random.default_rng()
1481 >>> x = np.linspace(-3, 3, 50)
1482 >>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(50)
1484 Now fit a smoothing cubic spline with a pre-defined internal knots.
1485 Here we make the knot vector (k+1)-regular by adding boundary knots:
1487 >>> from scipy.interpolate import make_lsq_spline, BSpline
1488 >>> t = [-1, 0, 1]
1489 >>> k = 3
1490 >>> t = np.r_[(x[0],)*(k+1),
1491 ... t,
1492 ... (x[-1],)*(k+1)]
1493 >>> spl = make_lsq_spline(x, y, t, k)
1495 For comparison, we also construct an interpolating spline for the same
1496 set of data:
1498 >>> from scipy.interpolate import make_interp_spline
1499 >>> spl_i = make_interp_spline(x, y)
1501 Plot both:
1503 >>> xs = np.linspace(-3, 3, 100)
1504 >>> plt.plot(x, y, 'ro', ms=5)
1505 >>> plt.plot(xs, spl(xs), 'g-', lw=3, label='LSQ spline')
1506 >>> plt.plot(xs, spl_i(xs), 'b-', lw=3, alpha=0.7, label='interp spline')
1507 >>> plt.legend(loc='best')
1508 >>> plt.show()
1510 **NaN handling**: If the input arrays contain ``nan`` values, the result is
1511 not useful since the underlying spline fitting routines cannot deal with
1512 ``nan``. A workaround is to use zero weights for not-a-number data points:
1514 >>> y[8] = np.nan
1515 >>> w = np.isnan(y)
1516 >>> y[w] = 0.
1517 >>> tck = make_lsq_spline(x, y, t, w=~w)
1519 Notice the need to replace a ``nan`` by a numerical value (precise value
1520 does not matter as long as the corresponding weight is zero.)
1522 """
1523 x = _as_float_array(x, check_finite)
1524 y = _as_float_array(y, check_finite)
1525 t = _as_float_array(t, check_finite)
1526 if w is not None:
1527 w = _as_float_array(w, check_finite)
1528 else:
1529 w = np.ones_like(x)
1530 k = operator.index(k)
1532 axis = normalize_axis_index(axis, y.ndim)
1534 y = np.moveaxis(y, axis, 0) # now internally interp axis is zero
1536 if x.ndim != 1 or np.any(x[1:] - x[:-1] <= 0):
1537 raise ValueError("Expect x to be a 1-D sorted array_like.")
1538 if x.shape[0] < k+1:
1539 raise ValueError("Need more x points.")
1540 if k < 0:
1541 raise ValueError("Expect non-negative k.")
1542 if t.ndim != 1 or np.any(t[1:] - t[:-1] < 0):
1543 raise ValueError("Expect t to be a 1-D sorted array_like.")
1544 if x.size != y.shape[0]:
1545 raise ValueError('Shapes of x {} and y {} are incompatible'
1546 .format(x.shape, y.shape))
1547 if k > 0 and np.any((x < t[k]) | (x > t[-k])):
1548 raise ValueError('Out of bounds w/ x = %s.' % x)
1549 if x.size != w.size:
1550 raise ValueError('Shapes of x {} and w {} are incompatible'
1551 .format(x.shape, w.shape))
1553 # number of coefficients
1554 n = t.size - k - 1
1556 # construct A.T @ A and rhs with A the collocation matrix, and
1557 # rhs = A.T @ y for solving the LSQ problem ``A.T @ A @ c = A.T @ y``
1558 lower = True
1559 extradim = prod(y.shape[1:])
1560 ab = np.zeros((k+1, n), dtype=np.float_, order='F')
1561 rhs = np.zeros((n, extradim), dtype=y.dtype, order='F')
1562 _bspl._norm_eq_lsq(x, t, k,
1563 y.reshape(-1, extradim),
1564 w,
1565 ab, rhs)
1566 rhs = rhs.reshape((n,) + y.shape[1:])
1568 # have observation matrix & rhs, can solve the LSQ problem
1569 cho_decomp = cholesky_banded(ab, overwrite_ab=True, lower=lower,
1570 check_finite=check_finite)
1571 c = cho_solve_banded((cho_decomp, lower), rhs, overwrite_b=True,
1572 check_finite=check_finite)
1574 c = np.ascontiguousarray(c)
1575 return BSpline.construct_fast(t, c, k, axis=axis)
1578#############################
1579# Smoothing spline helpers #
1580#############################
1582def _compute_optimal_gcv_parameter(X, wE, y, w):
1583 """
1584 Returns an optimal regularization parameter from the GCV criteria [1].
1586 Parameters
1587 ----------
1588 X : array, shape (5, n)
1589 5 bands of the design matrix ``X`` stored in LAPACK banded storage.
1590 wE : array, shape (5, n)
1591 5 bands of the penalty matrix :math:`W^{-1} E` stored in LAPACK banded
1592 storage.
1593 y : array, shape (n,)
1594 Ordinates.
1595 w : array, shape (n,)
1596 Vector of weights.
1598 Returns
1599 -------
1600 lam : float
1601 An optimal from the GCV criteria point of view regularization
1602 parameter.
1604 Notes
1605 -----
1606 No checks are performed.
1608 References
1609 ----------
1610 .. [1] G. Wahba, "Estimating the smoothing parameter" in Spline models
1611 for observational data, Philadelphia, Pennsylvania: Society for
1612 Industrial and Applied Mathematics, 1990, pp. 45-65.
1613 :doi:`10.1137/1.9781611970128`
1615 """
1617 def compute_banded_symmetric_XT_W_Y(X, w, Y):
1618 """
1619 Assuming that the product :math:`X^T W Y` is symmetric and both ``X``
1620 and ``Y`` are 5-banded, compute the unique bands of the product.
1622 Parameters
1623 ----------
1624 X : array, shape (5, n)
1625 5 bands of the matrix ``X`` stored in LAPACK banded storage.
1626 w : array, shape (n,)
1627 Array of weights
1628 Y : array, shape (5, n)
1629 5 bands of the matrix ``Y`` stored in LAPACK banded storage.
1631 Returns
1632 -------
1633 res : array, shape (4, n)
1634 The result of the product :math:`X^T Y` stored in the banded way.
1636 Notes
1637 -----
1638 As far as the matrices ``X`` and ``Y`` are 5-banded, their product
1639 :math:`X^T W Y` is 7-banded. It is also symmetric, so we can store only
1640 unique diagonals.
1642 """
1643 # compute W Y
1644 W_Y = np.copy(Y)
1646 W_Y[2] *= w
1647 for i in range(2):
1648 W_Y[i, 2 - i:] *= w[:-2 + i]
1649 W_Y[3 + i, :-1 - i] *= w[1 + i:]
1651 n = X.shape[1]
1652 res = np.zeros((4, n))
1653 for i in range(n):
1654 for j in range(min(n-i, 4)):
1655 res[-j-1, i + j] = sum(X[j:, i] * W_Y[:5-j, i + j])
1656 return res
1658 def compute_b_inv(A):
1659 """
1660 Inverse 3 central bands of matrix :math:`A=U^T D^{-1} U` assuming that
1661 ``U`` is a unit upper triangular banded matrix using an algorithm
1662 proposed in [1].
1664 Parameters
1665 ----------
1666 A : array, shape (4, n)
1667 Matrix to inverse, stored in LAPACK banded storage.
1669 Returns
1670 -------
1671 B : array, shape (4, n)
1672 3 unique bands of the symmetric matrix that is an inverse to ``A``.
1673 The first row is filled with zeros.
1675 Notes
1676 -----
1677 The algorithm is based on the cholesky decomposition and, therefore,
1678 in case matrix ``A`` is close to not positive defined, the function
1679 raises LinalgError.
1681 Both matrices ``A`` and ``B`` are stored in LAPACK banded storage.
1683 References
1684 ----------
1685 .. [1] M. F. Hutchinson and F. R. de Hoog, "Smoothing noisy data with
1686 spline functions," Numerische Mathematik, vol. 47, no. 1,
1687 pp. 99-106, 1985.
1688 :doi:`10.1007/BF01389878`
1690 """
1692 def find_b_inv_elem(i, j, U, D, B):
1693 rng = min(3, n - i - 1)
1694 rng_sum = 0.
1695 if j == 0:
1696 # use 2-nd formula from [1]
1697 for k in range(1, rng + 1):
1698 rng_sum -= U[-k - 1, i + k] * B[-k - 1, i + k]
1699 rng_sum += D[i]
1700 B[-1, i] = rng_sum
1701 else:
1702 # use 1-st formula from [1]
1703 for k in range(1, rng + 1):
1704 diag = abs(k - j)
1705 ind = i + min(k, j)
1706 rng_sum -= U[-k - 1, i + k] * B[-diag - 1, ind + diag]
1707 B[-j - 1, i + j] = rng_sum
1709 U = cholesky_banded(A)
1710 for i in range(2, 5):
1711 U[-i, i-1:] /= U[-1, :-i+1]
1712 D = 1. / (U[-1])**2
1713 U[-1] /= U[-1]
1715 n = U.shape[1]
1717 B = np.zeros(shape=(4, n))
1718 for i in range(n - 1, -1, -1):
1719 for j in range(min(3, n - i - 1), -1, -1):
1720 find_b_inv_elem(i, j, U, D, B)
1721 # the first row contains garbage and should be removed
1722 B[0] = [0.] * n
1723 return B
1725 def _gcv(lam, X, XtWX, wE, XtE):
1726 r"""
1727 Computes the generalized cross-validation criteria [1].
1729 Parameters
1730 ----------
1731 lam : float, (:math:`\lambda \geq 0`)
1732 Regularization parameter.
1733 X : array, shape (5, n)
1734 Matrix is stored in LAPACK banded storage.
1735 XtWX : array, shape (4, n)
1736 Product :math:`X^T W X` stored in LAPACK banded storage.
1737 wE : array, shape (5, n)
1738 Matrix :math:`W^{-1} E` stored in LAPACK banded storage.
1739 XtE : array, shape (4, n)
1740 Product :math:`X^T E` stored in LAPACK banded storage.
1742 Returns
1743 -------
1744 res : float
1745 Value of the GCV criteria with the regularization parameter
1746 :math:`\lambda`.
1748 Notes
1749 -----
1750 Criteria is computed from the formula (1.3.2) [3]:
1752 .. math:
1754 GCV(\lambda) = \dfrac{1}{n} \sum\limits_{k = 1}^{n} \dfrac{ \left(
1755 y_k - f_{\lambda}(x_k) \right)^2}{\left( 1 - \Tr{A}/n\right)^2}$.
1756 The criteria is discussed in section 1.3 [3].
1758 The numerator is computed using (2.2.4) [3] and the denominator is
1759 computed using an algorithm from [2] (see in the ``compute_b_inv``
1760 function).
1762 References
1763 ----------
1764 .. [1] G. Wahba, "Estimating the smoothing parameter" in Spline models
1765 for observational data, Philadelphia, Pennsylvania: Society for
1766 Industrial and Applied Mathematics, 1990, pp. 45-65.
1767 :doi:`10.1137/1.9781611970128`
1768 .. [2] M. F. Hutchinson and F. R. de Hoog, "Smoothing noisy data with
1769 spline functions," Numerische Mathematik, vol. 47, no. 1,
1770 pp. 99-106, 1985.
1771 :doi:`10.1007/BF01389878`
1772 .. [3] E. Zemlyanoy, "Generalized cross-validation smoothing splines",
1773 BSc thesis, 2022. Might be available (in Russian)
1774 `here <https://www.hse.ru/ba/am/students/diplomas/620910604>`_
1776 """
1777 # Compute the numerator from (2.2.4) [3]
1778 n = X.shape[1]
1779 c = solve_banded((2, 2), X + lam * wE, y)
1780 res = np.zeros(n)
1781 # compute ``W^{-1} E c`` with respect to banded-storage of ``E``
1782 tmp = wE * c
1783 for i in range(n):
1784 for j in range(max(0, i - n + 3), min(5, i + 3)):
1785 res[i] += tmp[j, i + 2 - j]
1786 numer = np.linalg.norm(lam * res)**2 / n
1788 # compute the denominator
1789 lhs = XtWX + lam * XtE
1790 try:
1791 b_banded = compute_b_inv(lhs)
1792 # compute the trace of the product b_banded @ XtX
1793 tr = b_banded * XtWX
1794 tr[:-1] *= 2
1795 # find the denominator
1796 denom = (1 - sum(sum(tr)) / n)**2
1797 except LinAlgError:
1798 # cholesky decomposition cannot be performed
1799 raise ValueError('Seems like the problem is ill-posed')
1801 res = numer / denom
1803 return res
1805 n = X.shape[1]
1807 XtWX = compute_banded_symmetric_XT_W_Y(X, w, X)
1808 XtE = compute_banded_symmetric_XT_W_Y(X, w, wE)
1810 def fun(lam):
1811 return _gcv(lam, X, XtWX, wE, XtE)
1813 gcv_est = minimize_scalar(fun, bounds=(0, n), method='Bounded')
1814 if gcv_est.success:
1815 return gcv_est.x
1816 raise ValueError(f"Unable to find minimum of the GCV "
1817 f"function: {gcv_est.message}")
1820def _coeff_of_divided_diff(x):
1821 """
1822 Returns the coefficients of the divided difference.
1824 Parameters
1825 ----------
1826 x : array, shape (n,)
1827 Array which is used for the computation of divided difference.
1829 Returns
1830 -------
1831 res : array_like, shape (n,)
1832 Coefficients of the divided difference.
1834 Notes
1835 -----
1836 Vector ``x`` should have unique elements, otherwise an error division by
1837 zero might be raised.
1839 No checks are performed.
1841 """
1842 n = x.shape[0]
1843 res = np.zeros(n)
1844 for i in range(n):
1845 pp = 1.
1846 for k in range(n):
1847 if k != i:
1848 pp *= (x[i] - x[k])
1849 res[i] = 1. / pp
1850 return res
1853def make_smoothing_spline(x, y, w=None, lam=None):
1854 r"""
1855 Compute the (coefficients of) smoothing cubic spline function using
1856 ``lam`` to control the tradeoff between the amount of smoothness of the
1857 curve and its proximity to the data. In case ``lam`` is None, using the
1858 GCV criteria [1] to find it.
1860 A smoothing spline is found as a solution to the regularized weighted
1861 linear regression problem:
1863 .. math::
1865 \sum\limits_{i=1}^n w_i\lvert y_i - f(x_i) \rvert^2 +
1866 \lambda\int\limits_{x_1}^{x_n} (f^{(2)}(u))^2 d u
1868 where :math:`f` is a spline function, :math:`w` is a vector of weights and
1869 :math:`\lambda` is a regularization parameter.
1871 If ``lam`` is None, we use the GCV criteria to find an optimal
1872 regularization parameter, otherwise we solve the regularized weighted
1873 linear regression problem with given parameter. The parameter controls
1874 the tradeoff in the following way: the larger the parameter becomes, the
1875 smoother the function gets.
1877 Parameters
1878 ----------
1879 x : array_like, shape (n,)
1880 Abscissas. `n` must be larger than 5.
1881 y : array_like, shape (n,)
1882 Ordinates. `n` must be larger than 5.
1883 w : array_like, shape (n,), optional
1884 Vector of weights. Default is ``np.ones_like(x)``.
1885 lam : float, (:math:`\lambda \geq 0`), optional
1886 Regularization parameter. If ``lam`` is None, then it is found from
1887 the GCV criteria. Default is None.
1889 Returns
1890 -------
1891 func : a BSpline object.
1892 A callable representing a spline in the B-spline basis
1893 as a solution of the problem of smoothing splines using
1894 the GCV criteria [1] in case ``lam`` is None, otherwise using the
1895 given parameter ``lam``.
1897 Notes
1898 -----
1899 This algorithm is a clean room reimplementation of the algorithm
1900 introduced by Woltring in FORTRAN [2]. The original version cannot be used
1901 in SciPy source code because of the license issues. The details of the
1902 reimplementation are discussed here (available only in Russian) [4].
1904 If the vector of weights ``w`` is None, we assume that all the points are
1905 equal in terms of weights, and vector of weights is vector of ones.
1907 Note that in weighted residual sum of squares, weights are not squared:
1908 :math:`\sum\limits_{i=1}^n w_i\lvert y_i - f(x_i) \rvert^2` while in
1909 ``splrep`` the sum is built from the squared weights.
1911 In cases when the initial problem is ill-posed (for example, the product
1912 :math:`X^T W X` where :math:`X` is a design matrix is not a positive
1913 defined matrix) a ValueError is raised.
1915 References
1916 ----------
1917 .. [1] G. Wahba, "Estimating the smoothing parameter" in Spline models for
1918 observational data, Philadelphia, Pennsylvania: Society for Industrial
1919 and Applied Mathematics, 1990, pp. 45-65.
1920 :doi:`10.1137/1.9781611970128`
1921 .. [2] H. J. Woltring, A Fortran package for generalized, cross-validatory
1922 spline smoothing and differentiation, Advances in Engineering
1923 Software, vol. 8, no. 2, pp. 104-113, 1986.
1924 :doi:`10.1016/0141-1195(86)90098-7`
1925 .. [3] T. Hastie, J. Friedman, and R. Tisbshirani, "Smoothing Splines" in
1926 The elements of Statistical Learning: Data Mining, Inference, and
1927 prediction, New York: Springer, 2017, pp. 241-249.
1928 :doi:`10.1007/978-0-387-84858-7`
1929 .. [4] E. Zemlyanoy, "Generalized cross-validation smoothing splines",
1930 BSc thesis, 2022.
1931 `<https://www.hse.ru/ba/am/students/diplomas/620910604>`_ (in
1932 Russian)
1934 Examples
1935 --------
1936 Generate some noisy data
1938 >>> import numpy as np
1939 >>> np.random.seed(1234)
1940 >>> n = 200
1941 >>> def func(x):
1942 ... return x**3 + x**2 * np.sin(4 * x)
1943 >>> x = np.sort(np.random.random_sample(n) * 4 - 2)
1944 >>> y = func(x) + np.random.normal(scale=1.5, size=n)
1946 Make a smoothing spline function
1948 >>> from scipy.interpolate import make_smoothing_spline
1949 >>> spl = make_smoothing_spline(x, y)
1951 Plot both
1953 >>> import matplotlib.pyplot as plt
1954 >>> grid = np.linspace(x[0], x[-1], 400)
1955 >>> plt.plot(grid, spl(grid), label='Spline')
1956 >>> plt.plot(grid, func(grid), label='Original function')
1957 >>> plt.scatter(x, y, marker='.')
1958 >>> plt.legend(loc='best')
1959 >>> plt.show()
1961 """
1963 x = np.ascontiguousarray(x, dtype=float)
1964 y = np.ascontiguousarray(y, dtype=float)
1966 if any(x[1:] - x[:-1] <= 0):
1967 raise ValueError('``x`` should be an ascending array')
1969 if x.ndim != 1 or y.ndim != 1 or x.shape[0] != y.shape[0]:
1970 raise ValueError('``x`` and ``y`` should be one dimensional and the'
1971 ' same size')
1973 if w is None:
1974 w = np.ones(len(x))
1975 else:
1976 w = np.ascontiguousarray(w)
1977 if any(w <= 0):
1978 raise ValueError('Invalid vector of weights')
1980 t = np.r_[[x[0]] * 3, x, [x[-1]] * 3]
1981 n = x.shape[0]
1983 if n <= 4:
1984 raise ValueError('``x`` and ``y`` length must be larger than 5')
1986 # It is known that the solution to the stated minimization problem exists
1987 # and is a natural cubic spline with vector of knots equal to the unique
1988 # elements of ``x`` [3], so we will solve the problem in the basis of
1989 # natural splines.
1991 # create design matrix in the B-spline basis
1992 X_bspl = BSpline.design_matrix(x, t, 3)
1993 # move from B-spline basis to the basis of natural splines using equations
1994 # (2.1.7) [4]
1995 # central elements
1996 X = np.zeros((5, n))
1997 for i in range(1, 4):
1998 X[i, 2: -2] = X_bspl[i: i - 4, 3: -3][np.diag_indices(n - 4)]
2000 # first elements
2001 X[1, 1] = X_bspl[0, 0]
2002 X[2, :2] = ((x[2] + x[1] - 2 * x[0]) * X_bspl[0, 0],
2003 X_bspl[1, 1] + X_bspl[1, 2])
2004 X[3, :2] = ((x[2] - x[0]) * X_bspl[1, 1], X_bspl[2, 2])
2006 # last elements
2007 X[1, -2:] = (X_bspl[-3, -3], (x[-1] - x[-3]) * X_bspl[-2, -2])
2008 X[2, -2:] = (X_bspl[-2, -3] + X_bspl[-2, -2],
2009 (2 * x[-1] - x[-2] - x[-3]) * X_bspl[-1, -1])
2010 X[3, -2] = X_bspl[-1, -1]
2012 # create penalty matrix and divide it by vector of weights: W^{-1} E
2013 wE = np.zeros((5, n))
2014 wE[2:, 0] = _coeff_of_divided_diff(x[:3]) / w[:3]
2015 wE[1:, 1] = _coeff_of_divided_diff(x[:4]) / w[:4]
2016 for j in range(2, n - 2):
2017 wE[:, j] = (x[j+2] - x[j-2]) * _coeff_of_divided_diff(x[j-2:j+3])\
2018 / w[j-2: j+3]
2020 wE[:-1, -2] = -_coeff_of_divided_diff(x[-4:]) / w[-4:]
2021 wE[:-2, -1] = _coeff_of_divided_diff(x[-3:]) / w[-3:]
2022 wE *= 6
2024 if lam is None:
2025 lam = _compute_optimal_gcv_parameter(X, wE, y, w)
2026 elif lam < 0.:
2027 raise ValueError('Regularization parameter should be non-negative')
2029 # solve the initial problem in the basis of natural splines
2030 c = solve_banded((2, 2), X + lam * wE, y)
2031 # move back to B-spline basis using equations (2.2.10) [4]
2032 c_ = np.r_[c[0] * (t[5] + t[4] - 2 * t[3]) + c[1],
2033 c[0] * (t[5] - t[3]) + c[1],
2034 c[1: -1],
2035 c[-1] * (t[-4] - t[-6]) + c[-2],
2036 c[-1] * (2 * t[-4] - t[-5] - t[-6]) + c[-2]]
2038 return BSpline.construct_fast(t, c_, 3)