Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_fitpack_impl.py: 6%
413 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2fitpack (dierckx in netlib) --- A Python-C wrapper to FITPACK (by P. Dierckx).
3 FITPACK is a collection of FORTRAN programs for curve and surface
4 fitting with splines and tensor product splines.
6See
7 https://web.archive.org/web/20010524124604/http://www.cs.kuleuven.ac.be:80/cwis/research/nalag/research/topics/fitpack.html
8or
9 http://www.netlib.org/dierckx/
11Copyright 2002 Pearu Peterson all rights reserved,
12Pearu Peterson <pearu@cens.ioc.ee>
13Permission to use, modify, and distribute this software is given under the
14terms of the SciPy (BSD style) license. See LICENSE.txt that came with
15this distribution for specifics.
17NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
19TODO: Make interfaces to the following fitpack functions:
20 For univariate splines: cocosp, concon, fourco, insert
21 For bivariate splines: profil, regrid, parsur, surev
22"""
24__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
25 'bisplrep', 'bisplev', 'insert', 'splder', 'splantider']
27import warnings
28import numpy as np
29from . import _fitpack
30from numpy import (atleast_1d, array, ones, zeros, sqrt, ravel, transpose,
31 empty, iinfo, asarray)
33# Try to replace _fitpack interface with
34# f2py-generated version
35from . import dfitpack
38dfitpack_int = dfitpack.types.intvar.dtype
41def _int_overflow(x, msg=None):
42 """Cast the value to an dfitpack_int and raise an OverflowError if the value
43 cannot fit.
44 """
45 if x > iinfo(dfitpack_int).max:
46 if msg is None:
47 msg = '%r cannot fit into an %r' % (x, dfitpack_int)
48 raise OverflowError(msg)
49 return dfitpack_int.type(x)
52_iermess = {
53 0: ["The spline has a residual sum of squares fp such that "
54 "abs(fp-s)/s<=0.001", None],
55 -1: ["The spline is an interpolating spline (fp=0)", None],
56 -2: ["The spline is weighted least-squares polynomial of degree k.\n"
57 "fp gives the upper bound fp0 for the smoothing factor s", None],
58 1: ["The required storage space exceeds the available storage space.\n"
59 "Probable causes: data (x,y) size is too small or smoothing parameter"
60 "\ns is too small (fp>s).", ValueError],
61 2: ["A theoretically impossible result when finding a smoothing spline\n"
62 "with fp = s. Probable cause: s too small. (abs(fp-s)/s>0.001)",
63 ValueError],
64 3: ["The maximal number of iterations (20) allowed for finding smoothing\n"
65 "spline with fp=s has been reached. Probable cause: s too small.\n"
66 "(abs(fp-s)/s>0.001)", ValueError],
67 10: ["Error on input data", ValueError],
68 'unknown': ["An error occurred", TypeError]
69}
71_iermess2 = {
72 0: ["The spline has a residual sum of squares fp such that "
73 "abs(fp-s)/s<=0.001", None],
74 -1: ["The spline is an interpolating spline (fp=0)", None],
75 -2: ["The spline is weighted least-squares polynomial of degree kx and ky."
76 "\nfp gives the upper bound fp0 for the smoothing factor s", None],
77 -3: ["Warning. The coefficients of the spline have been computed as the\n"
78 "minimal norm least-squares solution of a rank deficient system.",
79 None],
80 1: ["The required storage space exceeds the available storage space.\n"
81 "Probable causes: nxest or nyest too small or s is too small. (fp>s)",
82 ValueError],
83 2: ["A theoretically impossible result when finding a smoothing spline\n"
84 "with fp = s. Probable causes: s too small or badly chosen eps.\n"
85 "(abs(fp-s)/s>0.001)", ValueError],
86 3: ["The maximal number of iterations (20) allowed for finding smoothing\n"
87 "spline with fp=s has been reached. Probable cause: s too small.\n"
88 "(abs(fp-s)/s>0.001)", ValueError],
89 4: ["No more knots can be added because the number of B-spline\n"
90 "coefficients already exceeds the number of data points m.\n"
91 "Probable causes: either s or m too small. (fp>s)", ValueError],
92 5: ["No more knots can be added because the additional knot would\n"
93 "coincide with an old one. Probable cause: s too small or too large\n"
94 "a weight to an inaccurate data point. (fp>s)", ValueError],
95 10: ["Error on input data", ValueError],
96 11: ["rwrk2 too small, i.e., there is not enough workspace for computing\n"
97 "the minimal least-squares solution of a rank deficient system of\n"
98 "linear equations.", ValueError],
99 'unknown': ["An error occurred", TypeError]
100}
102_parcur_cache = {'t': array([], float), 'wrk': array([], float),
103 'iwrk': array([], dfitpack_int), 'u': array([], float),
104 'ub': 0, 'ue': 1}
107def splprep(x, w=None, u=None, ub=None, ue=None, k=3, task=0, s=None, t=None,
108 full_output=0, nest=None, per=0, quiet=1):
109 """
110 Find the B-spline representation of an N-D curve.
112 Given a list of N rank-1 arrays, `x`, which represent a curve in
113 N-dimensional space parametrized by `u`, find a smooth approximating
114 spline curve g(`u`). Uses the FORTRAN routine parcur from FITPACK.
116 Parameters
117 ----------
118 x : array_like
119 A list of sample vector arrays representing the curve.
120 w : array_like, optional
121 Strictly positive rank-1 array of weights the same length as `x[0]`.
122 The weights are used in computing the weighted least-squares spline
123 fit. If the errors in the `x` values have standard-deviation given by
124 the vector d, then `w` should be 1/d. Default is ``ones(len(x[0]))``.
125 u : array_like, optional
126 An array of parameter values. If not given, these values are
127 calculated automatically as ``M = len(x[0])``, where
129 v[0] = 0
131 v[i] = v[i-1] + distance(`x[i]`, `x[i-1]`)
133 u[i] = v[i] / v[M-1]
135 ub, ue : int, optional
136 The end-points of the parameters interval. Defaults to
137 u[0] and u[-1].
138 k : int, optional
139 Degree of the spline. Cubic splines are recommended.
140 Even values of `k` should be avoided especially with a small s-value.
141 ``1 <= k <= 5``, default is 3.
142 task : int, optional
143 If task==0 (default), find t and c for a given smoothing factor, s.
144 If task==1, find t and c for another value of the smoothing factor, s.
145 There must have been a previous call with task=0 or task=1
146 for the same set of data.
147 If task=-1 find the weighted least square spline for a given set of
148 knots, t.
149 s : float, optional
150 A smoothing condition. The amount of smoothness is determined by
151 satisfying the conditions: ``sum((w * (y - g))**2,axis=0) <= s``,
152 where g(x) is the smoothed interpolation of (x,y). The user can
153 use `s` to control the trade-off between closeness and smoothness
154 of fit. Larger `s` means more smoothing while smaller values of `s`
155 indicate less smoothing. Recommended values of `s` depend on the
156 weights, w. If the weights represent the inverse of the
157 standard-deviation of y, then a good `s` value should be found in
158 the range ``(m-sqrt(2*m),m+sqrt(2*m))``, where m is the number of
159 data points in x, y, and w.
160 t : int, optional
161 The knots needed for task=-1.
162 full_output : int, optional
163 If non-zero, then return optional outputs.
164 nest : int, optional
165 An over-estimate of the total number of knots of the spline to
166 help in determining the storage space. By default nest=m/2.
167 Always large enough is nest=m+k+1.
168 per : int, optional
169 If non-zero, data points are considered periodic with period
170 ``x[m-1] - x[0]`` and a smooth periodic spline approximation is
171 returned. Values of ``y[m-1]`` and ``w[m-1]`` are not used.
172 quiet : int, optional
173 Non-zero to suppress messages.
175 Returns
176 -------
177 tck : tuple
178 A tuple (t,c,k) containing the vector of knots, the B-spline
179 coefficients, and the degree of the spline.
180 u : array
181 An array of the values of the parameter.
182 fp : float
183 The weighted sum of squared residuals of the spline approximation.
184 ier : int
185 An integer flag about splrep success. Success is indicated
186 if ier<=0. If ier in [1,2,3] an error occurred but was not raised.
187 Otherwise an error is raised.
188 msg : str
189 A message corresponding to the integer flag, ier.
191 See Also
192 --------
193 splrep, splev, sproot, spalde, splint,
194 bisplrep, bisplev
195 UnivariateSpline, BivariateSpline
197 Notes
198 -----
199 See `splev` for evaluation of the spline and its derivatives.
200 The number of dimensions N must be smaller than 11.
202 References
203 ----------
204 .. [1] P. Dierckx, "Algorithms for smoothing data with periodic and
205 parametric splines, Computer Graphics and Image Processing",
206 20 (1982) 171-184.
207 .. [2] P. Dierckx, "Algorithms for smoothing data with periodic and
208 parametric splines", report tw55, Dept. Computer Science,
209 K.U.Leuven, 1981.
210 .. [3] P. Dierckx, "Curve and surface fitting with splines", Monographs on
211 Numerical Analysis, Oxford University Press, 1993.
213 """
214 if task <= 0:
215 _parcur_cache = {'t': array([], float), 'wrk': array([], float),
216 'iwrk': array([], dfitpack_int), 'u': array([], float),
217 'ub': 0, 'ue': 1}
218 x = atleast_1d(x)
219 idim, m = x.shape
220 if per:
221 for i in range(idim):
222 if x[i][0] != x[i][-1]:
223 if not quiet:
224 warnings.warn(RuntimeWarning('Setting x[%d][%d]=x[%d][0]' %
225 (i, m, i)))
226 x[i][-1] = x[i][0]
227 if not 0 < idim < 11:
228 raise TypeError('0 < idim < 11 must hold')
229 if w is None:
230 w = ones(m, float)
231 else:
232 w = atleast_1d(w)
233 ipar = (u is not None)
234 if ipar:
235 _parcur_cache['u'] = u
236 if ub is None:
237 _parcur_cache['ub'] = u[0]
238 else:
239 _parcur_cache['ub'] = ub
240 if ue is None:
241 _parcur_cache['ue'] = u[-1]
242 else:
243 _parcur_cache['ue'] = ue
244 else:
245 _parcur_cache['u'] = zeros(m, float)
246 if not (1 <= k <= 5):
247 raise TypeError('1 <= k= %d <=5 must hold' % k)
248 if not (-1 <= task <= 1):
249 raise TypeError('task must be -1, 0 or 1')
250 if (not len(w) == m) or (ipar == 1 and (not len(u) == m)):
251 raise TypeError('Mismatch of input dimensions')
252 if s is None:
253 s = m - sqrt(2*m)
254 if t is None and task == -1:
255 raise TypeError('Knots must be given for task=-1')
256 if t is not None:
257 _parcur_cache['t'] = atleast_1d(t)
258 n = len(_parcur_cache['t'])
259 if task == -1 and n < 2*k + 2:
260 raise TypeError('There must be at least 2*k+2 knots for task=-1')
261 if m <= k:
262 raise TypeError('m > k must hold')
263 if nest is None:
264 nest = m + 2*k
266 if (task >= 0 and s == 0) or (nest < 0):
267 if per:
268 nest = m + 2*k
269 else:
270 nest = m + k + 1
271 nest = max(nest, 2*k + 3)
272 u = _parcur_cache['u']
273 ub = _parcur_cache['ub']
274 ue = _parcur_cache['ue']
275 t = _parcur_cache['t']
276 wrk = _parcur_cache['wrk']
277 iwrk = _parcur_cache['iwrk']
278 t, c, o = _fitpack._parcur(ravel(transpose(x)), w, u, ub, ue, k,
279 task, ipar, s, t, nest, wrk, iwrk, per)
280 _parcur_cache['u'] = o['u']
281 _parcur_cache['ub'] = o['ub']
282 _parcur_cache['ue'] = o['ue']
283 _parcur_cache['t'] = t
284 _parcur_cache['wrk'] = o['wrk']
285 _parcur_cache['iwrk'] = o['iwrk']
286 ier = o['ier']
287 fp = o['fp']
288 n = len(t)
289 u = o['u']
290 c.shape = idim, n - k - 1
291 tcku = [t, list(c), k], u
292 if ier <= 0 and not quiet:
293 warnings.warn(RuntimeWarning(_iermess[ier][0] +
294 "\tk=%d n=%d m=%d fp=%f s=%f" %
295 (k, len(t), m, fp, s)))
296 if ier > 0 and not full_output:
297 if ier in [1, 2, 3]:
298 warnings.warn(RuntimeWarning(_iermess[ier][0]))
299 else:
300 try:
301 raise _iermess[ier][1](_iermess[ier][0])
302 except KeyError as e:
303 raise _iermess['unknown'][1](_iermess['unknown'][0]) from e
304 if full_output:
305 try:
306 return tcku, fp, ier, _iermess[ier][0]
307 except KeyError:
308 return tcku, fp, ier, _iermess['unknown'][0]
309 else:
310 return tcku
313_curfit_cache = {'t': array([], float), 'wrk': array([], float),
314 'iwrk': array([], dfitpack_int)}
317def splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None,
318 full_output=0, per=0, quiet=1):
319 """
320 Find the B-spline representation of 1-D curve.
322 Given the set of data points ``(x[i], y[i])`` determine a smooth spline
323 approximation of degree k on the interval ``xb <= x <= xe``.
325 Parameters
326 ----------
327 x, y : array_like
328 The data points defining a curve y = f(x).
329 w : array_like, optional
330 Strictly positive rank-1 array of weights the same length as x and y.
331 The weights are used in computing the weighted least-squares spline
332 fit. If the errors in the y values have standard-deviation given by the
333 vector d, then w should be 1/d. Default is ones(len(x)).
334 xb, xe : float, optional
335 The interval to fit. If None, these default to x[0] and x[-1]
336 respectively.
337 k : int, optional
338 The order of the spline fit. It is recommended to use cubic splines.
339 Even order splines should be avoided especially with small s values.
340 1 <= k <= 5
341 task : {1, 0, -1}, optional
342 If task==0 find t and c for a given smoothing factor, s.
344 If task==1 find t and c for another value of the smoothing factor, s.
345 There must have been a previous call with task=0 or task=1 for the same
346 set of data (t will be stored an used internally)
348 If task=-1 find the weighted least square spline for a given set of
349 knots, t. These should be interior knots as knots on the ends will be
350 added automatically.
351 s : float, optional
352 A smoothing condition. The amount of smoothness is determined by
353 satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s, where g(x)
354 is the smoothed interpolation of (x,y). The user can use s to control
355 the tradeoff between closeness and smoothness of fit. Larger s means
356 more smoothing while smaller values of s indicate less smoothing.
357 Recommended values of s depend on the weights, w. If the weights
358 represent the inverse of the standard-deviation of y, then a good s
359 value should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m is
360 the number of datapoints in x, y, and w. default : s=m-sqrt(2*m) if
361 weights are supplied. s = 0.0 (interpolating) if no weights are
362 supplied.
363 t : array_like, optional
364 The knots needed for task=-1. If given then task is automatically set
365 to -1.
366 full_output : bool, optional
367 If non-zero, then return optional outputs.
368 per : bool, optional
369 If non-zero, data points are considered periodic with period x[m-1] -
370 x[0] and a smooth periodic spline approximation is returned. Values of
371 y[m-1] and w[m-1] are not used.
372 quiet : bool, optional
373 Non-zero to suppress messages.
375 Returns
376 -------
377 tck : tuple
378 (t,c,k) a tuple containing the vector of knots, the B-spline
379 coefficients, and the degree of the spline.
380 fp : array, optional
381 The weighted sum of squared residuals of the spline approximation.
382 ier : int, optional
383 An integer flag about splrep success. Success is indicated if ier<=0.
384 If ier in [1,2,3] an error occurred but was not raised. Otherwise an
385 error is raised.
386 msg : str, optional
387 A message corresponding to the integer flag, ier.
389 See Also
390 --------
391 UnivariateSpline, BivariateSpline
392 splprep, splev, sproot, spalde, splint
393 bisplrep, bisplev
395 Notes
396 -----
397 See splev for evaluation of the spline and its derivatives. Uses the
398 FORTRAN routine curfit from FITPACK.
400 The user is responsible for assuring that the values of *x* are unique.
401 Otherwise, *splrep* will not return sensible results.
403 If provided, knots `t` must satisfy the Schoenberg-Whitney conditions,
404 i.e., there must be a subset of data points ``x[j]`` such that
405 ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``.
407 References
408 ----------
409 Based on algorithms described in [1]_, [2]_, [3]_, and [4]_:
411 .. [1] P. Dierckx, "An algorithm for smoothing, differentiation and
412 integration of experimental data using spline functions",
413 J.Comp.Appl.Maths 1 (1975) 165-184.
414 .. [2] P. Dierckx, "A fast algorithm for smoothing data on a rectangular
415 grid while using spline functions", SIAM J.Numer.Anal. 19 (1982)
416 1286-1304.
417 .. [3] P. Dierckx, "An improved algorithm for curve fitting with spline
418 functions", report tw54, Dept. Computer Science,K.U. Leuven, 1981.
419 .. [4] P. Dierckx, "Curve and surface fitting with splines", Monographs on
420 Numerical Analysis, Oxford University Press, 1993.
422 Examples
423 --------
425 >>> import matplotlib.pyplot as plt
426 >>> from scipy.interpolate import splev, splrep
427 >>> x = np.linspace(0, 10, 10)
428 >>> y = np.sin(x)
429 >>> tck = splrep(x, y)
430 >>> x2 = np.linspace(0, 10, 200)
431 >>> y2 = splev(x2, tck)
432 >>> plt.plot(x, y, 'o', x2, y2)
433 >>> plt.show()
435 """
436 if task <= 0:
437 _curfit_cache = {}
438 x, y = map(atleast_1d, [x, y])
439 m = len(x)
440 if w is None:
441 w = ones(m, float)
442 if s is None:
443 s = 0.0
444 else:
445 w = atleast_1d(w)
446 if s is None:
447 s = m - sqrt(2*m)
448 if not len(w) == m:
449 raise TypeError('len(w)=%d is not equal to m=%d' % (len(w), m))
450 if (m != len(y)) or (m != len(w)):
451 raise TypeError('Lengths of the first three arguments (x,y,w) must '
452 'be equal')
453 if not (1 <= k <= 5):
454 raise TypeError('Given degree of the spline (k=%d) is not supported. '
455 '(1<=k<=5)' % k)
456 if m <= k:
457 raise TypeError('m > k must hold')
458 if xb is None:
459 xb = x[0]
460 if xe is None:
461 xe = x[-1]
462 if not (-1 <= task <= 1):
463 raise TypeError('task must be -1, 0 or 1')
464 if t is not None:
465 task = -1
466 if task == -1:
467 if t is None:
468 raise TypeError('Knots must be given for task=-1')
469 numknots = len(t)
470 _curfit_cache['t'] = empty((numknots + 2*k + 2,), float)
471 _curfit_cache['t'][k+1:-k-1] = t
472 nest = len(_curfit_cache['t'])
473 elif task == 0:
474 if per:
475 nest = max(m + 2*k, 2*k + 3)
476 else:
477 nest = max(m + k + 1, 2*k + 3)
478 t = empty((nest,), float)
479 _curfit_cache['t'] = t
480 if task <= 0:
481 if per:
482 _curfit_cache['wrk'] = empty((m*(k + 1) + nest*(8 + 5*k),), float)
483 else:
484 _curfit_cache['wrk'] = empty((m*(k + 1) + nest*(7 + 3*k),), float)
485 _curfit_cache['iwrk'] = empty((nest,), dfitpack_int)
486 try:
487 t = _curfit_cache['t']
488 wrk = _curfit_cache['wrk']
489 iwrk = _curfit_cache['iwrk']
490 except KeyError as e:
491 raise TypeError("must call with task=1 only after"
492 " call with task=0,-1") from e
493 if not per:
494 n, c, fp, ier = dfitpack.curfit(task, x, y, w, t, wrk, iwrk,
495 xb, xe, k, s)
496 else:
497 n, c, fp, ier = dfitpack.percur(task, x, y, w, t, wrk, iwrk, k, s)
498 tck = (t[:n], c[:n], k)
499 if ier <= 0 and not quiet:
500 _mess = (_iermess[ier][0] + "\tk=%d n=%d m=%d fp=%f s=%f" %
501 (k, len(t), m, fp, s))
502 warnings.warn(RuntimeWarning(_mess))
503 if ier > 0 and not full_output:
504 if ier in [1, 2, 3]:
505 warnings.warn(RuntimeWarning(_iermess[ier][0]))
506 else:
507 try:
508 raise _iermess[ier][1](_iermess[ier][0])
509 except KeyError as e:
510 raise _iermess['unknown'][1](_iermess['unknown'][0]) from e
511 if full_output:
512 try:
513 return tck, fp, ier, _iermess[ier][0]
514 except KeyError:
515 return tck, fp, ier, _iermess['unknown'][0]
516 else:
517 return tck
520def splev(x, tck, der=0, ext=0):
521 """
522 Evaluate a B-spline or its derivatives.
524 Given the knots and coefficients of a B-spline representation, evaluate
525 the value of the smoothing polynomial and its derivatives. This is a
526 wrapper around the FORTRAN routines splev and splder of FITPACK.
528 Parameters
529 ----------
530 x : array_like
531 An array of points at which to return the value of the smoothed
532 spline or its derivatives. If `tck` was returned from `splprep`,
533 then the parameter values, u should be given.
534 tck : tuple
535 A sequence of length 3 returned by `splrep` or `splprep` containing
536 the knots, coefficients, and degree of the spline.
537 der : int, optional
538 The order of derivative of the spline to compute (must be less than
539 or equal to k).
540 ext : int, optional
541 Controls the value returned for elements of ``x`` not in the
542 interval defined by the knot sequence.
544 * if ext=0, return the extrapolated value.
545 * if ext=1, return 0
546 * if ext=2, raise a ValueError
547 * if ext=3, return the boundary value.
549 The default value is 0.
551 Returns
552 -------
553 y : ndarray or list of ndarrays
554 An array of values representing the spline function evaluated at
555 the points in ``x``. If `tck` was returned from `splprep`, then this
556 is a list of arrays representing the curve in N-D space.
558 See Also
559 --------
560 splprep, splrep, sproot, spalde, splint
561 bisplrep, bisplev
563 References
564 ----------
565 .. [1] C. de Boor, "On calculating with b-splines", J. Approximation
566 Theory, 6, p.50-62, 1972.
567 .. [2] M.G. Cox, "The numerical evaluation of b-splines", J. Inst. Maths
568 Applics, 10, p.134-149, 1972.
569 .. [3] P. Dierckx, "Curve and surface fitting with splines", Monographs
570 on Numerical Analysis, Oxford University Press, 1993.
572 """
573 t, c, k = tck
574 try:
575 c[0][0]
576 parametric = True
577 except Exception:
578 parametric = False
579 if parametric:
580 return list(map(lambda c, x=x, t=t, k=k, der=der:
581 splev(x, [t, c, k], der, ext), c))
582 else:
583 if not (0 <= der <= k):
584 raise ValueError("0<=der=%d<=k=%d must hold" % (der, k))
585 if ext not in (0, 1, 2, 3):
586 raise ValueError("ext = %s not in (0, 1, 2, 3) " % ext)
588 x = asarray(x)
589 shape = x.shape
590 x = atleast_1d(x).ravel()
591 y, ier = _fitpack._spl_(x, der, t, c, k, ext)
593 if ier == 10:
594 raise ValueError("Invalid input data")
595 if ier == 1:
596 raise ValueError("Found x value not in the domain")
597 if ier:
598 raise TypeError("An error occurred")
600 return y.reshape(shape)
603def splint(a, b, tck, full_output=0):
604 """
605 Evaluate the definite integral of a B-spline.
607 Given the knots and coefficients of a B-spline, evaluate the definite
608 integral of the smoothing polynomial between two given points.
610 Parameters
611 ----------
612 a, b : float
613 The end-points of the integration interval.
614 tck : tuple
615 A tuple (t,c,k) containing the vector of knots, the B-spline
616 coefficients, and the degree of the spline (see `splev`).
617 full_output : int, optional
618 Non-zero to return optional output.
620 Returns
621 -------
622 integral : float
623 The resulting integral.
624 wrk : ndarray
625 An array containing the integrals of the normalized B-splines
626 defined on the set of knots.
628 Notes
629 -----
630 splint silently assumes that the spline function is zero outside the data
631 interval (a, b).
633 See Also
634 --------
635 splprep, splrep, sproot, spalde, splev
636 bisplrep, bisplev
637 UnivariateSpline, BivariateSpline
639 References
640 ----------
641 .. [1] P.W. Gaffney, The calculation of indefinite integrals of b-splines",
642 J. Inst. Maths Applics, 17, p.37-41, 1976.
643 .. [2] P. Dierckx, "Curve and surface fitting with splines", Monographs
644 on Numerical Analysis, Oxford University Press, 1993.
646 """
647 t, c, k = tck
648 try:
649 c[0][0]
650 parametric = True
651 except Exception:
652 parametric = False
653 if parametric:
654 return list(map(lambda c, a=a, b=b, t=t, k=k:
655 splint(a, b, [t, c, k]), c))
656 else:
657 aint, wrk = dfitpack.splint(t, c, k, a, b)
658 if full_output:
659 return aint, wrk
660 else:
661 return aint
664def sproot(tck, mest=10):
665 """
666 Find the roots of a cubic B-spline.
668 Given the knots (>=8) and coefficients of a cubic B-spline return the
669 roots of the spline.
671 Parameters
672 ----------
673 tck : tuple
674 A tuple (t,c,k) containing the vector of knots,
675 the B-spline coefficients, and the degree of the spline.
676 The number of knots must be >= 8, and the degree must be 3.
677 The knots must be a montonically increasing sequence.
678 mest : int, optional
679 An estimate of the number of zeros (Default is 10).
681 Returns
682 -------
683 zeros : ndarray
684 An array giving the roots of the spline.
686 See Also
687 --------
688 splprep, splrep, splint, spalde, splev
689 bisplrep, bisplev
690 UnivariateSpline, BivariateSpline
693 References
694 ----------
695 .. [1] C. de Boor, "On calculating with b-splines", J. Approximation
696 Theory, 6, p.50-62, 1972.
697 .. [2] M.G. Cox, "The numerical evaluation of b-splines", J. Inst. Maths
698 Applics, 10, p.134-149, 1972.
699 .. [3] P. Dierckx, "Curve and surface fitting with splines", Monographs
700 on Numerical Analysis, Oxford University Press, 1993.
702 """
703 t, c, k = tck
704 if k != 3:
705 raise ValueError("sproot works only for cubic (k=3) splines")
706 try:
707 c[0][0]
708 parametric = True
709 except Exception:
710 parametric = False
711 if parametric:
712 return list(map(lambda c, t=t, k=k, mest=mest:
713 sproot([t, c, k], mest), c))
714 else:
715 if len(t) < 8:
716 raise TypeError("The number of knots %d>=8" % len(t))
717 z, m, ier = dfitpack.sproot(t, c, mest)
718 if ier == 10:
719 raise TypeError("Invalid input data. "
720 "t1<=..<=t4<t5<..<tn-3<=..<=tn must hold.")
721 if ier == 0:
722 return z[:m]
723 if ier == 1:
724 warnings.warn(RuntimeWarning("The number of zeros exceeds mest"))
725 return z[:m]
726 raise TypeError("Unknown error")
729def spalde(x, tck):
730 """
731 Evaluate all derivatives of a B-spline.
733 Given the knots and coefficients of a cubic B-spline compute all
734 derivatives up to order k at a point (or set of points).
736 Parameters
737 ----------
738 x : array_like
739 A point or a set of points at which to evaluate the derivatives.
740 Note that ``t(k) <= x <= t(n-k+1)`` must hold for each `x`.
741 tck : tuple
742 A tuple (t,c,k) containing the vector of knots,
743 the B-spline coefficients, and the degree of the spline.
745 Returns
746 -------
747 results : {ndarray, list of ndarrays}
748 An array (or a list of arrays) containing all derivatives
749 up to order k inclusive for each point `x`.
751 See Also
752 --------
753 splprep, splrep, splint, sproot, splev, bisplrep, bisplev,
754 UnivariateSpline, BivariateSpline
756 References
757 ----------
758 .. [1] de Boor C : On calculating with b-splines, J. Approximation Theory
759 6 (1972) 50-62.
760 .. [2] Cox M.G. : The numerical evaluation of b-splines, J. Inst. Maths
761 applics 10 (1972) 134-149.
762 .. [3] Dierckx P. : Curve and surface fitting with splines, Monographs on
763 Numerical Analysis, Oxford University Press, 1993.
765 """
766 t, c, k = tck
767 try:
768 c[0][0]
769 parametric = True
770 except Exception:
771 parametric = False
772 if parametric:
773 return list(map(lambda c, x=x, t=t, k=k:
774 spalde(x, [t, c, k]), c))
775 else:
776 x = atleast_1d(x)
777 if len(x) > 1:
778 return list(map(lambda x, tck=tck: spalde(x, tck), x))
779 d, ier = dfitpack.spalde(t, c, k+1, x[0])
780 if ier == 0:
781 return d
782 if ier == 10:
783 raise TypeError("Invalid input data. t(k)<=x<=t(n-k+1) must hold.")
784 raise TypeError("Unknown error")
786# def _curfit(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,
787# full_output=0,nest=None,per=0,quiet=1):
790_surfit_cache = {'tx': array([], float), 'ty': array([], float),
791 'wrk': array([], float), 'iwrk': array([], dfitpack_int)}
794def bisplrep(x, y, z, w=None, xb=None, xe=None, yb=None, ye=None,
795 kx=3, ky=3, task=0, s=None, eps=1e-16, tx=None, ty=None,
796 full_output=0, nxest=None, nyest=None, quiet=1):
797 """
798 Find a bivariate B-spline representation of a surface.
800 Given a set of data points (x[i], y[i], z[i]) representing a surface
801 z=f(x,y), compute a B-spline representation of the surface. Based on
802 the routine SURFIT from FITPACK.
804 Parameters
805 ----------
806 x, y, z : ndarray
807 Rank-1 arrays of data points.
808 w : ndarray, optional
809 Rank-1 array of weights. By default ``w=np.ones(len(x))``.
810 xb, xe : float, optional
811 End points of approximation interval in `x`.
812 By default ``xb = x.min(), xe=x.max()``.
813 yb, ye : float, optional
814 End points of approximation interval in `y`.
815 By default ``yb=y.min(), ye = y.max()``.
816 kx, ky : int, optional
817 The degrees of the spline (1 <= kx, ky <= 5).
818 Third order (kx=ky=3) is recommended.
819 task : int, optional
820 If task=0, find knots in x and y and coefficients for a given
821 smoothing factor, s.
822 If task=1, find knots and coefficients for another value of the
823 smoothing factor, s. bisplrep must have been previously called
824 with task=0 or task=1.
825 If task=-1, find coefficients for a given set of knots tx, ty.
826 s : float, optional
827 A non-negative smoothing factor. If weights correspond
828 to the inverse of the standard-deviation of the errors in z,
829 then a good s-value should be found in the range
830 ``(m-sqrt(2*m),m+sqrt(2*m))`` where m=len(x).
831 eps : float, optional
832 A threshold for determining the effective rank of an
833 over-determined linear system of equations (0 < eps < 1).
834 `eps` is not likely to need changing.
835 tx, ty : ndarray, optional
836 Rank-1 arrays of the knots of the spline for task=-1
837 full_output : int, optional
838 Non-zero to return optional outputs.
839 nxest, nyest : int, optional
840 Over-estimates of the total number of knots. If None then
841 ``nxest = max(kx+sqrt(m/2),2*kx+3)``,
842 ``nyest = max(ky+sqrt(m/2),2*ky+3)``.
843 quiet : int, optional
844 Non-zero to suppress printing of messages.
846 Returns
847 -------
848 tck : array_like
849 A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and
850 coefficients (c) of the bivariate B-spline representation of the
851 surface along with the degree of the spline.
852 fp : ndarray
853 The weighted sum of squared residuals of the spline approximation.
854 ier : int
855 An integer flag about splrep success. Success is indicated if
856 ier<=0. If ier in [1,2,3] an error occurred but was not raised.
857 Otherwise an error is raised.
858 msg : str
859 A message corresponding to the integer flag, ier.
861 See Also
862 --------
863 splprep, splrep, splint, sproot, splev
864 UnivariateSpline, BivariateSpline
866 Notes
867 -----
868 See `bisplev` to evaluate the value of the B-spline given its tck
869 representation.
871 If the input data is such that input dimensions have incommensurate
872 units and differ by many orders of magnitude, the interpolant may have
873 numerical artifacts. Consider rescaling the data before interpolation.
875 References
876 ----------
877 .. [1] Dierckx P.:An algorithm for surface fitting with spline functions
878 Ima J. Numer. Anal. 1 (1981) 267-283.
879 .. [2] Dierckx P.:An algorithm for surface fitting with spline functions
880 report tw50, Dept. Computer Science,K.U.Leuven, 1980.
881 .. [3] Dierckx P.:Curve and surface fitting with splines, Monographs on
882 Numerical Analysis, Oxford University Press, 1993.
884 Examples
885 --------
886 Examples are given :ref:`in the tutorial <tutorial-interpolate_2d_spline>`.
888 """
889 x, y, z = map(ravel, [x, y, z]) # ensure 1-d arrays.
890 m = len(x)
891 if not (m == len(y) == len(z)):
892 raise TypeError('len(x)==len(y)==len(z) must hold.')
893 if w is None:
894 w = ones(m, float)
895 else:
896 w = atleast_1d(w)
897 if not len(w) == m:
898 raise TypeError('len(w)=%d is not equal to m=%d' % (len(w), m))
899 if xb is None:
900 xb = x.min()
901 if xe is None:
902 xe = x.max()
903 if yb is None:
904 yb = y.min()
905 if ye is None:
906 ye = y.max()
907 if not (-1 <= task <= 1):
908 raise TypeError('task must be -1, 0 or 1')
909 if s is None:
910 s = m - sqrt(2*m)
911 if tx is None and task == -1:
912 raise TypeError('Knots_x must be given for task=-1')
913 if tx is not None:
914 _surfit_cache['tx'] = atleast_1d(tx)
915 nx = len(_surfit_cache['tx'])
916 if ty is None and task == -1:
917 raise TypeError('Knots_y must be given for task=-1')
918 if ty is not None:
919 _surfit_cache['ty'] = atleast_1d(ty)
920 ny = len(_surfit_cache['ty'])
921 if task == -1 and nx < 2*kx+2:
922 raise TypeError('There must be at least 2*kx+2 knots_x for task=-1')
923 if task == -1 and ny < 2*ky+2:
924 raise TypeError('There must be at least 2*ky+2 knots_x for task=-1')
925 if not ((1 <= kx <= 5) and (1 <= ky <= 5)):
926 raise TypeError('Given degree of the spline (kx,ky=%d,%d) is not '
927 'supported. (1<=k<=5)' % (kx, ky))
928 if m < (kx + 1)*(ky + 1):
929 raise TypeError('m >= (kx+1)(ky+1) must hold')
930 if nxest is None:
931 nxest = int(kx + sqrt(m/2))
932 if nyest is None:
933 nyest = int(ky + sqrt(m/2))
934 nxest, nyest = max(nxest, 2*kx + 3), max(nyest, 2*ky + 3)
935 if task >= 0 and s == 0:
936 nxest = int(kx + sqrt(3*m))
937 nyest = int(ky + sqrt(3*m))
938 if task == -1:
939 _surfit_cache['tx'] = atleast_1d(tx)
940 _surfit_cache['ty'] = atleast_1d(ty)
941 tx, ty = _surfit_cache['tx'], _surfit_cache['ty']
942 wrk = _surfit_cache['wrk']
943 u = nxest - kx - 1
944 v = nyest - ky - 1
945 km = max(kx, ky) + 1
946 ne = max(nxest, nyest)
947 bx, by = kx*v + ky + 1, ky*u + kx + 1
948 b1, b2 = bx, bx + v - ky
949 if bx > by:
950 b1, b2 = by, by + u - kx
951 msg = "Too many data points to interpolate"
952 lwrk1 = _int_overflow(u*v*(2 + b1 + b2) +
953 2*(u + v + km*(m + ne) + ne - kx - ky) + b2 + 1,
954 msg=msg)
955 lwrk2 = _int_overflow(u*v*(b2 + 1) + b2, msg=msg)
956 tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
957 task, s, eps, tx, ty, nxest, nyest,
958 wrk, lwrk1, lwrk2)
959 _curfit_cache['tx'] = tx
960 _curfit_cache['ty'] = ty
961 _curfit_cache['wrk'] = o['wrk']
962 ier, fp = o['ier'], o['fp']
963 tck = [tx, ty, c, kx, ky]
965 ierm = min(11, max(-3, ier))
966 if ierm <= 0 and not quiet:
967 _mess = (_iermess2[ierm][0] +
968 "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f" %
969 (kx, ky, len(tx), len(ty), m, fp, s))
970 warnings.warn(RuntimeWarning(_mess))
971 if ierm > 0 and not full_output:
972 if ier in [1, 2, 3, 4, 5]:
973 _mess = ("\n\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f" %
974 (kx, ky, len(tx), len(ty), m, fp, s))
975 warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))
976 else:
977 try:
978 raise _iermess2[ierm][1](_iermess2[ierm][0])
979 except KeyError as e:
980 raise _iermess2['unknown'][1](_iermess2['unknown'][0]) from e
981 if full_output:
982 try:
983 return tck, fp, ier, _iermess2[ierm][0]
984 except KeyError:
985 return tck, fp, ier, _iermess2['unknown'][0]
986 else:
987 return tck
990def bisplev(x, y, tck, dx=0, dy=0):
991 """
992 Evaluate a bivariate B-spline and its derivatives.
994 Return a rank-2 array of spline function values (or spline derivative
995 values) at points given by the cross-product of the rank-1 arrays `x` and
996 `y`. In special cases, return an array or just a float if either `x` or
997 `y` or both are floats. Based on BISPEV and PARDER from FITPACK.
999 Parameters
1000 ----------
1001 x, y : ndarray
1002 Rank-1 arrays specifying the domain over which to evaluate the
1003 spline or its derivative.
1004 tck : tuple
1005 A sequence of length 5 returned by `bisplrep` containing the knot
1006 locations, the coefficients, and the degree of the spline:
1007 [tx, ty, c, kx, ky].
1008 dx, dy : int, optional
1009 The orders of the partial derivatives in `x` and `y` respectively.
1011 Returns
1012 -------
1013 vals : ndarray
1014 The B-spline or its derivative evaluated over the set formed by
1015 the cross-product of `x` and `y`.
1017 See Also
1018 --------
1019 splprep, splrep, splint, sproot, splev
1020 UnivariateSpline, BivariateSpline
1022 Notes
1023 -----
1024 See `bisplrep` to generate the `tck` representation.
1026 References
1027 ----------
1028 .. [1] Dierckx P. : An algorithm for surface fitting
1029 with spline functions
1030 Ima J. Numer. Anal. 1 (1981) 267-283.
1031 .. [2] Dierckx P. : An algorithm for surface fitting
1032 with spline functions
1033 report tw50, Dept. Computer Science,K.U.Leuven, 1980.
1034 .. [3] Dierckx P. : Curve and surface fitting with splines,
1035 Monographs on Numerical Analysis, Oxford University Press, 1993.
1037 Examples
1038 --------
1039 Examples are given :ref:`in the tutorial <tutorial-interpolate_2d_spline>`.
1041 """
1042 tx, ty, c, kx, ky = tck
1043 if not (0 <= dx < kx):
1044 raise ValueError("0 <= dx = %d < kx = %d must hold" % (dx, kx))
1045 if not (0 <= dy < ky):
1046 raise ValueError("0 <= dy = %d < ky = %d must hold" % (dy, ky))
1047 x, y = map(atleast_1d, [x, y])
1048 if (len(x.shape) != 1) or (len(y.shape) != 1):
1049 raise ValueError("First two entries should be rank-1 arrays.")
1050 z, ier = _fitpack._bispev(tx, ty, c, kx, ky, x, y, dx, dy)
1051 if ier == 10:
1052 raise ValueError("Invalid input data")
1053 if ier:
1054 raise TypeError("An error occurred")
1055 z.shape = len(x), len(y)
1056 if len(z) > 1:
1057 return z
1058 if len(z[0]) > 1:
1059 return z[0]
1060 return z[0][0]
1063def dblint(xa, xb, ya, yb, tck):
1064 """Evaluate the integral of a spline over area [xa,xb] x [ya,yb].
1066 Parameters
1067 ----------
1068 xa, xb : float
1069 The end-points of the x integration interval.
1070 ya, yb : float
1071 The end-points of the y integration interval.
1072 tck : list [tx, ty, c, kx, ky]
1073 A sequence of length 5 returned by bisplrep containing the knot
1074 locations tx, ty, the coefficients c, and the degrees kx, ky
1075 of the spline.
1077 Returns
1078 -------
1079 integ : float
1080 The value of the resulting integral.
1081 """
1082 tx, ty, c, kx, ky = tck
1083 return dfitpack.dblint(tx, ty, c, kx, ky, xa, xb, ya, yb)
1086def insert(x, tck, m=1, per=0):
1087 """
1088 Insert knots into a B-spline.
1090 Given the knots and coefficients of a B-spline representation, create a
1091 new B-spline with a knot inserted `m` times at point `x`.
1092 This is a wrapper around the FORTRAN routine insert of FITPACK.
1094 Parameters
1095 ----------
1096 x (u) : array_like
1097 A 1-D point at which to insert a new knot(s). If `tck` was returned
1098 from ``splprep``, then the parameter values, u should be given.
1099 tck : tuple
1100 A tuple (t,c,k) returned by ``splrep`` or ``splprep`` containing
1101 the vector of knots, the B-spline coefficients,
1102 and the degree of the spline.
1103 m : int, optional
1104 The number of times to insert the given knot (its multiplicity).
1105 Default is 1.
1106 per : int, optional
1107 If non-zero, the input spline is considered periodic.
1109 Returns
1110 -------
1111 tck : tuple
1112 A tuple (t,c,k) containing the vector of knots, the B-spline
1113 coefficients, and the degree of the new spline.
1114 ``t(k+1) <= x <= t(n-k)``, where k is the degree of the spline.
1115 In case of a periodic spline (``per != 0``) there must be
1116 either at least k interior knots t(j) satisfying ``t(k+1)<t(j)<=x``
1117 or at least k interior knots t(j) satisfying ``x<=t(j)<t(n-k)``.
1119 Notes
1120 -----
1121 Based on algorithms from [1]_ and [2]_.
1123 References
1124 ----------
1125 .. [1] W. Boehm, "Inserting new knots into b-spline curves.",
1126 Computer Aided Design, 12, p.199-201, 1980.
1127 .. [2] P. Dierckx, "Curve and surface fitting with splines, Monographs on
1128 Numerical Analysis", Oxford University Press, 1993.
1130 """
1131 t, c, k = tck
1132 try:
1133 c[0][0]
1134 parametric = True
1135 except Exception:
1136 parametric = False
1137 if parametric:
1138 cc = []
1139 for c_vals in c:
1140 tt, cc_val, kk = insert(x, [t, c_vals, k], m)
1141 cc.append(cc_val)
1142 return (tt, cc, kk)
1143 else:
1144 tt, cc, ier = _fitpack._insert(per, t, c, k, x, m)
1145 if ier == 10:
1146 raise ValueError("Invalid input data")
1147 if ier:
1148 raise TypeError("An error occurred")
1149 return (tt, cc, k)
1152def splder(tck, n=1):
1153 """
1154 Compute the spline representation of the derivative of a given spline
1156 Parameters
1157 ----------
1158 tck : tuple of (t, c, k)
1159 Spline whose derivative to compute
1160 n : int, optional
1161 Order of derivative to evaluate. Default: 1
1163 Returns
1164 -------
1165 tck_der : tuple of (t2, c2, k2)
1166 Spline of order k2=k-n representing the derivative
1167 of the input spline.
1169 Notes
1170 -----
1172 .. versionadded:: 0.13.0
1174 See Also
1175 --------
1176 splantider, splev, spalde
1178 Examples
1179 --------
1180 This can be used for finding maxima of a curve:
1182 >>> from scipy.interpolate import splrep, splder, sproot
1183 >>> x = np.linspace(0, 10, 70)
1184 >>> y = np.sin(x)
1185 >>> spl = splrep(x, y, k=4)
1187 Now, differentiate the spline and find the zeros of the
1188 derivative. (NB: `sproot` only works for order 3 splines, so we
1189 fit an order 4 spline):
1191 >>> dspl = splder(spl)
1192 >>> sproot(dspl) / np.pi
1193 array([ 0.50000001, 1.5 , 2.49999998])
1195 This agrees well with roots :math:`\\pi/2 + n\\pi` of
1196 :math:`\\cos(x) = \\sin'(x)`.
1198 """
1199 if n < 0:
1200 return splantider(tck, -n)
1202 t, c, k = tck
1204 if n > k:
1205 raise ValueError(("Order of derivative (n = %r) must be <= "
1206 "order of spline (k = %r)") % (n, tck[2]))
1208 # Extra axes for the trailing dims of the `c` array:
1209 sh = (slice(None),) + ((None,)*len(c.shape[1:]))
1211 with np.errstate(invalid='raise', divide='raise'):
1212 try:
1213 for j in range(n):
1214 # See e.g. Schumaker, Spline Functions: Basic Theory, Chapter 5
1216 # Compute the denominator in the differentiation formula.
1217 # (and append traling dims, if necessary)
1218 dt = t[k+1:-1] - t[1:-k-1]
1219 dt = dt[sh]
1220 # Compute the new coefficients
1221 c = (c[1:-1-k] - c[:-2-k]) * k / dt
1222 # Pad coefficient array to same size as knots (FITPACK
1223 # convention)
1224 c = np.r_[c, np.zeros((k,) + c.shape[1:])]
1225 # Adjust knots
1226 t = t[1:-1]
1227 k -= 1
1228 except FloatingPointError as e:
1229 raise ValueError(("The spline has internal repeated knots "
1230 "and is not differentiable %d times") % n) from e
1232 return t, c, k
1235def splantider(tck, n=1):
1236 """
1237 Compute the spline for the antiderivative (integral) of a given spline.
1239 Parameters
1240 ----------
1241 tck : tuple of (t, c, k)
1242 Spline whose antiderivative to compute
1243 n : int, optional
1244 Order of antiderivative to evaluate. Default: 1
1246 Returns
1247 -------
1248 tck_ader : tuple of (t2, c2, k2)
1249 Spline of order k2=k+n representing the antiderivative of the input
1250 spline.
1252 See Also
1253 --------
1254 splder, splev, spalde
1256 Notes
1257 -----
1258 The `splder` function is the inverse operation of this function.
1259 Namely, ``splder(splantider(tck))`` is identical to `tck`, modulo
1260 rounding error.
1262 .. versionadded:: 0.13.0
1264 Examples
1265 --------
1266 >>> from scipy.interpolate import splrep, splder, splantider, splev
1267 >>> x = np.linspace(0, np.pi/2, 70)
1268 >>> y = 1 / np.sqrt(1 - 0.8*np.sin(x)**2)
1269 >>> spl = splrep(x, y)
1271 The derivative is the inverse operation of the antiderivative,
1272 although some floating point error accumulates:
1274 >>> splev(1.7, spl), splev(1.7, splder(splantider(spl)))
1275 (array(2.1565429877197317), array(2.1565429877201865))
1277 Antiderivative can be used to evaluate definite integrals:
1279 >>> ispl = splantider(spl)
1280 >>> splev(np.pi/2, ispl) - splev(0, ispl)
1281 2.2572053588768486
1283 This is indeed an approximation to the complete elliptic integral
1284 :math:`K(m) = \\int_0^{\\pi/2} [1 - m\\sin^2 x]^{-1/2} dx`:
1286 >>> from scipy.special import ellipk
1287 >>> ellipk(0.8)
1288 2.2572053268208538
1290 """
1291 if n < 0:
1292 return splder(tck, -n)
1294 t, c, k = tck
1296 # Extra axes for the trailing dims of the `c` array:
1297 sh = (slice(None),) + (None,)*len(c.shape[1:])
1299 for j in range(n):
1300 # This is the inverse set of operations to splder.
1302 # Compute the multiplier in the antiderivative formula.
1303 dt = t[k+1:] - t[:-k-1]
1304 dt = dt[sh]
1305 # Compute the new coefficients
1306 c = np.cumsum(c[:-k-1] * dt, axis=0) / (k + 1)
1307 c = np.r_[np.zeros((1,) + c.shape[1:]),
1308 c,
1309 [c[-1]] * (k+2)]
1310 # New knots
1311 t = np.r_[t[0], t, t[-1]]
1312 k += 1
1314 return t, c, k