_fitpack_py.py 32 KB

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