Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_bsplines.py: 8%
577 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import operator
3import numpy as np
4from numpy.core.multiarray import normalize_axis_index
5from scipy.linalg import (get_lapack_funcs, LinAlgError,
6 cholesky_banded, cho_solve_banded,
7 solve, solve_banded)
8from scipy.optimize import minimize_scalar
9from . import _bspl
10from . import _fitpack_impl
11from scipy._lib._util import prod
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)
879def _woodbury_algorithm(A, ur, ll, b, k):
880 '''
881 Solve a cyclic banded linear system with upper right
882 and lower blocks of size ``(k-1) / 2`` using
883 the Woodbury formula
885 Parameters
886 ----------
887 A : 2-D array, shape(k, n)
888 Matrix of diagonals of original matrix(see
889 ``solve_banded`` documentation).
890 ur : 2-D array, shape(bs, bs)
891 Upper right block matrix.
892 ll : 2-D array, shape(bs, bs)
893 Lower left block matrix.
894 b : 1-D array, shape(n,)
895 Vector of constant terms of the system of linear equations.
896 k : int
897 B-spline degree.
899 Returns
900 -------
901 c : 1-D array, shape(n,)
902 Solution of the original system of linear equations.
904 Notes
905 -----
906 This algorithm works only for systems with banded matrix A plus
907 a correction term U @ V.T, where the matrix U @ V.T gives upper right
908 and lower left block of A
909 The system is solved with the following steps:
910 1. New systems of linear equations are constructed:
911 A @ z_i = u_i,
912 u_i - columnn vector of U,
913 i = 1, ..., k - 1
914 2. Matrix Z is formed from vectors z_i:
915 Z = [ z_1 | z_2 | ... | z_{k - 1} ]
916 3. Matrix H = (1 + V.T @ Z)^{-1}
917 4. The system A' @ y = b is solved
918 5. x = y - Z @ (H @ V.T @ y)
919 Also, ``n`` should be greater than ``k``, otherwise corner block
920 elements will intersect with diagonals.
922 Examples
923 --------
924 Consider the case of n = 8, k = 5 (size of blocks - 2 x 2).
925 The matrix of a system: U: V:
926 x x x * * a b a b 0 0 0 0 1 0
927 x x x x * * c 0 c 0 0 0 0 0 1
928 x x x x x * * 0 0 0 0 0 0 0 0
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 d * * x x x x 0 0 d 0 1 0 0 0
932 e f * * x x x 0 0 e f 0 1 0 0
934 References
935 ----------
936 .. [1] William H. Press, Saul A. Teukolsky, William T. Vetterling
937 and Brian P. Flannery, Numerical Recipes, 2007, Section 2.7.3
939 '''
940 k_mod = k - k % 2
941 bs = int((k - 1) / 2) + (k + 1) % 2
943 n = A.shape[1] + 1
944 U = np.zeros((n - 1, k_mod))
945 VT = np.zeros((k_mod, n - 1)) # V transpose
947 # upper right block
948 U[:bs, :bs] = ur
949 VT[np.arange(bs), np.arange(bs) - bs] = 1
951 # lower left block
952 U[-bs:, -bs:] = ll
953 VT[np.arange(bs) - bs, np.arange(bs)] = 1
955 Z = solve_banded((bs, bs), A, U)
957 H = solve(np.identity(k_mod) + VT @ Z, np.identity(k_mod))
959 y = solve_banded((bs, bs), A, b)
960 c = y - Z @ (H @ (VT @ y))
962 return c
964def _periodic_knots(x, k):
965 '''
966 returns vector of nodes on circle
967 '''
968 xc = np.copy(x)
969 n = len(xc)
970 if k % 2 == 0:
971 dx = np.diff(xc)
972 xc[1: -1] -= dx[:-1] / 2
973 dx = np.diff(xc)
974 t = np.zeros(n + 2 * k)
975 t[k: -k] = xc
976 for i in range(0, k):
977 # filling first `k` elements in descending order
978 t[k - i - 1] = t[k - i] - dx[-(i % (n - 1)) - 1]
979 # filling last `k` elements in ascending order
980 t[-k + i] = t[-k + i - 1] + dx[i % (n - 1)]
981 return t
984def _make_interp_per_full_matr(x, y, t, k):
985 '''
986 Returns a solution of a system for B-spline interpolation with periodic
987 boundary conditions. First ``k - 1`` rows of matrix are condtions of
988 periodicity (continuity of ``k - 1`` derivatives at the boundary points).
989 Last ``n`` rows are interpolation conditions.
990 RHS is ``k - 1`` zeros and ``n`` ordinates in this case.
992 Parameters
993 ----------
994 x : 1-D array, shape (n,)
995 Values of x - coordinate of a given set of points.
996 y : 1-D array, shape (n,)
997 Values of y - coordinate of a given set of points.
998 t : 1-D array, shape(n+2*k,)
999 Vector of knots.
1000 k : int
1001 The maximum degree of spline
1003 Returns
1004 -------
1005 c : 1-D array, shape (n+k-1,)
1006 B-spline coefficients
1008 Notes
1009 -----
1010 ``t`` is supposed to be taken on circle.
1012 '''
1014 x, y, t = map(np.asarray, (x, y, t))
1016 n = x.size
1017 # LHS: the collocation matrix + derivatives at edges
1018 matr = np.zeros((n + k - 1, n + k - 1))
1020 # derivatives at x[0] and x[-1]:
1021 for i in range(k - 1):
1022 bb = _bspl.evaluate_all_bspl(t, k, x[0], k, nu=i + 1)
1023 matr[i, : k + 1] += bb
1024 bb = _bspl.evaluate_all_bspl(t, k, x[-1], n + k - 1, nu=i + 1)[:-1]
1025 matr[i, -k:] -= bb
1027 # collocation matrix
1028 for i in range(n):
1029 xval = x[i]
1030 # find interval
1031 if xval == t[k]:
1032 left = k
1033 else:
1034 left = np.searchsorted(t, xval) - 1
1036 # fill a row
1037 bb = _bspl.evaluate_all_bspl(t, k, xval, left)
1038 matr[i + k - 1, left-k:left+1] = bb
1040 # RHS
1041 b = np.r_[[0] * (k - 1), y]
1043 c = solve(matr, b)
1044 return c
1046def _make_periodic_spline(x, y, t, k, axis):
1047 '''
1048 Compute the (coefficients of) interpolating B-spline with periodic
1049 boundary conditions.
1051 Parameters
1052 ----------
1053 x : array_like, shape (n,)
1054 Abscissas.
1055 y : array_like, shape (n,)
1056 Ordinates.
1057 k : int
1058 B-spline degree.
1059 t : array_like, shape (n + 2 * k,).
1060 Knots taken on a circle, ``k`` on the left and ``k`` on the right
1061 of the vector ``x``.
1063 Returns
1064 -------
1065 b : a BSpline object of the degree ``k`` and with knots ``t``.
1067 Notes
1068 -----
1069 The original system is formed by ``n + k - 1`` equations where the first
1070 ``k - 1`` of them stand for the ``k - 1`` derivatives continuity on the
1071 edges while the other equations correspond to an interpolating case
1072 (matching all the input points). Due to a special form of knot vector, it
1073 can be proved that in the original system the first and last ``k``
1074 coefficients of a spline function are the same, respectively. It follows
1075 from the fact that all ``k - 1`` derivatives are equal term by term at ends
1076 and that the matrix of the original system of linear equations is
1077 non-degenerate. So, we can reduce the number of equations to ``n - 1``
1078 (first ``k - 1`` equations could be reduced). Another trick of this
1079 implementation is cyclic shift of values of B-splines due to equality of
1080 ``k`` unknown coefficients. With this we can receive matrix of the system
1081 with upper right and lower left blocks, and ``k`` diagonals. It allows
1082 to use Woodbury formula to optimize the computations.
1084 '''
1085 n = y.shape[0]
1087 extradim = prod(y.shape[1:])
1088 y_new = y.reshape(n, extradim)
1089 c = np.zeros((n + k - 1, extradim))
1091 # n <= k case is solved with full matrix
1092 if n <= k:
1093 for i in range(extradim):
1094 c[:, i] = _make_interp_per_full_matr(x, y_new[:, i], t, k)
1095 c = np.ascontiguousarray(c.reshape((n + k - 1,) + y.shape[1:]))
1096 return BSpline.construct_fast(t, c, k, extrapolate='periodic', axis=axis)
1098 nt = len(t) - k - 1
1100 # size of block elements
1101 kul = int(k / 2)
1103 # kl = ku = k
1104 ab = np.zeros((3 * k + 1, nt), dtype=np.float_, order='F')
1106 # upper right and lower left blocks
1107 ur = np.zeros((kul, kul))
1108 ll = np.zeros_like(ur)
1110 # `offset` is made to shift all the non-zero elements to the end of the
1111 # matrix
1112 _bspl._colloc(x, t, k, ab, offset=k)
1114 # remove zeros before the matrix
1115 ab = ab[-k - (k + 1) % 2:, :]
1117 # The least elements in rows (except repetitions) are diagonals
1118 # of block matrices. Upper right matrix is an upper triangular
1119 # matrix while lower left is a lower triangular one.
1120 for i in range(kul):
1121 ur += np.diag(ab[-i - 1, i: kul], k=i)
1122 ll += np.diag(ab[i, -kul - (k % 2): n - 1 + 2 * kul - i], k=-i)
1124 # remove elements that occur in the last point
1125 # (first and last points are equivalent)
1126 A = ab[:, kul: -k + kul]
1128 for i in range(extradim):
1129 cc = _woodbury_algorithm(A, ur, ll, y_new[:, i][:-1], k)
1130 c[:, i] = np.concatenate((cc[-kul:], cc, cc[:kul + k % 2]))
1131 c = np.ascontiguousarray(c.reshape((n + k - 1,) + y.shape[1:]))
1132 return BSpline.construct_fast(t, c, k, extrapolate='periodic', axis=axis)
1134def make_interp_spline(x, y, k=3, t=None, bc_type=None, axis=0,
1135 check_finite=True):
1136 """Compute the (coefficients of) interpolating B-spline.
1138 Parameters
1139 ----------
1140 x : array_like, shape (n,)
1141 Abscissas.
1142 y : array_like, shape (n, ...)
1143 Ordinates.
1144 k : int, optional
1145 B-spline degree. Default is cubic, ``k = 3``.
1146 t : array_like, shape (nt + k + 1,), optional.
1147 Knots.
1148 The number of knots needs to agree with the number of data points and
1149 the number of derivatives at the edges. Specifically, ``nt - n`` must
1150 equal ``len(deriv_l) + len(deriv_r)``.
1151 bc_type : 2-tuple or None
1152 Boundary conditions.
1153 Default is None, which means choosing the boundary conditions
1154 automatically. Otherwise, it must be a length-two tuple where the first
1155 element (``deriv_l``) sets the boundary conditions at ``x[0]`` and
1156 the second element (``deriv_r``) sets the boundary conditions at
1157 ``x[-1]``. Each of these must be an iterable of pairs
1158 ``(order, value)`` which gives the values of derivatives of specified
1159 orders at the given edge of the interpolation interval.
1160 Alternatively, the following string aliases are recognized:
1162 * ``"clamped"``: The first derivatives at the ends are zero. This is
1163 equivalent to ``bc_type=([(1, 0.0)], [(1, 0.0)])``.
1164 * ``"natural"``: The second derivatives at ends are zero. This is
1165 equivalent to ``bc_type=([(2, 0.0)], [(2, 0.0)])``.
1166 * ``"not-a-knot"`` (default): The first and second segments are the
1167 same polynomial. This is equivalent to having ``bc_type=None``.
1168 * ``"periodic"``: The values and the first ``k-1`` derivatives at the
1169 ends are equivalent.
1171 axis : int, optional
1172 Interpolation axis. Default is 0.
1173 check_finite : bool, optional
1174 Whether to check that the input arrays contain only finite numbers.
1175 Disabling may give a performance gain, but may result in problems
1176 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1177 Default is True.
1179 Returns
1180 -------
1181 b : a BSpline object of the degree ``k`` and with knots ``t``.
1183 Examples
1184 --------
1186 Use cubic interpolation on Chebyshev nodes:
1188 >>> import numpy as np
1189 >>> import matplotlib.pyplot as plt
1190 >>> def cheb_nodes(N):
1191 ... jj = 2.*np.arange(N) + 1
1192 ... x = np.cos(np.pi * jj / 2 / N)[::-1]
1193 ... return x
1195 >>> x = cheb_nodes(20)
1196 >>> y = np.sqrt(1 - x**2)
1198 >>> from scipy.interpolate import BSpline, make_interp_spline
1199 >>> b = make_interp_spline(x, y)
1200 >>> np.allclose(b(x), y)
1201 True
1203 Note that the default is a cubic spline with a not-a-knot boundary condition
1205 >>> b.k
1206 3
1208 Here we use a 'natural' spline, with zero 2nd derivatives at edges:
1210 >>> l, r = [(2, 0.0)], [(2, 0.0)]
1211 >>> b_n = make_interp_spline(x, y, bc_type=(l, r)) # or, bc_type="natural"
1212 >>> np.allclose(b_n(x), y)
1213 True
1214 >>> x0, x1 = x[0], x[-1]
1215 >>> np.allclose([b_n(x0, 2), b_n(x1, 2)], [0, 0])
1216 True
1218 Interpolation of parametric curves is also supported. As an example, we
1219 compute a discretization of a snail curve in polar coordinates
1221 >>> phi = np.linspace(0, 2.*np.pi, 40)
1222 >>> r = 0.3 + np.cos(phi)
1223 >>> x, y = r*np.cos(phi), r*np.sin(phi) # convert to Cartesian coordinates
1225 Build an interpolating curve, parameterizing it by the angle
1227 >>> spl = make_interp_spline(phi, np.c_[x, y])
1229 Evaluate the interpolant on a finer grid (note that we transpose the result
1230 to unpack it into a pair of x- and y-arrays)
1232 >>> phi_new = np.linspace(0, 2.*np.pi, 100)
1233 >>> x_new, y_new = spl(phi_new).T
1235 Plot the result
1237 >>> plt.plot(x, y, 'o')
1238 >>> plt.plot(x_new, y_new, '-')
1239 >>> plt.show()
1241 Build a B-spline curve with 2 dimensional y
1243 >>> x = np.linspace(0, 2*np.pi, 10)
1244 >>> y = np.array([np.sin(x), np.cos(x)])
1246 Periodic condition is satisfied because y coordinates of points on the ends
1247 are equivalent
1249 >>> ax = plt.axes(projection='3d')
1250 >>> xx = np.linspace(0, 2*np.pi, 100)
1251 >>> bspl = make_interp_spline(x, y, k=5, bc_type='periodic', axis=1)
1252 >>> ax.plot3D(xx, *bspl(xx))
1253 >>> ax.scatter3D(x, *y, color='red')
1254 >>> plt.show()
1256 See Also
1257 --------
1258 BSpline : base class representing the B-spline objects
1259 CubicSpline : a cubic spline in the polynomial basis
1260 make_lsq_spline : a similar factory function for spline fitting
1261 UnivariateSpline : a wrapper over FITPACK spline fitting routines
1262 splrep : a wrapper over FITPACK spline fitting routines
1264 """
1265 # convert string aliases for the boundary conditions
1266 if bc_type is None or bc_type == 'not-a-knot' or bc_type == 'periodic':
1267 deriv_l, deriv_r = None, None
1268 elif isinstance(bc_type, str):
1269 deriv_l, deriv_r = bc_type, bc_type
1270 else:
1271 try:
1272 deriv_l, deriv_r = bc_type
1273 except TypeError as e:
1274 raise ValueError("Unknown boundary condition: %s" % bc_type) from e
1276 y = np.asarray(y)
1278 axis = normalize_axis_index(axis, y.ndim)
1280 x = _as_float_array(x, check_finite)
1281 y = _as_float_array(y, check_finite)
1283 y = np.moveaxis(y, axis, 0) # now internally interp axis is zero
1285 # sanity check the input
1286 if bc_type == 'periodic' and not np.allclose(y[0], y[-1], atol=1e-15):
1287 raise ValueError("First and last points does not match while "
1288 "periodic case expected")
1289 if x.size != y.shape[0]:
1290 raise ValueError('Shapes of x {} and y {} are incompatible'
1291 .format(x.shape, y.shape))
1292 if np.any(x[1:] == x[:-1]):
1293 raise ValueError("Expect x to not have duplicates")
1294 if x.ndim != 1 or np.any(x[1:] < x[:-1]):
1295 raise ValueError("Expect x to be a 1D strictly increasing sequence.")
1297 # special-case k=0 right away
1298 if k == 0:
1299 if any(_ is not None for _ in (t, deriv_l, deriv_r)):
1300 raise ValueError("Too much info for k=0: t and bc_type can only "
1301 "be None.")
1302 t = np.r_[x, x[-1]]
1303 c = np.asarray(y)
1304 c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
1305 return BSpline.construct_fast(t, c, k, axis=axis)
1307 # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
1308 if k == 1 and t is None:
1309 if not (deriv_l is None and deriv_r is None):
1310 raise ValueError("Too much info for k=1: bc_type can only be None.")
1311 t = np.r_[x[0], x, x[-1]]
1312 c = np.asarray(y)
1313 c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
1314 return BSpline.construct_fast(t, c, k, axis=axis)
1316 k = operator.index(k)
1318 if bc_type == 'periodic' and t is not None:
1319 raise NotImplementedError("For periodic case t is constructed "
1320 "automatically and can not be passed manually")
1322 # come up with a sensible knot vector, if needed
1323 if t is None:
1324 if deriv_l is None and deriv_r is None:
1325 if bc_type == 'periodic':
1326 t = _periodic_knots(x, k)
1327 elif k == 2:
1328 # OK, it's a bit ad hoc: Greville sites + omit
1329 # 2nd and 2nd-to-last points, a la not-a-knot
1330 t = (x[1:] + x[:-1]) / 2.
1331 t = np.r_[(x[0],)*(k+1),
1332 t[1:-1],
1333 (x[-1],)*(k+1)]
1334 else:
1335 t = _not_a_knot(x, k)
1336 else:
1337 t = _augknt(x, k)
1339 t = _as_float_array(t, check_finite)
1341 if k < 0:
1342 raise ValueError("Expect non-negative k.")
1343 if t.ndim != 1 or np.any(t[1:] < t[:-1]):
1344 raise ValueError("Expect t to be a 1-D sorted array_like.")
1345 if t.size < x.size + k + 1:
1346 raise ValueError('Got %d knots, need at least %d.' %
1347 (t.size, x.size + k + 1))
1348 if (x[0] < t[k]) or (x[-1] > t[-k]):
1349 raise ValueError('Out of bounds w/ x = %s.' % x)
1351 if bc_type == 'periodic':
1352 return _make_periodic_spline(x, y, t, k, axis)
1354 # Here : deriv_l, r = [(nu, value), ...]
1355 deriv_l = _convert_string_aliases(deriv_l, y.shape[1:])
1356 deriv_l_ords, deriv_l_vals = _process_deriv_spec(deriv_l)
1357 nleft = deriv_l_ords.shape[0]
1359 deriv_r = _convert_string_aliases(deriv_r, y.shape[1:])
1360 deriv_r_ords, deriv_r_vals = _process_deriv_spec(deriv_r)
1361 nright = deriv_r_ords.shape[0]
1363 # have `n` conditions for `nt` coefficients; need nt-n derivatives
1364 n = x.size
1365 nt = t.size - k - 1
1367 if nt - n != nleft + nright:
1368 raise ValueError("The number of derivatives at boundaries does not "
1369 "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
1371 # bail out if the `y` array is zero-sized
1372 if y.size == 0:
1373 c = np.zeros((nt,) + y.shape[1:], dtype=float)
1374 return BSpline.construct_fast(t, c, k, axis=axis)
1376 # set up the LHS: the collocation matrix + derivatives at boundaries
1377 kl = ku = k
1378 ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
1379 _bspl._colloc(x, t, k, ab, offset=nleft)
1380 if nleft > 0:
1381 _bspl._handle_lhs_derivatives(t, k, x[0], ab, kl, ku, deriv_l_ords)
1382 if nright > 0:
1383 _bspl._handle_lhs_derivatives(t, k, x[-1], ab, kl, ku, deriv_r_ords,
1384 offset=nt-nright)
1386 # set up the RHS: values to interpolate (+ derivative values, if any)
1387 extradim = prod(y.shape[1:])
1388 rhs = np.empty((nt, extradim), dtype=y.dtype)
1389 if nleft > 0:
1390 rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
1391 rhs[nleft:nt - nright] = y.reshape(-1, extradim)
1392 if nright > 0:
1393 rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
1395 # solve Ab @ x = rhs; this is the relevant part of linalg.solve_banded
1396 if check_finite:
1397 ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
1398 gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
1399 lu, piv, c, info = gbsv(kl, ku, ab, rhs,
1400 overwrite_ab=True, overwrite_b=True)
1402 if info > 0:
1403 raise LinAlgError("Collocation matrix is singular.")
1404 elif info < 0:
1405 raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
1407 c = np.ascontiguousarray(c.reshape((nt,) + y.shape[1:]))
1408 return BSpline.construct_fast(t, c, k, axis=axis)
1411def make_lsq_spline(x, y, t, k=3, w=None, axis=0, check_finite=True):
1412 r"""Compute the (coefficients of) an LSQ (Least SQuared) based
1413 fitting B-spline.
1415 The result is a linear combination
1417 .. math::
1419 S(x) = \sum_j c_j B_j(x; t)
1421 of the B-spline basis elements, :math:`B_j(x; t)`, which minimizes
1423 .. math::
1425 \sum_{j} \left( w_j \times (S(x_j) - y_j) \right)^2
1427 Parameters
1428 ----------
1429 x : array_like, shape (m,)
1430 Abscissas.
1431 y : array_like, shape (m, ...)
1432 Ordinates.
1433 t : array_like, shape (n + k + 1,).
1434 Knots.
1435 Knots and data points must satisfy Schoenberg-Whitney conditions.
1436 k : int, optional
1437 B-spline degree. Default is cubic, ``k = 3``.
1438 w : array_like, shape (m,), optional
1439 Weights for spline fitting. Must be positive. If ``None``,
1440 then weights are all equal.
1441 Default is ``None``.
1442 axis : int, optional
1443 Interpolation axis. Default is zero.
1444 check_finite : bool, optional
1445 Whether to check that the input arrays contain only finite numbers.
1446 Disabling may give a performance gain, but may result in problems
1447 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1448 Default is True.
1450 Returns
1451 -------
1452 b : a BSpline object of the degree ``k`` with knots ``t``.
1454 Notes
1455 -----
1456 The number of data points must be larger than the spline degree ``k``.
1458 Knots ``t`` must satisfy the Schoenberg-Whitney conditions,
1459 i.e., there must be a subset of data points ``x[j]`` such that
1460 ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``.
1462 Examples
1463 --------
1464 Generate some noisy data:
1466 >>> import numpy as np
1467 >>> import matplotlib.pyplot as plt
1468 >>> rng = np.random.default_rng()
1469 >>> x = np.linspace(-3, 3, 50)
1470 >>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(50)
1472 Now fit a smoothing cubic spline with a pre-defined internal knots.
1473 Here we make the knot vector (k+1)-regular by adding boundary knots:
1475 >>> from scipy.interpolate import make_lsq_spline, BSpline
1476 >>> t = [-1, 0, 1]
1477 >>> k = 3
1478 >>> t = np.r_[(x[0],)*(k+1),
1479 ... t,
1480 ... (x[-1],)*(k+1)]
1481 >>> spl = make_lsq_spline(x, y, t, k)
1483 For comparison, we also construct an interpolating spline for the same
1484 set of data:
1486 >>> from scipy.interpolate import make_interp_spline
1487 >>> spl_i = make_interp_spline(x, y)
1489 Plot both:
1491 >>> xs = np.linspace(-3, 3, 100)
1492 >>> plt.plot(x, y, 'ro', ms=5)
1493 >>> plt.plot(xs, spl(xs), 'g-', lw=3, label='LSQ spline')
1494 >>> plt.plot(xs, spl_i(xs), 'b-', lw=3, alpha=0.7, label='interp spline')
1495 >>> plt.legend(loc='best')
1496 >>> plt.show()
1498 **NaN handling**: If the input arrays contain ``nan`` values, the result is
1499 not useful since the underlying spline fitting routines cannot deal with
1500 ``nan``. A workaround is to use zero weights for not-a-number data points:
1502 >>> y[8] = np.nan
1503 >>> w = np.isnan(y)
1504 >>> y[w] = 0.
1505 >>> tck = make_lsq_spline(x, y, t, w=~w)
1507 Notice the need to replace a ``nan`` by a numerical value (precise value
1508 does not matter as long as the corresponding weight is zero.)
1510 See Also
1511 --------
1512 BSpline : base class representing the B-spline objects
1513 make_interp_spline : a similar factory function for interpolating splines
1514 LSQUnivariateSpline : a FITPACK-based spline fitting routine
1515 splrep : a FITPACK-based fitting routine
1517 """
1518 x = _as_float_array(x, check_finite)
1519 y = _as_float_array(y, check_finite)
1520 t = _as_float_array(t, check_finite)
1521 if w is not None:
1522 w = _as_float_array(w, check_finite)
1523 else:
1524 w = np.ones_like(x)
1525 k = operator.index(k)
1527 axis = normalize_axis_index(axis, y.ndim)
1529 y = np.moveaxis(y, axis, 0) # now internally interp axis is zero
1531 if x.ndim != 1 or np.any(x[1:] - x[:-1] <= 0):
1532 raise ValueError("Expect x to be a 1-D sorted array_like.")
1533 if x.shape[0] < k+1:
1534 raise ValueError("Need more x points.")
1535 if k < 0:
1536 raise ValueError("Expect non-negative k.")
1537 if t.ndim != 1 or np.any(t[1:] - t[:-1] < 0):
1538 raise ValueError("Expect t to be a 1-D sorted array_like.")
1539 if x.size != y.shape[0]:
1540 raise ValueError('Shapes of x {} and y {} are incompatible'
1541 .format(x.shape, y.shape))
1542 if k > 0 and np.any((x < t[k]) | (x > t[-k])):
1543 raise ValueError('Out of bounds w/ x = %s.' % x)
1544 if x.size != w.size:
1545 raise ValueError('Shapes of x {} and w {} are incompatible'
1546 .format(x.shape, w.shape))
1548 # number of coefficients
1549 n = t.size - k - 1
1551 # construct A.T @ A and rhs with A the collocation matrix, and
1552 # rhs = A.T @ y for solving the LSQ problem ``A.T @ A @ c = A.T @ y``
1553 lower = True
1554 extradim = prod(y.shape[1:])
1555 ab = np.zeros((k+1, n), dtype=np.float_, order='F')
1556 rhs = np.zeros((n, extradim), dtype=y.dtype, order='F')
1557 _bspl._norm_eq_lsq(x, t, k,
1558 y.reshape(-1, extradim),
1559 w,
1560 ab, rhs)
1561 rhs = rhs.reshape((n,) + y.shape[1:])
1563 # have observation matrix & rhs, can solve the LSQ problem
1564 cho_decomp = cholesky_banded(ab, overwrite_ab=True, lower=lower,
1565 check_finite=check_finite)
1566 c = cho_solve_banded((cho_decomp, lower), rhs, overwrite_b=True,
1567 check_finite=check_finite)
1569 c = np.ascontiguousarray(c)
1570 return BSpline.construct_fast(t, c, k, axis=axis)
1573#############################
1574# Smoothing spline helpers #
1575#############################
1577def _compute_optimal_gcv_parameter(X, wE, y, w):
1578 """
1579 Returns an optimal regularization parameter from the GCV criteria [1].
1581 Parameters
1582 ----------
1583 X : array, shape (5, n)
1584 5 bands of the design matrix ``X`` stored in LAPACK banded storage.
1585 wE : array, shape (5, n)
1586 5 bands of the penalty matrix :math:`W^{-1} E` stored in LAPACK banded
1587 storage.
1588 y : array, shape (n,)
1589 Ordinates.
1590 w : array, shape (n,)
1591 Vector of weights.
1593 Returns
1594 -------
1595 lam : float
1596 An optimal from the GCV criteria point of view regularization
1597 parameter.
1599 Notes
1600 -----
1601 No checks are performed.
1603 References
1604 ----------
1605 .. [1] G. Wahba, "Estimating the smoothing parameter" in Spline models
1606 for observational data, Philadelphia, Pennsylvania: Society for
1607 Industrial and Applied Mathematics, 1990, pp. 45-65.
1608 :doi:`10.1137/1.9781611970128`
1610 """
1612 def compute_banded_symmetric_XT_W_Y(X, w, Y):
1613 """
1614 Assuming that the product :math:`X^T W Y` is symmetric and both ``X``
1615 and ``Y`` are 5-banded, compute the unique bands of the product.
1617 Parameters
1618 ----------
1619 X : array, shape (5, n)
1620 5 bands of the matrix ``X`` stored in LAPACK banded storage.
1621 w : array, shape (n,)
1622 Array of weights
1623 Y : array, shape (5, n)
1624 5 bands of the matrix ``Y`` stored in LAPACK banded storage.
1626 Returns
1627 -------
1628 res : array, shape (4, n)
1629 The result of the product :math:`X^T Y` stored in the banded way.
1631 Notes
1632 -----
1633 As far as the matrices ``X`` and ``Y`` are 5-banded, their product
1634 :math:`X^T W Y` is 7-banded. It is also symmetric, so we can store only
1635 unique diagonals.
1637 """
1638 # compute W Y
1639 W_Y = np.copy(Y)
1641 W_Y[2] *= w
1642 for i in range(2):
1643 W_Y[i, 2 - i:] *= w[:-2 + i]
1644 W_Y[3 + i, :-1 - i] *= w[1 + i:]
1646 n = X.shape[1]
1647 res = np.zeros((4, n))
1648 for i in range(n):
1649 for j in range(min(n-i, 4)):
1650 res[-j-1, i + j] = sum(X[j:, i] * W_Y[:5-j, i + j])
1651 return res
1653 def compute_b_inv(A):
1654 """
1655 Inverse 3 central bands of matrix :math:`A=U^T D^{-1} U` assuming that
1656 ``U`` is a unit upper triangular banded matrix using an algorithm
1657 proposed in [1].
1659 Parameters
1660 ----------
1661 A : array, shape (4, n)
1662 Matrix to inverse, stored in LAPACK banded storage.
1664 Returns
1665 -------
1666 B : array, shape (4, n)
1667 3 unique bands of the symmetric matrix that is an inverse to ``A``.
1668 The first row is filled with zeros.
1670 Notes
1671 -----
1672 The algorithm is based on the cholesky decomposition and, therefore,
1673 in case matrix ``A`` is close to not positive defined, the function
1674 raises LinalgError.
1676 Both matrices ``A`` and ``B`` are stored in LAPACK banded storage.
1678 References
1679 ----------
1680 .. [1] M. F. Hutchinson and F. R. de Hoog, "Smoothing noisy data with
1681 spline functions," Numerische Mathematik, vol. 47, no. 1,
1682 pp. 99-106, 1985.
1683 :doi:`10.1007/BF01389878`
1685 """
1687 def find_b_inv_elem(i, j, U, D, B):
1688 rng = min(3, n - i - 1)
1689 rng_sum = 0.
1690 if j == 0:
1691 # use 2-nd formula from [1]
1692 for k in range(1, rng + 1):
1693 rng_sum -= U[-k - 1, i + k] * B[-k - 1, i + k]
1694 rng_sum += D[i]
1695 B[-1, i] = rng_sum
1696 else:
1697 # use 1-st formula from [1]
1698 for k in range(1, rng + 1):
1699 diag = abs(k - j)
1700 ind = i + min(k, j)
1701 rng_sum -= U[-k - 1, i + k] * B[-diag - 1, ind + diag]
1702 B[-j - 1, i + j] = rng_sum
1704 U = cholesky_banded(A)
1705 for i in range(2, 5):
1706 U[-i, i-1:] /= U[-1, :-i+1]
1707 D = 1. / (U[-1])**2
1708 U[-1] /= U[-1]
1710 n = U.shape[1]
1712 B = np.zeros(shape=(4, n))
1713 for i in range(n - 1, -1, -1):
1714 for j in range(min(3, n - i - 1), -1, -1):
1715 find_b_inv_elem(i, j, U, D, B)
1716 # the first row contains garbage and should be removed
1717 B[0] = [0.] * n
1718 return B
1720 def _gcv(lam, X, XtWX, wE, XtE):
1721 r"""
1722 Computes the generalized cross-validation criteria [1].
1724 Parameters
1725 ----------
1726 lam : float, (:math:`\lambda \geq 0`)
1727 Regularization parameter.
1728 X : array, shape (5, n)
1729 Matrix is stored in LAPACK banded storage.
1730 XtWX : array, shape (4, n)
1731 Product :math:`X^T W X` stored in LAPACK banded storage.
1732 wE : array, shape (5, n)
1733 Matrix :math:`W^{-1} E` stored in LAPACK banded storage.
1734 XtE : array, shape (4, n)
1735 Product :math:`X^T E` stored in LAPACK banded storage.
1737 Returns
1738 -------
1739 res : float
1740 Value of the GCV criteria with the regularization parameter
1741 :math:`\lambda`.
1743 Notes
1744 -----
1745 Criteria is computed from the formula (1.3.2) [3]:
1747 .. math:
1749 GCV(\lambda) = \dfrac{1}{n} \sum\limits_{k = 1}^{n} \dfrac{ \left(
1750 y_k - f_{\lambda}(x_k) \right)^2}{\left( 1 - \Tr{A}/n\right)^2}$.
1751 The criteria is discussed in section 1.3 [3].
1753 The numerator is computed using (2.2.4) [3] and the denominator is
1754 computed using an algorithm from [2] (see in the ``compute_b_inv``
1755 function).
1757 References
1758 ----------
1759 .. [1] G. Wahba, "Estimating the smoothing parameter" in Spline models
1760 for observational data, Philadelphia, Pennsylvania: Society for
1761 Industrial and Applied Mathematics, 1990, pp. 45-65.
1762 :doi:`10.1137/1.9781611970128`
1763 .. [2] M. F. Hutchinson and F. R. de Hoog, "Smoothing noisy data with
1764 spline functions," Numerische Mathematik, vol. 47, no. 1,
1765 pp. 99-106, 1985.
1766 :doi:`10.1007/BF01389878`
1767 .. [3] E. Zemlyanoy, "Generalized cross-validation smoothing splines",
1768 BSc thesis, 2022. Might be available (in Russian)
1769 `here <https://www.hse.ru/ba/am/students/diplomas/620910604>`_
1771 """
1772 # Compute the numerator from (2.2.4) [3]
1773 n = X.shape[1]
1774 c = solve_banded((2, 2), X + lam * wE, y)
1775 res = np.zeros(n)
1776 # compute ``W^{-1} E c`` with respect to banded-storage of ``E``
1777 tmp = wE * c
1778 for i in range(n):
1779 for j in range(max(0, i - n + 3), min(5, i + 3)):
1780 res[i] += tmp[j, i + 2 - j]
1781 numer = np.linalg.norm(lam * res)**2 / n
1783 # compute the denominator
1784 lhs = XtWX + lam * XtE
1785 try:
1786 b_banded = compute_b_inv(lhs)
1787 # compute the trace of the product b_banded @ XtX
1788 tr = b_banded * XtWX
1789 tr[:-1] *= 2
1790 # find the denominator
1791 denom = (1 - sum(sum(tr)) / n)**2
1792 except LinAlgError:
1793 # cholesky decomposition cannot be performed
1794 raise ValueError('Seems like the problem is ill-posed')
1796 res = numer / denom
1798 return res
1800 n = X.shape[1]
1802 XtWX = compute_banded_symmetric_XT_W_Y(X, w, X)
1803 XtE = compute_banded_symmetric_XT_W_Y(X, w, wE)
1805 def fun(lam):
1806 return _gcv(lam, X, XtWX, wE, XtE)
1808 gcv_est = minimize_scalar(fun, bounds=(0, n), method='Bounded')
1809 if gcv_est.success:
1810 return gcv_est.x
1811 raise ValueError(f"Unable to find minimum of the GCV "
1812 f"function: {gcv_est.message}")
1815def _coeff_of_divided_diff(x):
1816 """
1817 Returns the coefficients of the divided difference.
1819 Parameters
1820 ----------
1821 x : array, shape (n,)
1822 Array which is used for the computation of divided difference.
1824 Returns
1825 -------
1826 res : array_like, shape (n,)
1827 Coefficients of the divided difference.
1829 Notes
1830 -----
1831 Vector ``x`` should have unique elements, otherwise an error division by
1832 zero might be raised.
1834 No checks are performed.
1836 """
1837 n = x.shape[0]
1838 res = np.zeros(n)
1839 for i in range(n):
1840 pp = 1.
1841 for k in range(n):
1842 if k != i:
1843 pp *= (x[i] - x[k])
1844 res[i] = 1. / pp
1845 return res
1848def make_smoothing_spline(x, y, w=None, lam=None):
1849 r"""
1850 Compute the (coefficients of) smoothing cubic spline function using
1851 ``lam`` to control the tradeoff between the amount of smoothness of the
1852 curve and its proximity to the data. In case ``lam`` is None, using the
1853 GCV criteria [1] to find it.
1855 A smoothing spline is found as a solution to the regularized weighted
1856 linear regression problem:
1858 .. math::
1860 \sum\limits_{i=1}^n w_i\lvert y_i - f(x_i) \rvert^2 +
1861 \lambda\int\limits_{x_1}^{x_n} (f^{(2)}(u))^2 d u
1863 where :math:`f` is a spline function, :math:`w` is a vector of weights and
1864 :math:`\lambda` is a regularization parameter.
1866 If ``lam`` is None, we use the GCV criteria to find an optimal
1867 regularization parameter, otherwise we solve the regularized weighted
1868 linear regression problem with given parameter. The parameter controls
1869 the tradeoff in the following way: the larger the parameter becomes, the
1870 smoother the function gets.
1872 Parameters
1873 ----------
1874 x : array_like, shape (n,)
1875 Abscissas.
1876 y : array_like, shape (n,)
1877 Ordinates.
1878 w : array_like, shape (n,), optional
1879 Vector of weights. Default is ``np.ones_like(x)``.
1880 lam : float, (:math:`\lambda \geq 0`), optional
1881 Regularization parameter. If ``lam`` is None, then it is found from
1882 the GCV criteria. Default is None.
1884 Returns
1885 -------
1886 func : a BSpline object.
1887 A callable representing a spline in the B-spline basis
1888 as a solution of the problem of smoothing splines using
1889 the GCV criteria [1] in case ``lam`` is None, otherwise using the
1890 given parameter ``lam``.
1892 Notes
1893 -----
1894 This algorithm is a clean room reimplementation of the algorithm
1895 introduced by Woltring in FORTRAN [2]. The original version cannot be used
1896 in SciPy source code because of the license issues. The details of the
1897 reimplementation are discussed here (available only in Russian) [4].
1899 If the vector of weights ``w`` is None, we assume that all the points are
1900 equal in terms of weights, and vector of weights is vector of ones.
1902 Note that in weighted residual sum of squares, weights are not squared:
1903 :math:`\sum\limits_{i=1}^n w_i\lvert y_i - f(x_i) \rvert^2` while in
1904 ``splrep`` the sum is built from the squared weights.
1906 In cases when the initial problem is ill-posed (for example, the product
1907 :math:`X^T W X` where :math:`X` is a design matrix is not a positive
1908 defined matrix) a ValueError is raised.
1910 References
1911 ----------
1912 .. [1] G. Wahba, "Estimating the smoothing parameter" in Spline models for
1913 observational data, Philadelphia, Pennsylvania: Society for Industrial
1914 and Applied Mathematics, 1990, pp. 45-65.
1915 :doi:`10.1137/1.9781611970128`
1916 .. [2] H. J. Woltring, A Fortran package for generalized, cross-validatory
1917 spline smoothing and differentiation, Advances in Engineering
1918 Software, vol. 8, no. 2, pp. 104-113, 1986.
1919 :doi:`10.1016/0141-1195(86)90098-7`
1920 .. [3] T. Hastie, J. Friedman, and R. Tisbshirani, "Smoothing Splines" in
1921 The elements of Statistical Learning: Data Mining, Inference, and
1922 prediction, New York: Springer, 2017, pp. 241-249.
1923 :doi:`10.1007/978-0-387-84858-7`
1924 .. [4] E. Zemlyanoy, "Generalized cross-validation smoothing splines",
1925 BSc thesis, 2022.
1926 `<https://www.hse.ru/ba/am/students/diplomas/620910604>`_ (in
1927 Russian)
1929 Examples
1930 --------
1931 Generate some noisy data
1933 >>> import numpy as np
1934 >>> np.random.seed(1234)
1935 >>> n = 200
1936 >>> def func(x):
1937 ... return x**3 + x**2 * np.sin(4 * x)
1938 >>> x = np.sort(np.random.random_sample(n) * 4 - 2)
1939 >>> y = func(x) + np.random.normal(scale=1.5, size=n)
1941 Make a smoothing spline function
1943 >>> from scipy.interpolate import make_smoothing_spline
1944 >>> spl = make_smoothing_spline(x, y)
1946 Plot both
1948 >>> import matplotlib.pyplot as plt
1949 >>> grid = np.linspace(x[0], x[-1], 400)
1950 >>> plt.plot(grid, spl(grid), label='Spline')
1951 >>> plt.plot(grid, func(grid), label='Original function')
1952 >>> plt.scatter(x, y, marker='.')
1953 >>> plt.legend(loc='best')
1954 >>> plt.show()
1956 """
1958 x = np.ascontiguousarray(x, dtype=float)
1959 y = np.ascontiguousarray(y, dtype=float)
1961 if any(x[1:] - x[:-1] <= 0):
1962 raise ValueError('``x`` should be an ascending array')
1964 if x.ndim != 1 or y.ndim != 1 or x.shape[0] != y.shape[0]:
1965 raise ValueError('``x`` and ``y`` should be one dimensional and the'
1966 ' same size')
1968 if w is None:
1969 w = np.ones(len(x))
1970 else:
1971 w = np.ascontiguousarray(w)
1972 if any(w <= 0):
1973 raise ValueError('Invalid vector of weights')
1975 t = np.r_[[x[0]] * 3, x, [x[-1]] * 3]
1976 n = x.shape[0]
1978 # It is known that the solution to the stated minimization problem exists
1979 # and is a natural cubic spline with vector of knots equal to the unique
1980 # elements of ``x`` [3], so we will solve the problem in the basis of
1981 # natural splines.
1983 # create design matrix in the B-spline basis
1984 X_bspl = BSpline.design_matrix(x, t, 3)
1985 # move from B-spline basis to the basis of natural splines using equations
1986 # (2.1.7) [4]
1987 # central elements
1988 X = np.zeros((5, n))
1989 for i in range(1, 4):
1990 X[i, 2: -2] = X_bspl[i: i - 4, 3: -3][np.diag_indices(n - 4)]
1992 # first elements
1993 X[1, 1] = X_bspl[0, 0]
1994 X[2, :2] = ((x[2] + x[1] - 2 * x[0]) * X_bspl[0, 0],
1995 X_bspl[1, 1] + X_bspl[1, 2])
1996 X[3, :2] = ((x[2] - x[0]) * X_bspl[1, 1], X_bspl[2, 2])
1998 # last elements
1999 X[1, -2:] = (X_bspl[-3, -3], (x[-1] - x[-3]) * X_bspl[-2, -2])
2000 X[2, -2:] = (X_bspl[-2, -3] + X_bspl[-2, -2],
2001 (2 * x[-1] - x[-2] - x[-3]) * X_bspl[-1, -1])
2002 X[3, -2] = X_bspl[-1, -1]
2004 # create penalty matrix and divide it by vector of weights: W^{-1} E
2005 wE = np.zeros((5, n))
2006 wE[2:, 0] = _coeff_of_divided_diff(x[:3]) / w[:3]
2007 wE[1:, 1] = _coeff_of_divided_diff(x[:4]) / w[:4]
2008 for j in range(2, n - 2):
2009 wE[:, j] = (x[j+2] - x[j-2]) * _coeff_of_divided_diff(x[j-2:j+3])\
2010 / w[j-2: j+3]
2012 wE[:-1, -2] = -_coeff_of_divided_diff(x[-4:]) / w[-4:]
2013 wE[:-2, -1] = _coeff_of_divided_diff(x[-3:]) / w[-3:]
2014 wE *= 6
2016 if lam is None:
2017 lam = _compute_optimal_gcv_parameter(X, wE, y, w)
2018 elif lam < 0.:
2019 raise ValueError('Regularization parameter should be non-negative')
2021 # solve the initial problem in the basis of natural splines
2022 c = solve_banded((2, 2), X + lam * wE, y)
2023 # move back to B-spline basis using equations (2.2.10) [4]
2024 c_ = np.r_[c[0] * (t[5] + t[4] - 2 * t[3]) + c[1],
2025 c[0] * (t[5] - t[3]) + c[1],
2026 c[1: -1],
2027 c[-1] * (t[-4] - t[-6]) + c[-2],
2028 c[-1] * (2 * t[-4] - t[-5] - t[-6]) + c[-2]]
2030 return BSpline.construct_fast(t, c_, 3)