_fitpack_py.py 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788
  1. __all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
  2. 'bisplrep', 'bisplev', 'insert', 'splder', 'splantider']
  3. import numpy as np
  4. # These are in the API for fitpack even if not used in fitpack.py itself.
  5. from ._fitpack_impl import bisplrep, bisplev, dblint
  6. from . import _fitpack_impl as _impl
  7. from ._bsplines import BSpline
  8. def splprep(x, w=None, u=None, ub=None, ue=None, k=3, task=0, s=None, t=None,
  9. full_output=0, nest=None, per=0, quiet=1):
  10. """
  11. Find the B-spline representation of an N-D curve.
  12. Given a list of N rank-1 arrays, `x`, which represent a curve in
  13. N-D space parametrized by `u`, find a smooth approximating
  14. spline curve g(`u`). Uses the FORTRAN routine parcur from FITPACK.
  15. Parameters
  16. ----------
  17. x : array_like
  18. A list of sample vector arrays representing the curve.
  19. w : array_like, optional
  20. Strictly positive rank-1 array of weights the same length as `x[0]`.
  21. The weights are used in computing the weighted least-squares spline
  22. fit. If the errors in the `x` values have standard-deviation given by
  23. the vector d, then `w` should be 1/d. Default is ``ones(len(x[0]))``.
  24. u : array_like, optional
  25. An array of parameter values. If not given, these values are
  26. calculated automatically as ``M = len(x[0])``, where
  27. v[0] = 0
  28. v[i] = v[i-1] + distance(`x[i]`, `x[i-1]`)
  29. u[i] = v[i] / v[M-1]
  30. ub, ue : int, optional
  31. The end-points of the parameters interval. Defaults to
  32. u[0] and u[-1].
  33. k : int, optional
  34. Degree of the spline. Cubic splines are recommended.
  35. Even values of `k` should be avoided especially with a small s-value.
  36. ``1 <= k <= 5``, default is 3.
  37. task : int, optional
  38. If task==0 (default), find t and c for a given smoothing factor, s.
  39. If task==1, find t and c for another value of the smoothing factor, s.
  40. There must have been a previous call with task=0 or task=1
  41. for the same set of data.
  42. If task=-1 find the weighted least square spline for a given set of
  43. knots, t.
  44. s : float, optional
  45. A smoothing condition. The amount of smoothness is determined by
  46. satisfying the conditions: ``sum((w * (y - g))**2,axis=0) <= s``,
  47. where g(x) is the smoothed interpolation of (x,y). The user can
  48. use `s` to control the trade-off between closeness and smoothness
  49. of fit. Larger `s` means more smoothing while smaller values of `s`
  50. indicate less smoothing. Recommended values of `s` depend on the
  51. weights, w. If the weights represent the inverse of the
  52. standard-deviation of y, then a good `s` value should be found in
  53. the range ``(m-sqrt(2*m),m+sqrt(2*m))``, where m is the number of
  54. data points in x, y, and w.
  55. t : int, optional
  56. The knots needed for task=-1.
  57. full_output : int, optional
  58. If non-zero, then return optional outputs.
  59. nest : int, optional
  60. An over-estimate of the total number of knots of the spline to
  61. help in determining the storage space. By default nest=m/2.
  62. Always large enough is nest=m+k+1.
  63. per : int, optional
  64. If non-zero, data points are considered periodic with period
  65. ``x[m-1] - x[0]`` and a smooth periodic spline approximation is
  66. returned. Values of ``y[m-1]`` and ``w[m-1]`` are not used.
  67. quiet : int, optional
  68. Non-zero to suppress messages.
  69. Returns
  70. -------
  71. tck : tuple
  72. (t,c,k) a tuple containing the vector of knots, the B-spline
  73. coefficients, and the degree of the spline.
  74. u : array
  75. An array of the values of the parameter.
  76. fp : float
  77. The weighted sum of squared residuals of the spline approximation.
  78. ier : int
  79. An integer flag about splrep success. Success is indicated
  80. if ier<=0. If ier in [1,2,3] an error occurred but was not raised.
  81. Otherwise an error is raised.
  82. msg : str
  83. A message corresponding to the integer flag, ier.
  84. See Also
  85. --------
  86. splrep, splev, sproot, spalde, splint,
  87. bisplrep, bisplev
  88. UnivariateSpline, BivariateSpline
  89. BSpline
  90. make_interp_spline
  91. Notes
  92. -----
  93. See `splev` for evaluation of the spline and its derivatives.
  94. The number of dimensions N must be smaller than 11.
  95. The number of coefficients in the `c` array is ``k+1`` less than the number
  96. of knots, ``len(t)``. This is in contrast with `splrep`, which zero-pads
  97. the array of coefficients to have the same length as the array of knots.
  98. These additional coefficients are ignored by evaluation routines, `splev`
  99. and `BSpline`.
  100. References
  101. ----------
  102. .. [1] P. Dierckx, "Algorithms for smoothing data with periodic and
  103. parametric splines, Computer Graphics and Image Processing",
  104. 20 (1982) 171-184.
  105. .. [2] P. Dierckx, "Algorithms for smoothing data with periodic and
  106. parametric splines", report tw55, Dept. Computer Science,
  107. K.U.Leuven, 1981.
  108. .. [3] P. Dierckx, "Curve and surface fitting with splines", Monographs on
  109. Numerical Analysis, Oxford University Press, 1993.
  110. Examples
  111. --------
  112. Generate a discretization of a limacon curve in the polar coordinates:
  113. >>> import numpy as np
  114. >>> phi = np.linspace(0, 2.*np.pi, 40)
  115. >>> r = 0.5 + np.cos(phi) # polar coords
  116. >>> x, y = r * np.cos(phi), r * np.sin(phi) # convert to cartesian
  117. And interpolate:
  118. >>> from scipy.interpolate import splprep, splev
  119. >>> tck, u = splprep([x, y], s=0)
  120. >>> new_points = splev(u, tck)
  121. Notice that (i) we force interpolation by using `s=0`,
  122. (ii) the parameterization, ``u``, is generated automatically.
  123. Now plot the result:
  124. >>> import matplotlib.pyplot as plt
  125. >>> fig, ax = plt.subplots()
  126. >>> ax.plot(x, y, 'ro')
  127. >>> ax.plot(new_points[0], new_points[1], 'r-')
  128. >>> plt.show()
  129. """
  130. res = _impl.splprep(x, w, u, ub, ue, k, task, s, t, full_output, nest, per,
  131. quiet)
  132. return res
  133. def splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None,
  134. full_output=0, per=0, quiet=1):
  135. """
  136. Find the B-spline representation of a 1-D curve.
  137. Given the set of data points ``(x[i], y[i])`` determine a smooth spline
  138. approximation of degree k on the interval ``xb <= x <= xe``.
  139. Parameters
  140. ----------
  141. x, y : array_like
  142. The data points defining a curve y = f(x).
  143. w : array_like, optional
  144. Strictly positive rank-1 array of weights the same length as x and y.
  145. The weights are used in computing the weighted least-squares spline
  146. fit. If the errors in the y values have standard-deviation given by the
  147. vector d, then w should be 1/d. Default is ones(len(x)).
  148. xb, xe : float, optional
  149. The interval to fit. If None, these default to x[0] and x[-1]
  150. respectively.
  151. k : int, optional
  152. The degree of the spline fit. It is recommended to use cubic splines.
  153. Even values of k should be avoided especially with small s values.
  154. 1 <= k <= 5
  155. task : {1, 0, -1}, optional
  156. If task==0 find t and c for a given smoothing factor, s.
  157. If task==1 find t and c for another value of the smoothing factor, s.
  158. There must have been a previous call with task=0 or task=1 for the same
  159. set of data (t will be stored an used internally)
  160. If task=-1 find the weighted least square spline for a given set of
  161. knots, t. These should be interior knots as knots on the ends will be
  162. added automatically.
  163. s : float, optional
  164. A smoothing condition. The amount of smoothness is determined by
  165. satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where g(x)
  166. is the smoothed interpolation of (x,y). The user can use s to control
  167. the tradeoff between closeness and smoothness of fit. Larger s means
  168. more smoothing while smaller values of s indicate less smoothing.
  169. Recommended values of s depend on the weights, w. If the weights
  170. represent the inverse of the standard-deviation of y, then a good s
  171. value should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m is
  172. the number of datapoints in x, y, and w. default : s=m-sqrt(2*m) if
  173. weights are supplied. s = 0.0 (interpolating) if no weights are
  174. supplied.
  175. t : array_like, optional
  176. The knots needed for task=-1. If given then task is automatically set
  177. to -1.
  178. full_output : bool, optional
  179. If non-zero, then return optional outputs.
  180. per : bool, optional
  181. If non-zero, data points are considered periodic with period x[m-1] -
  182. x[0] and a smooth periodic spline approximation is returned. Values of
  183. y[m-1] and w[m-1] are not used.
  184. quiet : bool, optional
  185. Non-zero to suppress messages.
  186. Returns
  187. -------
  188. tck : tuple
  189. A tuple (t,c,k) containing the vector of knots, the B-spline
  190. coefficients, and the degree of the spline.
  191. fp : array, optional
  192. The weighted sum of squared residuals of the spline approximation.
  193. ier : int, optional
  194. An integer flag about splrep success. Success is indicated if ier<=0.
  195. If ier in [1,2,3] an error occurred but was not raised. Otherwise an
  196. error is raised.
  197. msg : str, optional
  198. A message corresponding to the integer flag, ier.
  199. See Also
  200. --------
  201. UnivariateSpline, BivariateSpline
  202. splprep, splev, sproot, spalde, splint
  203. bisplrep, bisplev
  204. BSpline
  205. make_interp_spline
  206. Notes
  207. -----
  208. See `splev` for evaluation of the spline and its derivatives. Uses the
  209. FORTRAN routine ``curfit`` from FITPACK.
  210. The user is responsible for assuring that the values of `x` are unique.
  211. Otherwise, `splrep` will not return sensible results.
  212. If provided, knots `t` must satisfy the Schoenberg-Whitney conditions,
  213. i.e., there must be a subset of data points ``x[j]`` such that
  214. ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``.
  215. This routine zero-pads the coefficients array ``c`` to have the same length
  216. as the array of knots ``t`` (the trailing ``k + 1`` coefficients are ignored
  217. by the evaluation routines, `splev` and `BSpline`.) This is in contrast with
  218. `splprep`, which does not zero-pad the coefficients.
  219. References
  220. ----------
  221. Based on algorithms described in [1]_, [2]_, [3]_, and [4]_:
  222. .. [1] P. Dierckx, "An algorithm for smoothing, differentiation and
  223. integration of experimental data using spline functions",
  224. J.Comp.Appl.Maths 1 (1975) 165-184.
  225. .. [2] P. Dierckx, "A fast algorithm for smoothing data on a rectangular
  226. grid while using spline functions", SIAM J.Numer.Anal. 19 (1982)
  227. 1286-1304.
  228. .. [3] P. Dierckx, "An improved algorithm for curve fitting with spline
  229. functions", report tw54, Dept. Computer Science,K.U. Leuven, 1981.
  230. .. [4] P. Dierckx, "Curve and surface fitting with splines", Monographs on
  231. Numerical Analysis, Oxford University Press, 1993.
  232. Examples
  233. --------
  234. You can interpolate 1-D points with a B-spline curve.
  235. Further examples are given in
  236. :ref:`in the tutorial <tutorial-interpolate_splXXX>`.
  237. >>> import numpy as np
  238. >>> import matplotlib.pyplot as plt
  239. >>> from scipy.interpolate import splev, splrep
  240. >>> x = np.linspace(0, 10, 10)
  241. >>> y = np.sin(x)
  242. >>> spl = splrep(x, y)
  243. >>> x2 = np.linspace(0, 10, 200)
  244. >>> y2 = splev(x2, spl)
  245. >>> plt.plot(x, y, 'o', x2, y2)
  246. >>> plt.show()
  247. """
  248. res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
  249. return res
  250. def splev(x, tck, der=0, ext=0):
  251. """
  252. Evaluate a B-spline or its derivatives.
  253. Given the knots and coefficients of a B-spline representation, evaluate
  254. the value of the smoothing polynomial and its derivatives. This is a
  255. wrapper around the FORTRAN routines splev and splder of FITPACK.
  256. Parameters
  257. ----------
  258. x : array_like
  259. An array of points at which to return the value of the smoothed
  260. spline or its derivatives. If `tck` was returned from `splprep`,
  261. then the parameter values, u should be given.
  262. tck : 3-tuple or a BSpline object
  263. If a tuple, then it should be a sequence of length 3 returned by
  264. `splrep` or `splprep` containing the knots, coefficients, and degree
  265. of the spline. (Also see Notes.)
  266. der : int, optional
  267. The order of derivative of the spline to compute (must be less than
  268. or equal to k, the degree of the spline).
  269. ext : int, optional
  270. Controls the value returned for elements of ``x`` not in the
  271. interval defined by the knot sequence.
  272. * if ext=0, return the extrapolated value.
  273. * if ext=1, return 0
  274. * if ext=2, raise a ValueError
  275. * if ext=3, return the boundary value.
  276. The default value is 0.
  277. Returns
  278. -------
  279. y : ndarray or list of ndarrays
  280. An array of values representing the spline function evaluated at
  281. the points in `x`. If `tck` was returned from `splprep`, then this
  282. is a list of arrays representing the curve in an N-D space.
  283. Notes
  284. -----
  285. Manipulating the tck-tuples directly is not recommended. In new code,
  286. prefer using `BSpline` objects.
  287. See Also
  288. --------
  289. splprep, splrep, sproot, spalde, splint
  290. bisplrep, bisplev
  291. BSpline
  292. References
  293. ----------
  294. .. [1] C. de Boor, "On calculating with b-splines", J. Approximation
  295. Theory, 6, p.50-62, 1972.
  296. .. [2] M. G. Cox, "The numerical evaluation of b-splines", J. Inst. Maths
  297. Applics, 10, p.134-149, 1972.
  298. .. [3] P. Dierckx, "Curve and surface fitting with splines", Monographs
  299. on Numerical Analysis, Oxford University Press, 1993.
  300. Examples
  301. --------
  302. Examples are given :ref:`in the tutorial <tutorial-interpolate_splXXX>`.
  303. """
  304. if isinstance(tck, BSpline):
  305. if tck.c.ndim > 1:
  306. mesg = ("Calling splev() with BSpline objects with c.ndim > 1 is "
  307. "not allowed. Use BSpline.__call__(x) instead.")
  308. raise ValueError(mesg)
  309. # remap the out-of-bounds behavior
  310. try:
  311. extrapolate = {0: True, }[ext]
  312. except KeyError as e:
  313. raise ValueError("Extrapolation mode %s is not supported "
  314. "by BSpline." % ext) from e
  315. return tck(x, der, extrapolate=extrapolate)
  316. else:
  317. return _impl.splev(x, tck, der, ext)
  318. def splint(a, b, tck, full_output=0):
  319. """
  320. Evaluate the definite integral of a B-spline between two given points.
  321. Parameters
  322. ----------
  323. a, b : float
  324. The end-points of the integration interval.
  325. tck : tuple or a BSpline instance
  326. If a tuple, then it should be a sequence of length 3, containing the
  327. vector of knots, the B-spline coefficients, and the degree of the
  328. spline (see `splev`).
  329. full_output : int, optional
  330. Non-zero to return optional output.
  331. Returns
  332. -------
  333. integral : float
  334. The resulting integral.
  335. wrk : ndarray
  336. An array containing the integrals of the normalized B-splines
  337. defined on the set of knots.
  338. (Only returned if `full_output` is non-zero)
  339. Notes
  340. -----
  341. `splint` silently assumes that the spline function is zero outside the data
  342. interval (`a`, `b`).
  343. Manipulating the tck-tuples directly is not recommended. In new code,
  344. prefer using the `BSpline` objects.
  345. See Also
  346. --------
  347. splprep, splrep, sproot, spalde, splev
  348. bisplrep, bisplev
  349. BSpline
  350. References
  351. ----------
  352. .. [1] P.W. Gaffney, The calculation of indefinite integrals of b-splines",
  353. J. Inst. Maths Applics, 17, p.37-41, 1976.
  354. .. [2] P. Dierckx, "Curve and surface fitting with splines", Monographs
  355. on Numerical Analysis, Oxford University Press, 1993.
  356. Examples
  357. --------
  358. Examples are given :ref:`in the tutorial <tutorial-interpolate_splXXX>`.
  359. """
  360. if isinstance(tck, BSpline):
  361. if tck.c.ndim > 1:
  362. mesg = ("Calling splint() with BSpline objects with c.ndim > 1 is "
  363. "not allowed. Use BSpline.integrate() instead.")
  364. raise ValueError(mesg)
  365. if full_output != 0:
  366. mesg = ("full_output = %s is not supported. Proceeding as if "
  367. "full_output = 0" % full_output)
  368. return tck.integrate(a, b, extrapolate=False)
  369. else:
  370. return _impl.splint(a, b, tck, full_output)
  371. def sproot(tck, mest=10):
  372. """
  373. Find the roots of a cubic B-spline.
  374. Given the knots (>=8) and coefficients of a cubic B-spline return the
  375. roots of the spline.
  376. Parameters
  377. ----------
  378. tck : tuple or a BSpline object
  379. If a tuple, then it should be a sequence of length 3, containing the
  380. vector of knots, the B-spline coefficients, and the degree of the
  381. spline.
  382. The number of knots must be >= 8, and the degree must be 3.
  383. The knots must be a montonically increasing sequence.
  384. mest : int, optional
  385. An estimate of the number of zeros (Default is 10).
  386. Returns
  387. -------
  388. zeros : ndarray
  389. An array giving the roots of the spline.
  390. Notes
  391. -----
  392. Manipulating the tck-tuples directly is not recommended. In new code,
  393. prefer using the `BSpline` objects.
  394. See Also
  395. --------
  396. splprep, splrep, splint, spalde, splev
  397. bisplrep, bisplev
  398. BSpline
  399. References
  400. ----------
  401. .. [1] C. de Boor, "On calculating with b-splines", J. Approximation
  402. Theory, 6, p.50-62, 1972.
  403. .. [2] M. G. Cox, "The numerical evaluation of b-splines", J. Inst. Maths
  404. Applics, 10, p.134-149, 1972.
  405. .. [3] P. Dierckx, "Curve and surface fitting with splines", Monographs
  406. on Numerical Analysis, Oxford University Press, 1993.
  407. Examples
  408. --------
  409. For some data, this method may miss a root. This happens when one of
  410. the spline knots (which FITPACK places automatically) happens to
  411. coincide with the true root. A workaround is to convert to `PPoly`,
  412. which uses a different root-finding algorithm.
  413. For example,
  414. >>> x = [1.96, 1.97, 1.98, 1.99, 2.00, 2.01, 2.02, 2.03, 2.04, 2.05]
  415. >>> y = [-6.365470e-03, -4.790580e-03, -3.204320e-03, -1.607270e-03,
  416. ... 4.440892e-16, 1.616930e-03, 3.243000e-03, 4.877670e-03,
  417. ... 6.520430e-03, 8.170770e-03]
  418. >>> from scipy.interpolate import splrep, sproot, PPoly
  419. >>> tck = splrep(x, y, s=0)
  420. >>> sproot(tck)
  421. array([], dtype=float64)
  422. Converting to a PPoly object does find the roots at `x=2`:
  423. >>> ppoly = PPoly.from_spline(tck)
  424. >>> ppoly.roots(extrapolate=False)
  425. array([2.])
  426. Further examples are given :ref:`in the tutorial
  427. <tutorial-interpolate_splXXX>`.
  428. """
  429. if isinstance(tck, BSpline):
  430. if tck.c.ndim > 1:
  431. mesg = ("Calling sproot() with BSpline objects with c.ndim > 1 is "
  432. "not allowed.")
  433. raise ValueError(mesg)
  434. t, c, k = tck.tck
  435. # _impl.sproot expects the interpolation axis to be last, so roll it.
  436. # NB: This transpose is a no-op if c is 1D.
  437. sh = tuple(range(c.ndim))
  438. c = c.transpose(sh[1:] + (0,))
  439. return _impl.sproot((t, c, k), mest)
  440. else:
  441. return _impl.sproot(tck, mest)
  442. def spalde(x, tck):
  443. """
  444. Evaluate all derivatives of a B-spline.
  445. Given the knots and coefficients of a cubic B-spline compute all
  446. derivatives up to order k at a point (or set of points).
  447. Parameters
  448. ----------
  449. x : array_like
  450. A point or a set of points at which to evaluate the derivatives.
  451. Note that ``t(k) <= x <= t(n-k+1)`` must hold for each `x`.
  452. tck : tuple
  453. A tuple ``(t, c, k)``, containing the vector of knots, the B-spline
  454. coefficients, and the degree of the spline (see `splev`).
  455. Returns
  456. -------
  457. results : {ndarray, list of ndarrays}
  458. An array (or a list of arrays) containing all derivatives
  459. up to order k inclusive for each point `x`.
  460. See Also
  461. --------
  462. splprep, splrep, splint, sproot, splev, bisplrep, bisplev,
  463. BSpline
  464. References
  465. ----------
  466. .. [1] C. de Boor: On calculating with b-splines, J. Approximation Theory
  467. 6 (1972) 50-62.
  468. .. [2] M. G. Cox : The numerical evaluation of b-splines, J. Inst. Maths
  469. applics 10 (1972) 134-149.
  470. .. [3] P. Dierckx : Curve and surface fitting with splines, Monographs on
  471. Numerical Analysis, Oxford University Press, 1993.
  472. Examples
  473. --------
  474. Examples are given :ref:`in the tutorial <tutorial-interpolate_splXXX>`.
  475. """
  476. if isinstance(tck, BSpline):
  477. raise TypeError("spalde does not accept BSpline instances.")
  478. else:
  479. return _impl.spalde(x, tck)
  480. def insert(x, tck, m=1, per=0):
  481. """
  482. Insert knots into a B-spline.
  483. Given the knots and coefficients of a B-spline representation, create a
  484. new B-spline with a knot inserted `m` times at point `x`.
  485. This is a wrapper around the FORTRAN routine insert of FITPACK.
  486. Parameters
  487. ----------
  488. x (u) : array_like
  489. A 1-D point at which to insert a new knot(s). If `tck` was returned
  490. from ``splprep``, then the parameter values, u should be given.
  491. tck : a `BSpline` instance or a tuple
  492. If tuple, then it is expected to be a tuple (t,c,k) containing
  493. the vector of knots, the B-spline coefficients, and the degree of
  494. the spline.
  495. m : int, optional
  496. The number of times to insert the given knot (its multiplicity).
  497. Default is 1.
  498. per : int, optional
  499. If non-zero, the input spline is considered periodic.
  500. Returns
  501. -------
  502. BSpline instance or a tuple
  503. A new B-spline with knots t, coefficients c, and degree k.
  504. ``t(k+1) <= x <= t(n-k)``, where k is the degree of the spline.
  505. In case of a periodic spline (``per != 0``) there must be
  506. either at least k interior knots t(j) satisfying ``t(k+1)<t(j)<=x``
  507. or at least k interior knots t(j) satisfying ``x<=t(j)<t(n-k)``.
  508. A tuple is returned iff the input argument `tck` is a tuple, otherwise
  509. a BSpline object is constructed and returned.
  510. Notes
  511. -----
  512. Based on algorithms from [1]_ and [2]_.
  513. Manipulating the tck-tuples directly is not recommended. In new code,
  514. prefer using the `BSpline` objects.
  515. References
  516. ----------
  517. .. [1] W. Boehm, "Inserting new knots into b-spline curves.",
  518. Computer Aided Design, 12, p.199-201, 1980.
  519. .. [2] P. Dierckx, "Curve and surface fitting with splines, Monographs on
  520. Numerical Analysis", Oxford University Press, 1993.
  521. Examples
  522. --------
  523. You can insert knots into a B-spline.
  524. >>> from scipy.interpolate import splrep, insert
  525. >>> import numpy as np
  526. >>> x = np.linspace(0, 10, 5)
  527. >>> y = np.sin(x)
  528. >>> tck = splrep(x, y)
  529. >>> tck[0]
  530. array([ 0., 0., 0., 0., 5., 10., 10., 10., 10.])
  531. A knot is inserted:
  532. >>> tck_inserted = insert(3, tck)
  533. >>> tck_inserted[0]
  534. array([ 0., 0., 0., 0., 3., 5., 10., 10., 10., 10.])
  535. Some knots are inserted:
  536. >>> tck_inserted2 = insert(8, tck, m=3)
  537. >>> tck_inserted2[0]
  538. array([ 0., 0., 0., 0., 5., 8., 8., 8., 10., 10., 10., 10.])
  539. """
  540. if isinstance(tck, BSpline):
  541. t, c, k = tck.tck
  542. # FITPACK expects the interpolation axis to be last, so roll it over
  543. # NB: if c array is 1D, transposes are no-ops
  544. sh = tuple(range(c.ndim))
  545. c = c.transpose(sh[1:] + (0,))
  546. t_, c_, k_ = _impl.insert(x, (t, c, k), m, per)
  547. # and roll the last axis back
  548. c_ = np.asarray(c_)
  549. c_ = c_.transpose((sh[-1],) + sh[:-1])
  550. return BSpline(t_, c_, k_)
  551. else:
  552. return _impl.insert(x, tck, m, per)
  553. def splder(tck, n=1):
  554. """
  555. Compute the spline representation of the derivative of a given spline
  556. Parameters
  557. ----------
  558. tck : BSpline instance or a tuple of (t, c, k)
  559. Spline whose derivative to compute
  560. n : int, optional
  561. Order of derivative to evaluate. Default: 1
  562. Returns
  563. -------
  564. `BSpline` instance or tuple
  565. Spline of order k2=k-n representing the derivative
  566. of the input spline.
  567. A tuple is returned iff the input argument `tck` is a tuple, otherwise
  568. a BSpline object is constructed and returned.
  569. Notes
  570. -----
  571. .. versionadded:: 0.13.0
  572. See Also
  573. --------
  574. splantider, splev, spalde
  575. BSpline
  576. Examples
  577. --------
  578. This can be used for finding maxima of a curve:
  579. >>> from scipy.interpolate import splrep, splder, sproot
  580. >>> import numpy as np
  581. >>> x = np.linspace(0, 10, 70)
  582. >>> y = np.sin(x)
  583. >>> spl = splrep(x, y, k=4)
  584. Now, differentiate the spline and find the zeros of the
  585. derivative. (NB: `sproot` only works for order 3 splines, so we
  586. fit an order 4 spline):
  587. >>> dspl = splder(spl)
  588. >>> sproot(dspl) / np.pi
  589. array([ 0.50000001, 1.5 , 2.49999998])
  590. This agrees well with roots :math:`\\pi/2 + n\\pi` of
  591. :math:`\\cos(x) = \\sin'(x)`.
  592. """
  593. if isinstance(tck, BSpline):
  594. return tck.derivative(n)
  595. else:
  596. return _impl.splder(tck, n)
  597. def splantider(tck, n=1):
  598. """
  599. Compute the spline for the antiderivative (integral) of a given spline.
  600. Parameters
  601. ----------
  602. tck : BSpline instance or a tuple of (t, c, k)
  603. Spline whose antiderivative to compute
  604. n : int, optional
  605. Order of antiderivative to evaluate. Default: 1
  606. Returns
  607. -------
  608. BSpline instance or a tuple of (t2, c2, k2)
  609. Spline of order k2=k+n representing the antiderivative of the input
  610. spline.
  611. A tuple is returned iff the input argument `tck` is a tuple, otherwise
  612. a BSpline object is constructed and returned.
  613. See Also
  614. --------
  615. splder, splev, spalde
  616. BSpline
  617. Notes
  618. -----
  619. The `splder` function is the inverse operation of this function.
  620. Namely, ``splder(splantider(tck))`` is identical to `tck`, modulo
  621. rounding error.
  622. .. versionadded:: 0.13.0
  623. Examples
  624. --------
  625. >>> from scipy.interpolate import splrep, splder, splantider, splev
  626. >>> import numpy as np
  627. >>> x = np.linspace(0, np.pi/2, 70)
  628. >>> y = 1 / np.sqrt(1 - 0.8*np.sin(x)**2)
  629. >>> spl = splrep(x, y)
  630. The derivative is the inverse operation of the antiderivative,
  631. although some floating point error accumulates:
  632. >>> splev(1.7, spl), splev(1.7, splder(splantider(spl)))
  633. (array(2.1565429877197317), array(2.1565429877201865))
  634. Antiderivative can be used to evaluate definite integrals:
  635. >>> ispl = splantider(spl)
  636. >>> splev(np.pi/2, ispl) - splev(0, ispl)
  637. 2.2572053588768486
  638. This is indeed an approximation to the complete elliptic integral
  639. :math:`K(m) = \\int_0^{\\pi/2} [1 - m\\sin^2 x]^{-1/2} dx`:
  640. >>> from scipy.special import ellipk
  641. >>> ellipk(0.8)
  642. 2.2572053268208538
  643. """
  644. if isinstance(tck, BSpline):
  645. return tck.antiderivative(n)
  646. else:
  647. return _impl.splantider(tck, n)