Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_interpolate.py: 17%
784 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1__all__ = ['interp1d', 'interp2d', 'lagrange', 'PPoly', 'BPoly', 'NdPPoly']
3from math import prod
5import numpy as np
6from numpy import (array, transpose, searchsorted, atleast_1d, atleast_2d,
7 ravel, poly1d, asarray, intp)
9import scipy.special as spec
10from scipy.special import comb
12from . import _fitpack_py
13from . import dfitpack
14from ._polyint import _Interpolator1D
15from . import _ppoly
16from .interpnd import _ndim_coords_from_arrays
17from ._bsplines import make_interp_spline, BSpline
20def lagrange(x, w):
21 r"""
22 Return a Lagrange interpolating polynomial.
24 Given two 1-D arrays `x` and `w,` returns the Lagrange interpolating
25 polynomial through the points ``(x, w)``.
27 Warning: This implementation is numerically unstable. Do not expect to
28 be able to use more than about 20 points even if they are chosen optimally.
30 Parameters
31 ----------
32 x : array_like
33 `x` represents the x-coordinates of a set of datapoints.
34 w : array_like
35 `w` represents the y-coordinates of a set of datapoints, i.e., f(`x`).
37 Returns
38 -------
39 lagrange : `numpy.poly1d` instance
40 The Lagrange interpolating polynomial.
42 Examples
43 --------
44 Interpolate :math:`f(x) = x^3` by 3 points.
46 >>> import numpy as np
47 >>> from scipy.interpolate import lagrange
48 >>> x = np.array([0, 1, 2])
49 >>> y = x**3
50 >>> poly = lagrange(x, y)
52 Since there are only 3 points, Lagrange polynomial has degree 2. Explicitly,
53 it is given by
55 .. math::
57 \begin{aligned}
58 L(x) &= 1\times \frac{x (x - 2)}{-1} + 8\times \frac{x (x-1)}{2} \\
59 &= x (-2 + 3x)
60 \end{aligned}
62 >>> from numpy.polynomial.polynomial import Polynomial
63 >>> Polynomial(poly.coef[::-1]).coef
64 array([ 0., -2., 3.])
66 >>> import matplotlib.pyplot as plt
67 >>> x_new = np.arange(0, 2.1, 0.1)
68 >>> plt.scatter(x, y, label='data')
69 >>> plt.plot(x_new, Polynomial(poly.coef[::-1])(x_new), label='Polynomial')
70 >>> plt.plot(x_new, 3*x_new**2 - 2*x_new + 0*x_new,
71 ... label=r"$3 x^2 - 2 x$", linestyle='-.')
72 >>> plt.legend()
73 >>> plt.show()
75 """
77 M = len(x)
78 p = poly1d(0.0)
79 for j in range(M):
80 pt = poly1d(w[j])
81 for k in range(M):
82 if k == j:
83 continue
84 fac = x[j]-x[k]
85 pt *= poly1d([1.0, -x[k]])/fac
86 p += pt
87 return p
90# !! Need to find argument for keeping initialize. If it isn't
91# !! found, get rid of it!
94dep_mesg = """\
95`interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy 1.13.0.
97For legacy code, nearly bug-for-bug compatible replacements are
98`RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
99scattered 2D data.
101In new code, for regular grids use `RegularGridInterpolator` instead.
102For scattered data, prefer `LinearNDInterpolator` or
103`CloughTocher2DInterpolator`.
105For more details see
106`https://scipy.github.io/devdocs/notebooks/interp_transition_guide.html`
107"""
109class interp2d:
110 """
111 interp2d(x, y, z, kind='linear', copy=True, bounds_error=False,
112 fill_value=None)
114 .. deprecated:: 1.10.0
116 `interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy
117 1.13.0.
119 For legacy code, nearly bug-for-bug compatible replacements are
120 `RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
121 scattered 2D data.
123 In new code, for regular grids use `RegularGridInterpolator` instead.
124 For scattered data, prefer `LinearNDInterpolator` or
125 `CloughTocher2DInterpolator`.
127 For more details see
128 `https://scipy.github.io/devdocs/notebooks/interp_transition_guide.html
129 <https://scipy.github.io/devdocs/notebooks/interp_transition_guide.html>`_
132 Interpolate over a 2-D grid.
134 `x`, `y` and `z` are arrays of values used to approximate some function
135 f: ``z = f(x, y)`` which returns a scalar value `z`. This class returns a
136 function whose call method uses spline interpolation to find the value
137 of new points.
139 If `x` and `y` represent a regular grid, consider using
140 `RectBivariateSpline`.
142 If `z` is a vector value, consider using `interpn`.
144 Note that calling `interp2d` with NaNs present in input values, or with
145 decreasing values in `x` an `y` results in undefined behaviour.
147 Methods
148 -------
149 __call__
151 Parameters
152 ----------
153 x, y : array_like
154 Arrays defining the data point coordinates.
155 The data point coordinates need to be sorted by increasing order.
157 If the points lie on a regular grid, `x` can specify the column
158 coordinates and `y` the row coordinates, for example::
160 >>> x = [0,1,2]; y = [0,3]; z = [[1,2,3], [4,5,6]]
162 Otherwise, `x` and `y` must specify the full coordinates for each
163 point, for example::
165 >>> x = [0,1,2,0,1,2]; y = [0,0,0,3,3,3]; z = [1,4,2,5,3,6]
167 If `x` and `y` are multidimensional, they are flattened before use.
168 z : array_like
169 The values of the function to interpolate at the data points. If
170 `z` is a multidimensional array, it is flattened before use assuming
171 Fortran-ordering (order='F'). The length of a flattened `z` array
172 is either len(`x`)*len(`y`) if `x` and `y` specify the column and
173 row coordinates or ``len(z) == len(x) == len(y)`` if `x` and `y`
174 specify coordinates for each point.
175 kind : {'linear', 'cubic', 'quintic'}, optional
176 The kind of spline interpolation to use. Default is 'linear'.
177 copy : bool, optional
178 If True, the class makes internal copies of x, y and z.
179 If False, references may be used. The default is to copy.
180 bounds_error : bool, optional
181 If True, when interpolated values are requested outside of the
182 domain of the input data (x,y), a ValueError is raised.
183 If False, then `fill_value` is used.
184 fill_value : number, optional
185 If provided, the value to use for points outside of the
186 interpolation domain. If omitted (None), values outside
187 the domain are extrapolated via nearest-neighbor extrapolation.
189 See Also
190 --------
191 RectBivariateSpline :
192 Much faster 2-D interpolation if your input data is on a grid
193 bisplrep, bisplev :
194 Spline interpolation based on FITPACK
195 BivariateSpline : a more recent wrapper of the FITPACK routines
196 interp1d : 1-D version of this function
197 RegularGridInterpolator : interpolation on a regular or rectilinear grid
198 in arbitrary dimensions.
199 interpn : Multidimensional interpolation on regular grids (wraps
200 `RegularGridInterpolator` and `RectBivariateSpline`).
202 Notes
203 -----
204 The minimum number of data points required along the interpolation
205 axis is ``(k+1)**2``, with k=1 for linear, k=3 for cubic and k=5 for
206 quintic interpolation.
208 The interpolator is constructed by `bisplrep`, with a smoothing factor
209 of 0. If more control over smoothing is needed, `bisplrep` should be
210 used directly.
212 The coordinates of the data points to interpolate `xnew` and `ynew`
213 have to be sorted by ascending order.
214 `interp2d` is legacy and is not
215 recommended for use in new code. New code should use
216 `RegularGridInterpolator` instead.
218 Examples
219 --------
220 Construct a 2-D grid and interpolate on it:
222 >>> import numpy as np
223 >>> from scipy import interpolate
224 >>> x = np.arange(-5.01, 5.01, 0.25)
225 >>> y = np.arange(-5.01, 5.01, 0.25)
226 >>> xx, yy = np.meshgrid(x, y)
227 >>> z = np.sin(xx**2+yy**2)
228 >>> f = interpolate.interp2d(x, y, z, kind='cubic')
230 Now use the obtained interpolation function and plot the result:
232 >>> import matplotlib.pyplot as plt
233 >>> xnew = np.arange(-5.01, 5.01, 1e-2)
234 >>> ynew = np.arange(-5.01, 5.01, 1e-2)
235 >>> znew = f(xnew, ynew)
236 >>> plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
237 >>> plt.show()
238 """
240 @np.deprecate(old_name='interp2d', message=dep_mesg)
241 def __init__(self, x, y, z, kind='linear', copy=True, bounds_error=False,
242 fill_value=None):
243 x = ravel(x)
244 y = ravel(y)
245 z = asarray(z)
247 rectangular_grid = (z.size == len(x) * len(y))
248 if rectangular_grid:
249 if z.ndim == 2:
250 if z.shape != (len(y), len(x)):
251 raise ValueError("When on a regular grid with x.size = m "
252 "and y.size = n, if z.ndim == 2, then z "
253 "must have shape (n, m)")
254 if not np.all(x[1:] >= x[:-1]):
255 j = np.argsort(x)
256 x = x[j]
257 z = z[:, j]
258 if not np.all(y[1:] >= y[:-1]):
259 j = np.argsort(y)
260 y = y[j]
261 z = z[j, :]
262 z = ravel(z.T)
263 else:
264 z = ravel(z)
265 if len(x) != len(y):
266 raise ValueError(
267 "x and y must have equal lengths for non rectangular grid")
268 if len(z) != len(x):
269 raise ValueError(
270 "Invalid length for input z for non rectangular grid")
272 interpolation_types = {'linear': 1, 'cubic': 3, 'quintic': 5}
273 try:
274 kx = ky = interpolation_types[kind]
275 except KeyError as e:
276 raise ValueError(
277 f"Unsupported interpolation type {repr(kind)}, must be "
278 f"either of {', '.join(map(repr, interpolation_types))}."
279 ) from e
281 if not rectangular_grid:
282 # TODO: surfit is really not meant for interpolation!
283 self.tck = _fitpack_py.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
284 else:
285 nx, tx, ny, ty, c, fp, ier = dfitpack.regrid_smth(
286 x, y, z, None, None, None, None,
287 kx=kx, ky=ky, s=0.0)
288 self.tck = (tx[:nx], ty[:ny], c[:(nx - kx - 1) * (ny - ky - 1)],
289 kx, ky)
291 self.bounds_error = bounds_error
292 self.fill_value = fill_value
293 self.x, self.y, self.z = (array(a, copy=copy) for a in (x, y, z))
295 self.x_min, self.x_max = np.amin(x), np.amax(x)
296 self.y_min, self.y_max = np.amin(y), np.amax(y)
298 @np.deprecate(old_name='interp2d', message=dep_mesg)
299 def __call__(self, x, y, dx=0, dy=0, assume_sorted=False):
300 """Interpolate the function.
302 Parameters
303 ----------
304 x : 1-D array
305 x-coordinates of the mesh on which to interpolate.
306 y : 1-D array
307 y-coordinates of the mesh on which to interpolate.
308 dx : int >= 0, < kx
309 Order of partial derivatives in x.
310 dy : int >= 0, < ky
311 Order of partial derivatives in y.
312 assume_sorted : bool, optional
313 If False, values of `x` and `y` can be in any order and they are
314 sorted first.
315 If True, `x` and `y` have to be arrays of monotonically
316 increasing values.
318 Returns
319 -------
320 z : 2-D array with shape (len(y), len(x))
321 The interpolated values.
322 """
324 x = atleast_1d(x)
325 y = atleast_1d(y)
327 if x.ndim != 1 or y.ndim != 1:
328 raise ValueError("x and y should both be 1-D arrays")
330 if not assume_sorted:
331 x = np.sort(x, kind="mergesort")
332 y = np.sort(y, kind="mergesort")
334 if self.bounds_error or self.fill_value is not None:
335 out_of_bounds_x = (x < self.x_min) | (x > self.x_max)
336 out_of_bounds_y = (y < self.y_min) | (y > self.y_max)
338 any_out_of_bounds_x = np.any(out_of_bounds_x)
339 any_out_of_bounds_y = np.any(out_of_bounds_y)
341 if self.bounds_error and (any_out_of_bounds_x or any_out_of_bounds_y):
342 raise ValueError("Values out of range; x must be in %r, y in %r"
343 % ((self.x_min, self.x_max),
344 (self.y_min, self.y_max)))
346 z = _fitpack_py.bisplev(x, y, self.tck, dx, dy)
347 z = atleast_2d(z)
348 z = transpose(z)
350 if self.fill_value is not None:
351 if any_out_of_bounds_x:
352 z[:, out_of_bounds_x] = self.fill_value
353 if any_out_of_bounds_y:
354 z[out_of_bounds_y, :] = self.fill_value
356 if len(z) == 1:
357 z = z[0]
358 return array(z)
361def _check_broadcast_up_to(arr_from, shape_to, name):
362 """Helper to check that arr_from broadcasts up to shape_to"""
363 shape_from = arr_from.shape
364 if len(shape_to) >= len(shape_from):
365 for t, f in zip(shape_to[::-1], shape_from[::-1]):
366 if f != 1 and f != t:
367 break
368 else: # all checks pass, do the upcasting that we need later
369 if arr_from.size != 1 and arr_from.shape != shape_to:
370 arr_from = np.ones(shape_to, arr_from.dtype) * arr_from
371 return arr_from.ravel()
372 # at least one check failed
373 raise ValueError('%s argument must be able to broadcast up '
374 'to shape %s but had shape %s'
375 % (name, shape_to, shape_from))
378def _do_extrapolate(fill_value):
379 """Helper to check if fill_value == "extrapolate" without warnings"""
380 return (isinstance(fill_value, str) and
381 fill_value == 'extrapolate')
384class interp1d(_Interpolator1D):
385 """
386 Interpolate a 1-D function.
388 .. legacy:: class
390 `x` and `y` are arrays of values used to approximate some function f:
391 ``y = f(x)``. This class returns a function whose call method uses
392 interpolation to find the value of new points.
394 Parameters
395 ----------
396 x : (npoints, ) array_like
397 A 1-D array of real values.
398 y : (..., npoints, ...) array_like
399 A N-D array of real values. The length of `y` along the interpolation
400 axis must be equal to the length of `x`. Use the ``axis`` parameter
401 to select correct axis. Unlike other interpolators, the default
402 interpolation axis is the last axis of `y`.
403 kind : str or int, optional
404 Specifies the kind of interpolation as a string or as an integer
405 specifying the order of the spline interpolator to use.
406 The string has to be one of 'linear', 'nearest', 'nearest-up', 'zero',
407 'slinear', 'quadratic', 'cubic', 'previous', or 'next'. 'zero',
408 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation of
409 zeroth, first, second or third order; 'previous' and 'next' simply
410 return the previous or next value of the point; 'nearest-up' and
411 'nearest' differ when interpolating half-integers (e.g. 0.5, 1.5)
412 in that 'nearest-up' rounds up and 'nearest' rounds down. Default
413 is 'linear'.
414 axis : int, optional
415 Axis in the ``y`` array corresponding to the x-coordinate values. Unlike
416 other interpolators, defaults to ``axis=-1``.
417 copy : bool, optional
418 If True, the class makes internal copies of x and y.
419 If False, references to `x` and `y` are used. The default is to copy.
420 bounds_error : bool, optional
421 If True, a ValueError is raised any time interpolation is attempted on
422 a value outside of the range of x (where extrapolation is
423 necessary). If False, out of bounds values are assigned `fill_value`.
424 By default, an error is raised unless ``fill_value="extrapolate"``.
425 fill_value : array-like or (array-like, array_like) or "extrapolate", optional
426 - if a ndarray (or float), this value will be used to fill in for
427 requested points outside of the data range. If not provided, then
428 the default is NaN. The array-like must broadcast properly to the
429 dimensions of the non-interpolation axes.
430 - If a two-element tuple, then the first element is used as a
431 fill value for ``x_new < x[0]`` and the second element is used for
432 ``x_new > x[-1]``. Anything that is not a 2-element tuple (e.g.,
433 list or ndarray, regardless of shape) is taken to be a single
434 array-like argument meant to be used for both bounds as
435 ``below, above = fill_value, fill_value``. Using a two-element tuple
436 or ndarray requires ``bounds_error=False``.
438 .. versionadded:: 0.17.0
439 - If "extrapolate", then points outside the data range will be
440 extrapolated.
442 .. versionadded:: 0.17.0
443 assume_sorted : bool, optional
444 If False, values of `x` can be in any order and they are sorted first.
445 If True, `x` has to be an array of monotonically increasing values.
447 Attributes
448 ----------
449 fill_value
451 Methods
452 -------
453 __call__
455 See Also
456 --------
457 splrep, splev
458 Spline interpolation/smoothing based on FITPACK.
459 UnivariateSpline : An object-oriented wrapper of the FITPACK routines.
460 interp2d : 2-D interpolation
462 Notes
463 -----
464 Calling `interp1d` with NaNs present in input values results in
465 undefined behaviour.
467 Input values `x` and `y` must be convertible to `float` values like
468 `int` or `float`.
470 If the values in `x` are not unique, the resulting behavior is
471 undefined and specific to the choice of `kind`, i.e., changing
472 `kind` will change the behavior for duplicates.
475 Examples
476 --------
477 >>> import numpy as np
478 >>> import matplotlib.pyplot as plt
479 >>> from scipy import interpolate
480 >>> x = np.arange(0, 10)
481 >>> y = np.exp(-x/3.0)
482 >>> f = interpolate.interp1d(x, y)
484 >>> xnew = np.arange(0, 9, 0.1)
485 >>> ynew = f(xnew) # use interpolation function returned by `interp1d`
486 >>> plt.plot(x, y, 'o', xnew, ynew, '-')
487 >>> plt.show()
488 """
490 def __init__(self, x, y, kind='linear', axis=-1,
491 copy=True, bounds_error=None, fill_value=np.nan,
492 assume_sorted=False):
493 """ Initialize a 1-D linear interpolation class."""
494 _Interpolator1D.__init__(self, x, y, axis=axis)
496 self.bounds_error = bounds_error # used by fill_value setter
497 self.copy = copy
499 if kind in ['zero', 'slinear', 'quadratic', 'cubic']:
500 order = {'zero': 0, 'slinear': 1,
501 'quadratic': 2, 'cubic': 3}[kind]
502 kind = 'spline'
503 elif isinstance(kind, int):
504 order = kind
505 kind = 'spline'
506 elif kind not in ('linear', 'nearest', 'nearest-up', 'previous',
507 'next'):
508 raise NotImplementedError("%s is unsupported: Use fitpack "
509 "routines for other types." % kind)
510 x = array(x, copy=self.copy)
511 y = array(y, copy=self.copy)
513 if not assume_sorted:
514 ind = np.argsort(x, kind="mergesort")
515 x = x[ind]
516 y = np.take(y, ind, axis=axis)
518 if x.ndim != 1:
519 raise ValueError("the x array must have exactly one dimension.")
520 if y.ndim == 0:
521 raise ValueError("the y array must have at least one dimension.")
523 # Force-cast y to a floating-point type, if it's not yet one
524 if not issubclass(y.dtype.type, np.inexact):
525 y = y.astype(np.float_)
527 # Backward compatibility
528 self.axis = axis % y.ndim
530 # Interpolation goes internally along the first axis
531 self.y = y
532 self._y = self._reshape_yi(self.y)
533 self.x = x
534 del y, x # clean up namespace to prevent misuse; use attributes
535 self._kind = kind
537 # Adjust to interpolation kind; store reference to *unbound*
538 # interpolation methods, in order to avoid circular references to self
539 # stored in the bound instance methods, and therefore delayed garbage
540 # collection. See: https://docs.python.org/reference/datamodel.html
541 if kind in ('linear', 'nearest', 'nearest-up', 'previous', 'next'):
542 # Make a "view" of the y array that is rotated to the interpolation
543 # axis.
544 minval = 1
545 if kind == 'nearest':
546 # Do division before addition to prevent possible integer
547 # overflow
548 self._side = 'left'
549 self.x_bds = self.x / 2.0
550 self.x_bds = self.x_bds[1:] + self.x_bds[:-1]
552 self._call = self.__class__._call_nearest
553 elif kind == 'nearest-up':
554 # Do division before addition to prevent possible integer
555 # overflow
556 self._side = 'right'
557 self.x_bds = self.x / 2.0
558 self.x_bds = self.x_bds[1:] + self.x_bds[:-1]
560 self._call = self.__class__._call_nearest
561 elif kind == 'previous':
562 # Side for np.searchsorted and index for clipping
563 self._side = 'left'
564 self._ind = 0
565 # Move x by one floating point value to the left
566 self._x_shift = np.nextafter(self.x, -np.inf)
567 self._call = self.__class__._call_previousnext
568 if _do_extrapolate(fill_value):
569 self._check_and_update_bounds_error_for_extrapolation()
570 # assume y is sorted by x ascending order here.
571 fill_value = (np.nan, np.take(self.y, -1, axis))
572 elif kind == 'next':
573 self._side = 'right'
574 self._ind = 1
575 # Move x by one floating point value to the right
576 self._x_shift = np.nextafter(self.x, np.inf)
577 self._call = self.__class__._call_previousnext
578 if _do_extrapolate(fill_value):
579 self._check_and_update_bounds_error_for_extrapolation()
580 # assume y is sorted by x ascending order here.
581 fill_value = (np.take(self.y, 0, axis), np.nan)
582 else:
583 # Check if we can delegate to numpy.interp (2x-10x faster).
584 np_types = (np.float_, np.int_)
585 cond = self.x.dtype in np_types and self.y.dtype in np_types
586 cond = cond and self.y.ndim == 1
587 cond = cond and not _do_extrapolate(fill_value)
589 if cond:
590 self._call = self.__class__._call_linear_np
591 else:
592 self._call = self.__class__._call_linear
593 else:
594 minval = order + 1
596 rewrite_nan = False
597 xx, yy = self.x, self._y
598 if order > 1:
599 # Quadratic or cubic spline. If input contains even a single
600 # nan, then the output is all nans. We cannot just feed data
601 # with nans to make_interp_spline because it calls LAPACK.
602 # So, we make up a bogus x and y with no nans and use it
603 # to get the correct shape of the output, which we then fill
604 # with nans.
605 # For slinear or zero order spline, we just pass nans through.
606 mask = np.isnan(self.x)
607 if mask.any():
608 sx = self.x[~mask]
609 if sx.size == 0:
610 raise ValueError("`x` array is all-nan")
611 xx = np.linspace(np.nanmin(self.x),
612 np.nanmax(self.x),
613 len(self.x))
614 rewrite_nan = True
615 if np.isnan(self._y).any():
616 yy = np.ones_like(self._y)
617 rewrite_nan = True
619 self._spline = make_interp_spline(xx, yy, k=order,
620 check_finite=False)
621 if rewrite_nan:
622 self._call = self.__class__._call_nan_spline
623 else:
624 self._call = self.__class__._call_spline
626 if len(self.x) < minval:
627 raise ValueError("x and y arrays must have at "
628 "least %d entries" % minval)
630 self.fill_value = fill_value # calls the setter, can modify bounds_err
632 @property
633 def fill_value(self):
634 """The fill value."""
635 # backwards compat: mimic a public attribute
636 return self._fill_value_orig
638 @fill_value.setter
639 def fill_value(self, fill_value):
640 # extrapolation only works for nearest neighbor and linear methods
641 if _do_extrapolate(fill_value):
642 self._check_and_update_bounds_error_for_extrapolation()
643 self._extrapolate = True
644 else:
645 broadcast_shape = (self.y.shape[:self.axis] +
646 self.y.shape[self.axis + 1:])
647 if len(broadcast_shape) == 0:
648 broadcast_shape = (1,)
649 # it's either a pair (_below_range, _above_range) or a single value
650 # for both above and below range
651 if isinstance(fill_value, tuple) and len(fill_value) == 2:
652 below_above = [np.asarray(fill_value[0]),
653 np.asarray(fill_value[1])]
654 names = ('fill_value (below)', 'fill_value (above)')
655 for ii in range(2):
656 below_above[ii] = _check_broadcast_up_to(
657 below_above[ii], broadcast_shape, names[ii])
658 else:
659 fill_value = np.asarray(fill_value)
660 below_above = [_check_broadcast_up_to(
661 fill_value, broadcast_shape, 'fill_value')] * 2
662 self._fill_value_below, self._fill_value_above = below_above
663 self._extrapolate = False
664 if self.bounds_error is None:
665 self.bounds_error = True
666 # backwards compat: fill_value was a public attr; make it writeable
667 self._fill_value_orig = fill_value
669 def _check_and_update_bounds_error_for_extrapolation(self):
670 if self.bounds_error:
671 raise ValueError("Cannot extrapolate and raise "
672 "at the same time.")
673 self.bounds_error = False
675 def _call_linear_np(self, x_new):
676 # Note that out-of-bounds values are taken care of in self._evaluate
677 return np.interp(x_new, self.x, self.y)
679 def _call_linear(self, x_new):
680 # 2. Find where in the original data, the values to interpolate
681 # would be inserted.
682 # Note: If x_new[n] == x[m], then m is returned by searchsorted.
683 x_new_indices = searchsorted(self.x, x_new)
685 # 3. Clip x_new_indices so that they are within the range of
686 # self.x indices and at least 1. Removes mis-interpolation
687 # of x_new[n] = x[0]
688 x_new_indices = x_new_indices.clip(1, len(self.x)-1).astype(int)
690 # 4. Calculate the slope of regions that each x_new value falls in.
691 lo = x_new_indices - 1
692 hi = x_new_indices
694 x_lo = self.x[lo]
695 x_hi = self.x[hi]
696 y_lo = self._y[lo]
697 y_hi = self._y[hi]
699 # Note that the following two expressions rely on the specifics of the
700 # broadcasting semantics.
701 slope = (y_hi - y_lo) / (x_hi - x_lo)[:, None]
703 # 5. Calculate the actual value for each entry in x_new.
704 y_new = slope*(x_new - x_lo)[:, None] + y_lo
706 return y_new
708 def _call_nearest(self, x_new):
709 """ Find nearest neighbor interpolated y_new = f(x_new)."""
711 # 2. Find where in the averaged data the values to interpolate
712 # would be inserted.
713 # Note: use side='left' (right) to searchsorted() to define the
714 # halfway point to be nearest to the left (right) neighbor
715 x_new_indices = searchsorted(self.x_bds, x_new, side=self._side)
717 # 3. Clip x_new_indices so that they are within the range of x indices.
718 x_new_indices = x_new_indices.clip(0, len(self.x)-1).astype(intp)
720 # 4. Calculate the actual value for each entry in x_new.
721 y_new = self._y[x_new_indices]
723 return y_new
725 def _call_previousnext(self, x_new):
726 """Use previous/next neighbor of x_new, y_new = f(x_new)."""
728 # 1. Get index of left/right value
729 x_new_indices = searchsorted(self._x_shift, x_new, side=self._side)
731 # 2. Clip x_new_indices so that they are within the range of x indices.
732 x_new_indices = x_new_indices.clip(1-self._ind,
733 len(self.x)-self._ind).astype(intp)
735 # 3. Calculate the actual value for each entry in x_new.
736 y_new = self._y[x_new_indices+self._ind-1]
738 return y_new
740 def _call_spline(self, x_new):
741 return self._spline(x_new)
743 def _call_nan_spline(self, x_new):
744 out = self._spline(x_new)
745 out[...] = np.nan
746 return out
748 def _evaluate(self, x_new):
749 # 1. Handle values in x_new that are outside of x. Throw error,
750 # or return a list of mask array indicating the outofbounds values.
751 # The behavior is set by the bounds_error variable.
752 x_new = asarray(x_new)
753 y_new = self._call(self, x_new)
754 if not self._extrapolate:
755 below_bounds, above_bounds = self._check_bounds(x_new)
756 if len(y_new) > 0:
757 # Note fill_value must be broadcast up to the proper size
758 # and flattened to work here
759 y_new[below_bounds] = self._fill_value_below
760 y_new[above_bounds] = self._fill_value_above
761 return y_new
763 def _check_bounds(self, x_new):
764 """Check the inputs for being in the bounds of the interpolated data.
766 Parameters
767 ----------
768 x_new : array
770 Returns
771 -------
772 out_of_bounds : bool array
773 The mask on x_new of values that are out of the bounds.
774 """
776 # If self.bounds_error is True, we raise an error if any x_new values
777 # fall outside the range of x. Otherwise, we return an array indicating
778 # which values are outside the boundary region.
779 below_bounds = x_new < self.x[0]
780 above_bounds = x_new > self.x[-1]
782 if self.bounds_error and below_bounds.any():
783 below_bounds_value = x_new[np.argmax(below_bounds)]
784 raise ValueError("A value ({}) in x_new is below "
785 "the interpolation range's minimum value ({})."
786 .format(below_bounds_value, self.x[0]))
787 if self.bounds_error and above_bounds.any():
788 above_bounds_value = x_new[np.argmax(above_bounds)]
789 raise ValueError("A value ({}) in x_new is above "
790 "the interpolation range's maximum value ({})."
791 .format(above_bounds_value, self.x[-1]))
793 # !! Should we emit a warning if some values are out of bounds?
794 # !! matlab does not.
795 return below_bounds, above_bounds
798class _PPolyBase:
799 """Base class for piecewise polynomials."""
800 __slots__ = ('c', 'x', 'extrapolate', 'axis')
802 def __init__(self, c, x, extrapolate=None, axis=0):
803 self.c = np.asarray(c)
804 self.x = np.ascontiguousarray(x, dtype=np.float64)
806 if extrapolate is None:
807 extrapolate = True
808 elif extrapolate != 'periodic':
809 extrapolate = bool(extrapolate)
810 self.extrapolate = extrapolate
812 if self.c.ndim < 2:
813 raise ValueError("Coefficients array must be at least "
814 "2-dimensional.")
816 if not (0 <= axis < self.c.ndim - 1):
817 raise ValueError("axis=%s must be between 0 and %s" %
818 (axis, self.c.ndim-1))
820 self.axis = axis
821 if axis != 0:
822 # move the interpolation axis to be the first one in self.c
823 # More specifically, the target shape for self.c is (k, m, ...),
824 # and axis !=0 means that we have c.shape (..., k, m, ...)
825 # ^
826 # axis
827 # So we roll two of them.
828 self.c = np.moveaxis(self.c, axis+1, 0)
829 self.c = np.moveaxis(self.c, axis+1, 0)
831 if self.x.ndim != 1:
832 raise ValueError("x must be 1-dimensional")
833 if self.x.size < 2:
834 raise ValueError("at least 2 breakpoints are needed")
835 if self.c.ndim < 2:
836 raise ValueError("c must have at least 2 dimensions")
837 if self.c.shape[0] == 0:
838 raise ValueError("polynomial must be at least of order 0")
839 if self.c.shape[1] != self.x.size-1:
840 raise ValueError("number of coefficients != len(x)-1")
841 dx = np.diff(self.x)
842 if not (np.all(dx >= 0) or np.all(dx <= 0)):
843 raise ValueError("`x` must be strictly increasing or decreasing.")
845 dtype = self._get_dtype(self.c.dtype)
846 self.c = np.ascontiguousarray(self.c, dtype=dtype)
848 def _get_dtype(self, dtype):
849 if np.issubdtype(dtype, np.complexfloating) \
850 or np.issubdtype(self.c.dtype, np.complexfloating):
851 return np.complex_
852 else:
853 return np.float_
855 @classmethod
856 def construct_fast(cls, c, x, extrapolate=None, axis=0):
857 """
858 Construct the piecewise polynomial without making checks.
860 Takes the same parameters as the constructor. Input arguments
861 ``c`` and ``x`` must be arrays of the correct shape and type. The
862 ``c`` array can only be of dtypes float and complex, and ``x``
863 array must have dtype float.
864 """
865 self = object.__new__(cls)
866 self.c = c
867 self.x = x
868 self.axis = axis
869 if extrapolate is None:
870 extrapolate = True
871 self.extrapolate = extrapolate
872 return self
874 def _ensure_c_contiguous(self):
875 """
876 c and x may be modified by the user. The Cython code expects
877 that they are C contiguous.
878 """
879 if not self.x.flags.c_contiguous:
880 self.x = self.x.copy()
881 if not self.c.flags.c_contiguous:
882 self.c = self.c.copy()
884 def extend(self, c, x):
885 """
886 Add additional breakpoints and coefficients to the polynomial.
888 Parameters
889 ----------
890 c : ndarray, size (k, m, ...)
891 Additional coefficients for polynomials in intervals. Note that
892 the first additional interval will be formed using one of the
893 ``self.x`` end points.
894 x : ndarray, size (m,)
895 Additional breakpoints. Must be sorted in the same order as
896 ``self.x`` and either to the right or to the left of the current
897 breakpoints.
898 """
900 c = np.asarray(c)
901 x = np.asarray(x)
903 if c.ndim < 2:
904 raise ValueError("invalid dimensions for c")
905 if x.ndim != 1:
906 raise ValueError("invalid dimensions for x")
907 if x.shape[0] != c.shape[1]:
908 raise ValueError("Shapes of x {} and c {} are incompatible"
909 .format(x.shape, c.shape))
910 if c.shape[2:] != self.c.shape[2:] or c.ndim != self.c.ndim:
911 raise ValueError("Shapes of c {} and self.c {} are incompatible"
912 .format(c.shape, self.c.shape))
914 if c.size == 0:
915 return
917 dx = np.diff(x)
918 if not (np.all(dx >= 0) or np.all(dx <= 0)):
919 raise ValueError("`x` is not sorted.")
921 if self.x[-1] >= self.x[0]:
922 if not x[-1] >= x[0]:
923 raise ValueError("`x` is in the different order "
924 "than `self.x`.")
926 if x[0] >= self.x[-1]:
927 action = 'append'
928 elif x[-1] <= self.x[0]:
929 action = 'prepend'
930 else:
931 raise ValueError("`x` is neither on the left or on the right "
932 "from `self.x`.")
933 else:
934 if not x[-1] <= x[0]:
935 raise ValueError("`x` is in the different order "
936 "than `self.x`.")
938 if x[0] <= self.x[-1]:
939 action = 'append'
940 elif x[-1] >= self.x[0]:
941 action = 'prepend'
942 else:
943 raise ValueError("`x` is neither on the left or on the right "
944 "from `self.x`.")
946 dtype = self._get_dtype(c.dtype)
948 k2 = max(c.shape[0], self.c.shape[0])
949 c2 = np.zeros((k2, self.c.shape[1] + c.shape[1]) + self.c.shape[2:],
950 dtype=dtype)
952 if action == 'append':
953 c2[k2-self.c.shape[0]:, :self.c.shape[1]] = self.c
954 c2[k2-c.shape[0]:, self.c.shape[1]:] = c
955 self.x = np.r_[self.x, x]
956 elif action == 'prepend':
957 c2[k2-self.c.shape[0]:, :c.shape[1]] = c
958 c2[k2-c.shape[0]:, c.shape[1]:] = self.c
959 self.x = np.r_[x, self.x]
961 self.c = c2
963 def __call__(self, x, nu=0, extrapolate=None):
964 """
965 Evaluate the piecewise polynomial or its derivative.
967 Parameters
968 ----------
969 x : array_like
970 Points to evaluate the interpolant at.
971 nu : int, optional
972 Order of derivative to evaluate. Must be non-negative.
973 extrapolate : {bool, 'periodic', None}, optional
974 If bool, determines whether to extrapolate to out-of-bounds points
975 based on first and last intervals, or to return NaNs.
976 If 'periodic', periodic extrapolation is used.
977 If None (default), use `self.extrapolate`.
979 Returns
980 -------
981 y : array_like
982 Interpolated values. Shape is determined by replacing
983 the interpolation axis in the original array with the shape of x.
985 Notes
986 -----
987 Derivatives are evaluated piecewise for each polynomial
988 segment, even if the polynomial is not differentiable at the
989 breakpoints. The polynomial intervals are considered half-open,
990 ``[a, b)``, except for the last interval which is closed
991 ``[a, b]``.
992 """
993 if extrapolate is None:
994 extrapolate = self.extrapolate
995 x = np.asarray(x)
996 x_shape, x_ndim = x.shape, x.ndim
997 x = np.ascontiguousarray(x.ravel(), dtype=np.float_)
999 # With periodic extrapolation we map x to the segment
1000 # [self.x[0], self.x[-1]].
1001 if extrapolate == 'periodic':
1002 x = self.x[0] + (x - self.x[0]) % (self.x[-1] - self.x[0])
1003 extrapolate = False
1005 out = np.empty((len(x), prod(self.c.shape[2:])), dtype=self.c.dtype)
1006 self._ensure_c_contiguous()
1007 self._evaluate(x, nu, extrapolate, out)
1008 out = out.reshape(x_shape + self.c.shape[2:])
1009 if self.axis != 0:
1010 # transpose to move the calculated values to the interpolation axis
1011 l = list(range(out.ndim))
1012 l = l[x_ndim:x_ndim+self.axis] + l[:x_ndim] + l[x_ndim+self.axis:]
1013 out = out.transpose(l)
1014 return out
1017class PPoly(_PPolyBase):
1018 """
1019 Piecewise polynomial in terms of coefficients and breakpoints
1021 The polynomial between ``x[i]`` and ``x[i + 1]`` is written in the
1022 local power basis::
1024 S = sum(c[m, i] * (xp - x[i])**(k-m) for m in range(k+1))
1026 where ``k`` is the degree of the polynomial.
1028 Parameters
1029 ----------
1030 c : ndarray, shape (k, m, ...)
1031 Polynomial coefficients, order `k` and `m` intervals.
1032 x : ndarray, shape (m+1,)
1033 Polynomial breakpoints. Must be sorted in either increasing or
1034 decreasing order.
1035 extrapolate : bool or 'periodic', optional
1036 If bool, determines whether to extrapolate to out-of-bounds points
1037 based on first and last intervals, or to return NaNs. If 'periodic',
1038 periodic extrapolation is used. Default is True.
1039 axis : int, optional
1040 Interpolation axis. Default is zero.
1042 Attributes
1043 ----------
1044 x : ndarray
1045 Breakpoints.
1046 c : ndarray
1047 Coefficients of the polynomials. They are reshaped
1048 to a 3-D array with the last dimension representing
1049 the trailing dimensions of the original coefficient array.
1050 axis : int
1051 Interpolation axis.
1053 Methods
1054 -------
1055 __call__
1056 derivative
1057 antiderivative
1058 integrate
1059 solve
1060 roots
1061 extend
1062 from_spline
1063 from_bernstein_basis
1064 construct_fast
1066 See also
1067 --------
1068 BPoly : piecewise polynomials in the Bernstein basis
1070 Notes
1071 -----
1072 High-order polynomials in the power basis can be numerically
1073 unstable. Precision problems can start to appear for orders
1074 larger than 20-30.
1075 """
1077 def _evaluate(self, x, nu, extrapolate, out):
1078 _ppoly.evaluate(self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1079 self.x, x, nu, bool(extrapolate), out)
1081 def derivative(self, nu=1):
1082 """
1083 Construct a new piecewise polynomial representing the derivative.
1085 Parameters
1086 ----------
1087 nu : int, optional
1088 Order of derivative to evaluate. Default is 1, i.e., compute the
1089 first derivative. If negative, the antiderivative is returned.
1091 Returns
1092 -------
1093 pp : PPoly
1094 Piecewise polynomial of order k2 = k - n representing the derivative
1095 of this polynomial.
1097 Notes
1098 -----
1099 Derivatives are evaluated piecewise for each polynomial
1100 segment, even if the polynomial is not differentiable at the
1101 breakpoints. The polynomial intervals are considered half-open,
1102 ``[a, b)``, except for the last interval which is closed
1103 ``[a, b]``.
1104 """
1105 if nu < 0:
1106 return self.antiderivative(-nu)
1108 # reduce order
1109 if nu == 0:
1110 c2 = self.c.copy()
1111 else:
1112 c2 = self.c[:-nu, :].copy()
1114 if c2.shape[0] == 0:
1115 # derivative of order 0 is zero
1116 c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)
1118 # multiply by the correct rising factorials
1119 factor = spec.poch(np.arange(c2.shape[0], 0, -1), nu)
1120 c2 *= factor[(slice(None),) + (None,)*(c2.ndim-1)]
1122 # construct a compatible polynomial
1123 return self.construct_fast(c2, self.x, self.extrapolate, self.axis)
1125 def antiderivative(self, nu=1):
1126 """
1127 Construct a new piecewise polynomial representing the antiderivative.
1129 Antiderivative is also the indefinite integral of the function,
1130 and derivative is its inverse operation.
1132 Parameters
1133 ----------
1134 nu : int, optional
1135 Order of antiderivative to evaluate. Default is 1, i.e., compute
1136 the first integral. If negative, the derivative is returned.
1138 Returns
1139 -------
1140 pp : PPoly
1141 Piecewise polynomial of order k2 = k + n representing
1142 the antiderivative of this polynomial.
1144 Notes
1145 -----
1146 The antiderivative returned by this function is continuous and
1147 continuously differentiable to order n-1, up to floating point
1148 rounding error.
1150 If antiderivative is computed and ``self.extrapolate='periodic'``,
1151 it will be set to False for the returned instance. This is done because
1152 the antiderivative is no longer periodic and its correct evaluation
1153 outside of the initially given x interval is difficult.
1154 """
1155 if nu <= 0:
1156 return self.derivative(-nu)
1158 c = np.zeros((self.c.shape[0] + nu, self.c.shape[1]) + self.c.shape[2:],
1159 dtype=self.c.dtype)
1160 c[:-nu] = self.c
1162 # divide by the correct rising factorials
1163 factor = spec.poch(np.arange(self.c.shape[0], 0, -1), nu)
1164 c[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]
1166 # fix continuity of added degrees of freedom
1167 self._ensure_c_contiguous()
1168 _ppoly.fix_continuity(c.reshape(c.shape[0], c.shape[1], -1),
1169 self.x, nu - 1)
1171 if self.extrapolate == 'periodic':
1172 extrapolate = False
1173 else:
1174 extrapolate = self.extrapolate
1176 # construct a compatible polynomial
1177 return self.construct_fast(c, self.x, extrapolate, self.axis)
1179 def integrate(self, a, b, extrapolate=None):
1180 """
1181 Compute a definite integral over a piecewise polynomial.
1183 Parameters
1184 ----------
1185 a : float
1186 Lower integration bound
1187 b : float
1188 Upper integration bound
1189 extrapolate : {bool, 'periodic', None}, optional
1190 If bool, determines whether to extrapolate to out-of-bounds points
1191 based on first and last intervals, or to return NaNs.
1192 If 'periodic', periodic extrapolation is used.
1193 If None (default), use `self.extrapolate`.
1195 Returns
1196 -------
1197 ig : array_like
1198 Definite integral of the piecewise polynomial over [a, b]
1199 """
1200 if extrapolate is None:
1201 extrapolate = self.extrapolate
1203 # Swap integration bounds if needed
1204 sign = 1
1205 if b < a:
1206 a, b = b, a
1207 sign = -1
1209 range_int = np.empty((prod(self.c.shape[2:]),), dtype=self.c.dtype)
1210 self._ensure_c_contiguous()
1212 # Compute the integral.
1213 if extrapolate == 'periodic':
1214 # Split the integral into the part over period (can be several
1215 # of them) and the remaining part.
1217 xs, xe = self.x[0], self.x[-1]
1218 period = xe - xs
1219 interval = b - a
1220 n_periods, left = divmod(interval, period)
1222 if n_periods > 0:
1223 _ppoly.integrate(
1224 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1225 self.x, xs, xe, False, out=range_int)
1226 range_int *= n_periods
1227 else:
1228 range_int.fill(0)
1230 # Map a to [xs, xe], b is always a + left.
1231 a = xs + (a - xs) % period
1232 b = a + left
1234 # If b <= xe then we need to integrate over [a, b], otherwise
1235 # over [a, xe] and from xs to what is remained.
1236 remainder_int = np.empty_like(range_int)
1237 if b <= xe:
1238 _ppoly.integrate(
1239 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1240 self.x, a, b, False, out=remainder_int)
1241 range_int += remainder_int
1242 else:
1243 _ppoly.integrate(
1244 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1245 self.x, a, xe, False, out=remainder_int)
1246 range_int += remainder_int
1248 _ppoly.integrate(
1249 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1250 self.x, xs, xs + left + a - xe, False, out=remainder_int)
1251 range_int += remainder_int
1252 else:
1253 _ppoly.integrate(
1254 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1255 self.x, a, b, bool(extrapolate), out=range_int)
1257 # Return
1258 range_int *= sign
1259 return range_int.reshape(self.c.shape[2:])
1261 def solve(self, y=0., discontinuity=True, extrapolate=None):
1262 """
1263 Find real solutions of the equation ``pp(x) == y``.
1265 Parameters
1266 ----------
1267 y : float, optional
1268 Right-hand side. Default is zero.
1269 discontinuity : bool, optional
1270 Whether to report sign changes across discontinuities at
1271 breakpoints as roots.
1272 extrapolate : {bool, 'periodic', None}, optional
1273 If bool, determines whether to return roots from the polynomial
1274 extrapolated based on first and last intervals, 'periodic' works
1275 the same as False. If None (default), use `self.extrapolate`.
1277 Returns
1278 -------
1279 roots : ndarray
1280 Roots of the polynomial(s).
1282 If the PPoly object describes multiple polynomials, the
1283 return value is an object array whose each element is an
1284 ndarray containing the roots.
1286 Notes
1287 -----
1288 This routine works only on real-valued polynomials.
1290 If the piecewise polynomial contains sections that are
1291 identically zero, the root list will contain the start point
1292 of the corresponding interval, followed by a ``nan`` value.
1294 If the polynomial is discontinuous across a breakpoint, and
1295 there is a sign change across the breakpoint, this is reported
1296 if the `discont` parameter is True.
1298 Examples
1299 --------
1301 Finding roots of ``[x**2 - 1, (x - 1)**2]`` defined on intervals
1302 ``[-2, 1], [1, 2]``:
1304 >>> import numpy as np
1305 >>> from scipy.interpolate import PPoly
1306 >>> pp = PPoly(np.array([[1, -4, 3], [1, 0, 0]]).T, [-2, 1, 2])
1307 >>> pp.solve()
1308 array([-1., 1.])
1309 """
1310 if extrapolate is None:
1311 extrapolate = self.extrapolate
1313 self._ensure_c_contiguous()
1315 if np.issubdtype(self.c.dtype, np.complexfloating):
1316 raise ValueError("Root finding is only for "
1317 "real-valued polynomials")
1319 y = float(y)
1320 r = _ppoly.real_roots(self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1321 self.x, y, bool(discontinuity),
1322 bool(extrapolate))
1323 if self.c.ndim == 2:
1324 return r[0]
1325 else:
1326 r2 = np.empty(prod(self.c.shape[2:]), dtype=object)
1327 # this for-loop is equivalent to ``r2[...] = r``, but that's broken
1328 # in NumPy 1.6.0
1329 for ii, root in enumerate(r):
1330 r2[ii] = root
1332 return r2.reshape(self.c.shape[2:])
1334 def roots(self, discontinuity=True, extrapolate=None):
1335 """
1336 Find real roots of the piecewise polynomial.
1338 Parameters
1339 ----------
1340 discontinuity : bool, optional
1341 Whether to report sign changes across discontinuities at
1342 breakpoints as roots.
1343 extrapolate : {bool, 'periodic', None}, optional
1344 If bool, determines whether to return roots from the polynomial
1345 extrapolated based on first and last intervals, 'periodic' works
1346 the same as False. If None (default), use `self.extrapolate`.
1348 Returns
1349 -------
1350 roots : ndarray
1351 Roots of the polynomial(s).
1353 If the PPoly object describes multiple polynomials, the
1354 return value is an object array whose each element is an
1355 ndarray containing the roots.
1357 See Also
1358 --------
1359 PPoly.solve
1360 """
1361 return self.solve(0, discontinuity, extrapolate)
1363 @classmethod
1364 def from_spline(cls, tck, extrapolate=None):
1365 """
1366 Construct a piecewise polynomial from a spline
1368 Parameters
1369 ----------
1370 tck
1371 A spline, as returned by `splrep` or a BSpline object.
1372 extrapolate : bool or 'periodic', optional
1373 If bool, determines whether to extrapolate to out-of-bounds points
1374 based on first and last intervals, or to return NaNs.
1375 If 'periodic', periodic extrapolation is used. Default is True.
1377 Examples
1378 --------
1379 Construct an interpolating spline and convert it to a `PPoly` instance
1381 >>> import numpy as np
1382 >>> from scipy.interpolate import splrep, PPoly
1383 >>> x = np.linspace(0, 1, 11)
1384 >>> y = np.sin(2*np.pi*x)
1385 >>> tck = splrep(x, y, s=0)
1386 >>> p = PPoly.from_spline(tck)
1387 >>> isinstance(p, PPoly)
1388 True
1390 Note that this function only supports 1D splines out of the box.
1392 If the ``tck`` object represents a parametric spline (e.g. constructed
1393 by `splprep` or a `BSpline` with ``c.ndim > 1``), you will need to loop
1394 over the dimensions manually.
1396 >>> from scipy.interpolate import splprep, splev
1397 >>> t = np.linspace(0, 1, 11)
1398 >>> x = np.sin(2*np.pi*t)
1399 >>> y = np.cos(2*np.pi*t)
1400 >>> (t, c, k), u = splprep([x, y], s=0)
1402 Note that ``c`` is a list of two arrays of length 11.
1404 >>> unew = np.arange(0, 1.01, 0.01)
1405 >>> out = splev(unew, (t, c, k))
1407 To convert this spline to the power basis, we convert each
1408 component of the list of b-spline coefficients, ``c``, into the
1409 corresponding cubic polynomial.
1411 >>> polys = [PPoly.from_spline((t, cj, k)) for cj in c]
1412 >>> polys[0].c.shape
1413 (4, 14)
1415 Note that the coefficients of the polynomials `polys` are in the
1416 power basis and their dimensions reflect just that: here 4 is the order
1417 (degree+1), and 14 is the number of intervals---which is nothing but
1418 the length of the knot array of the original `tck` minus one.
1420 Optionally, we can stack the components into a single `PPoly` along
1421 the third dimension:
1423 >>> cc = np.dstack([p.c for p in polys]) # has shape = (4, 14, 2)
1424 >>> poly = PPoly(cc, polys[0].x)
1425 >>> np.allclose(poly(unew).T, # note the transpose to match `splev`
1426 ... out, atol=1e-15)
1427 True
1429 """
1430 if isinstance(tck, BSpline):
1431 t, c, k = tck.tck
1432 if extrapolate is None:
1433 extrapolate = tck.extrapolate
1434 else:
1435 t, c, k = tck
1437 cvals = np.empty((k + 1, len(t)-1), dtype=c.dtype)
1438 for m in range(k, -1, -1):
1439 y = _fitpack_py.splev(t[:-1], tck, der=m)
1440 cvals[k - m, :] = y/spec.gamma(m+1)
1442 return cls.construct_fast(cvals, t, extrapolate)
1444 @classmethod
1445 def from_bernstein_basis(cls, bp, extrapolate=None):
1446 """
1447 Construct a piecewise polynomial in the power basis
1448 from a polynomial in Bernstein basis.
1450 Parameters
1451 ----------
1452 bp : BPoly
1453 A Bernstein basis polynomial, as created by BPoly
1454 extrapolate : bool or 'periodic', optional
1455 If bool, determines whether to extrapolate to out-of-bounds points
1456 based on first and last intervals, or to return NaNs.
1457 If 'periodic', periodic extrapolation is used. Default is True.
1458 """
1459 if not isinstance(bp, BPoly):
1460 raise TypeError(".from_bernstein_basis only accepts BPoly instances. "
1461 "Got %s instead." % type(bp))
1463 dx = np.diff(bp.x)
1464 k = bp.c.shape[0] - 1 # polynomial order
1466 rest = (None,)*(bp.c.ndim-2)
1468 c = np.zeros_like(bp.c)
1469 for a in range(k+1):
1470 factor = (-1)**a * comb(k, a) * bp.c[a]
1471 for s in range(a, k+1):
1472 val = comb(k-a, s-a) * (-1)**s
1473 c[k-s] += factor * val / dx[(slice(None),)+rest]**s
1475 if extrapolate is None:
1476 extrapolate = bp.extrapolate
1478 return cls.construct_fast(c, bp.x, extrapolate, bp.axis)
1481class BPoly(_PPolyBase):
1482 """Piecewise polynomial in terms of coefficients and breakpoints.
1484 The polynomial between ``x[i]`` and ``x[i + 1]`` is written in the
1485 Bernstein polynomial basis::
1487 S = sum(c[a, i] * b(a, k; x) for a in range(k+1)),
1489 where ``k`` is the degree of the polynomial, and::
1491 b(a, k; x) = binom(k, a) * t**a * (1 - t)**(k - a),
1493 with ``t = (x - x[i]) / (x[i+1] - x[i])`` and ``binom`` is the binomial
1494 coefficient.
1496 Parameters
1497 ----------
1498 c : ndarray, shape (k, m, ...)
1499 Polynomial coefficients, order `k` and `m` intervals
1500 x : ndarray, shape (m+1,)
1501 Polynomial breakpoints. Must be sorted in either increasing or
1502 decreasing order.
1503 extrapolate : bool, optional
1504 If bool, determines whether to extrapolate to out-of-bounds points
1505 based on first and last intervals, or to return NaNs. If 'periodic',
1506 periodic extrapolation is used. Default is True.
1507 axis : int, optional
1508 Interpolation axis. Default is zero.
1510 Attributes
1511 ----------
1512 x : ndarray
1513 Breakpoints.
1514 c : ndarray
1515 Coefficients of the polynomials. They are reshaped
1516 to a 3-D array with the last dimension representing
1517 the trailing dimensions of the original coefficient array.
1518 axis : int
1519 Interpolation axis.
1521 Methods
1522 -------
1523 __call__
1524 extend
1525 derivative
1526 antiderivative
1527 integrate
1528 construct_fast
1529 from_power_basis
1530 from_derivatives
1532 See also
1533 --------
1534 PPoly : piecewise polynomials in the power basis
1536 Notes
1537 -----
1538 Properties of Bernstein polynomials are well documented in the literature,
1539 see for example [1]_ [2]_ [3]_.
1541 References
1542 ----------
1543 .. [1] https://en.wikipedia.org/wiki/Bernstein_polynomial
1545 .. [2] Kenneth I. Joy, Bernstein polynomials,
1546 http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf
1548 .. [3] E. H. Doha, A. H. Bhrawy, and M. A. Saker, Boundary Value Problems,
1549 vol 2011, article ID 829546, :doi:`10.1155/2011/829543`.
1551 Examples
1552 --------
1553 >>> from scipy.interpolate import BPoly
1554 >>> x = [0, 1]
1555 >>> c = [[1], [2], [3]]
1556 >>> bp = BPoly(c, x)
1558 This creates a 2nd order polynomial
1560 .. math::
1562 B(x) = 1 \\times b_{0, 2}(x) + 2 \\times b_{1, 2}(x) + 3 \\times b_{2, 2}(x) \\\\
1563 = 1 \\times (1-x)^2 + 2 \\times 2 x (1 - x) + 3 \\times x^2
1565 """
1567 def _evaluate(self, x, nu, extrapolate, out):
1568 _ppoly.evaluate_bernstein(
1569 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1570 self.x, x, nu, bool(extrapolate), out)
1572 def derivative(self, nu=1):
1573 """
1574 Construct a new piecewise polynomial representing the derivative.
1576 Parameters
1577 ----------
1578 nu : int, optional
1579 Order of derivative to evaluate. Default is 1, i.e., compute the
1580 first derivative. If negative, the antiderivative is returned.
1582 Returns
1583 -------
1584 bp : BPoly
1585 Piecewise polynomial of order k - nu representing the derivative of
1586 this polynomial.
1588 """
1589 if nu < 0:
1590 return self.antiderivative(-nu)
1592 if nu > 1:
1593 bp = self
1594 for k in range(nu):
1595 bp = bp.derivative()
1596 return bp
1598 # reduce order
1599 if nu == 0:
1600 c2 = self.c.copy()
1601 else:
1602 # For a polynomial
1603 # B(x) = \sum_{a=0}^{k} c_a b_{a, k}(x),
1604 # we use the fact that
1605 # b'_{a, k} = k ( b_{a-1, k-1} - b_{a, k-1} ),
1606 # which leads to
1607 # B'(x) = \sum_{a=0}^{k-1} (c_{a+1} - c_a) b_{a, k-1}
1608 #
1609 # finally, for an interval [y, y + dy] with dy != 1,
1610 # we need to correct for an extra power of dy
1612 rest = (None,)*(self.c.ndim-2)
1614 k = self.c.shape[0] - 1
1615 dx = np.diff(self.x)[(None, slice(None))+rest]
1616 c2 = k * np.diff(self.c, axis=0) / dx
1618 if c2.shape[0] == 0:
1619 # derivative of order 0 is zero
1620 c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)
1622 # construct a compatible polynomial
1623 return self.construct_fast(c2, self.x, self.extrapolate, self.axis)
1625 def antiderivative(self, nu=1):
1626 """
1627 Construct a new piecewise polynomial representing the antiderivative.
1629 Parameters
1630 ----------
1631 nu : int, optional
1632 Order of antiderivative to evaluate. Default is 1, i.e., compute
1633 the first integral. If negative, the derivative is returned.
1635 Returns
1636 -------
1637 bp : BPoly
1638 Piecewise polynomial of order k + nu representing the
1639 antiderivative of this polynomial.
1641 Notes
1642 -----
1643 If antiderivative is computed and ``self.extrapolate='periodic'``,
1644 it will be set to False for the returned instance. This is done because
1645 the antiderivative is no longer periodic and its correct evaluation
1646 outside of the initially given x interval is difficult.
1647 """
1648 if nu <= 0:
1649 return self.derivative(-nu)
1651 if nu > 1:
1652 bp = self
1653 for k in range(nu):
1654 bp = bp.antiderivative()
1655 return bp
1657 # Construct the indefinite integrals on individual intervals
1658 c, x = self.c, self.x
1659 k = c.shape[0]
1660 c2 = np.zeros((k+1,) + c.shape[1:], dtype=c.dtype)
1662 c2[1:, ...] = np.cumsum(c, axis=0) / k
1663 delta = x[1:] - x[:-1]
1664 c2 *= delta[(None, slice(None)) + (None,)*(c.ndim-2)]
1666 # Now fix continuity: on the very first interval, take the integration
1667 # constant to be zero; on an interval [x_j, x_{j+1}) with j>0,
1668 # the integration constant is then equal to the jump of the `bp` at x_j.
1669 # The latter is given by the coefficient of B_{n+1, n+1}
1670 # *on the previous interval* (other B. polynomials are zero at the
1671 # breakpoint). Finally, use the fact that BPs form a partition of unity.
1672 c2[:,1:] += np.cumsum(c2[k, :], axis=0)[:-1]
1674 if self.extrapolate == 'periodic':
1675 extrapolate = False
1676 else:
1677 extrapolate = self.extrapolate
1679 return self.construct_fast(c2, x, extrapolate, axis=self.axis)
1681 def integrate(self, a, b, extrapolate=None):
1682 """
1683 Compute a definite integral over a piecewise polynomial.
1685 Parameters
1686 ----------
1687 a : float
1688 Lower integration bound
1689 b : float
1690 Upper integration bound
1691 extrapolate : {bool, 'periodic', None}, optional
1692 Whether to extrapolate to out-of-bounds points based on first
1693 and last intervals, or to return NaNs. If 'periodic', periodic
1694 extrapolation is used. If None (default), use `self.extrapolate`.
1696 Returns
1697 -------
1698 array_like
1699 Definite integral of the piecewise polynomial over [a, b]
1701 """
1702 # XXX: can probably use instead the fact that
1703 # \int_0^{1} B_{j, n}(x) \dx = 1/(n+1)
1704 ib = self.antiderivative()
1705 if extrapolate is None:
1706 extrapolate = self.extrapolate
1708 # ib.extrapolate shouldn't be 'periodic', it is converted to
1709 # False for 'periodic. in antiderivative() call.
1710 if extrapolate != 'periodic':
1711 ib.extrapolate = extrapolate
1713 if extrapolate == 'periodic':
1714 # Split the integral into the part over period (can be several
1715 # of them) and the remaining part.
1717 # For simplicity and clarity convert to a <= b case.
1718 if a <= b:
1719 sign = 1
1720 else:
1721 a, b = b, a
1722 sign = -1
1724 xs, xe = self.x[0], self.x[-1]
1725 period = xe - xs
1726 interval = b - a
1727 n_periods, left = divmod(interval, period)
1728 res = n_periods * (ib(xe) - ib(xs))
1730 # Map a and b to [xs, xe].
1731 a = xs + (a - xs) % period
1732 b = a + left
1734 # If b <= xe then we need to integrate over [a, b], otherwise
1735 # over [a, xe] and from xs to what is remained.
1736 if b <= xe:
1737 res += ib(b) - ib(a)
1738 else:
1739 res += ib(xe) - ib(a) + ib(xs + left + a - xe) - ib(xs)
1741 return sign * res
1742 else:
1743 return ib(b) - ib(a)
1745 def extend(self, c, x):
1746 k = max(self.c.shape[0], c.shape[0])
1747 self.c = self._raise_degree(self.c, k - self.c.shape[0])
1748 c = self._raise_degree(c, k - c.shape[0])
1749 return _PPolyBase.extend(self, c, x)
1750 extend.__doc__ = _PPolyBase.extend.__doc__
1752 @classmethod
1753 def from_power_basis(cls, pp, extrapolate=None):
1754 """
1755 Construct a piecewise polynomial in Bernstein basis
1756 from a power basis polynomial.
1758 Parameters
1759 ----------
1760 pp : PPoly
1761 A piecewise polynomial in the power basis
1762 extrapolate : bool or 'periodic', optional
1763 If bool, determines whether to extrapolate to out-of-bounds points
1764 based on first and last intervals, or to return NaNs.
1765 If 'periodic', periodic extrapolation is used. Default is True.
1766 """
1767 if not isinstance(pp, PPoly):
1768 raise TypeError(".from_power_basis only accepts PPoly instances. "
1769 "Got %s instead." % type(pp))
1771 dx = np.diff(pp.x)
1772 k = pp.c.shape[0] - 1 # polynomial order
1774 rest = (None,)*(pp.c.ndim-2)
1776 c = np.zeros_like(pp.c)
1777 for a in range(k+1):
1778 factor = pp.c[a] / comb(k, k-a) * dx[(slice(None),)+rest]**(k-a)
1779 for j in range(k-a, k+1):
1780 c[j] += factor * comb(j, k-a)
1782 if extrapolate is None:
1783 extrapolate = pp.extrapolate
1785 return cls.construct_fast(c, pp.x, extrapolate, pp.axis)
1787 @classmethod
1788 def from_derivatives(cls, xi, yi, orders=None, extrapolate=None):
1789 """Construct a piecewise polynomial in the Bernstein basis,
1790 compatible with the specified values and derivatives at breakpoints.
1792 Parameters
1793 ----------
1794 xi : array_like
1795 sorted 1-D array of x-coordinates
1796 yi : array_like or list of array_likes
1797 ``yi[i][j]`` is the ``j``th derivative known at ``xi[i]``
1798 orders : None or int or array_like of ints. Default: None.
1799 Specifies the degree of local polynomials. If not None, some
1800 derivatives are ignored.
1801 extrapolate : bool or 'periodic', optional
1802 If bool, determines whether to extrapolate to out-of-bounds points
1803 based on first and last intervals, or to return NaNs.
1804 If 'periodic', periodic extrapolation is used. Default is True.
1806 Notes
1807 -----
1808 If ``k`` derivatives are specified at a breakpoint ``x``, the
1809 constructed polynomial is exactly ``k`` times continuously
1810 differentiable at ``x``, unless the ``order`` is provided explicitly.
1811 In the latter case, the smoothness of the polynomial at
1812 the breakpoint is controlled by the ``order``.
1814 Deduces the number of derivatives to match at each end
1815 from ``order`` and the number of derivatives available. If
1816 possible it uses the same number of derivatives from
1817 each end; if the number is odd it tries to take the
1818 extra one from y2. In any case if not enough derivatives
1819 are available at one end or another it draws enough to
1820 make up the total from the other end.
1822 If the order is too high and not enough derivatives are available,
1823 an exception is raised.
1825 Examples
1826 --------
1828 >>> from scipy.interpolate import BPoly
1829 >>> BPoly.from_derivatives([0, 1], [[1, 2], [3, 4]])
1831 Creates a polynomial `f(x)` of degree 3, defined on `[0, 1]`
1832 such that `f(0) = 1, df/dx(0) = 2, f(1) = 3, df/dx(1) = 4`
1834 >>> BPoly.from_derivatives([0, 1, 2], [[0, 1], [0], [2]])
1836 Creates a piecewise polynomial `f(x)`, such that
1837 `f(0) = f(1) = 0`, `f(2) = 2`, and `df/dx(0) = 1`.
1838 Based on the number of derivatives provided, the order of the
1839 local polynomials is 2 on `[0, 1]` and 1 on `[1, 2]`.
1840 Notice that no restriction is imposed on the derivatives at
1841 ``x = 1`` and ``x = 2``.
1843 Indeed, the explicit form of the polynomial is::
1845 f(x) = | x * (1 - x), 0 <= x < 1
1846 | 2 * (x - 1), 1 <= x <= 2
1848 So that f'(1-0) = -1 and f'(1+0) = 2
1850 """
1851 xi = np.asarray(xi)
1852 if len(xi) != len(yi):
1853 raise ValueError("xi and yi need to have the same length")
1854 if np.any(xi[1:] - xi[:1] <= 0):
1855 raise ValueError("x coordinates are not in increasing order")
1857 # number of intervals
1858 m = len(xi) - 1
1860 # global poly order is k-1, local orders are <=k and can vary
1861 try:
1862 k = max(len(yi[i]) + len(yi[i+1]) for i in range(m))
1863 except TypeError as e:
1864 raise ValueError(
1865 "Using a 1-D array for y? Please .reshape(-1, 1)."
1866 ) from e
1868 if orders is None:
1869 orders = [None] * m
1870 else:
1871 if isinstance(orders, (int, np.integer)):
1872 orders = [orders] * m
1873 k = max(k, max(orders))
1875 if any(o <= 0 for o in orders):
1876 raise ValueError("Orders must be positive.")
1878 c = []
1879 for i in range(m):
1880 y1, y2 = yi[i], yi[i+1]
1881 if orders[i] is None:
1882 n1, n2 = len(y1), len(y2)
1883 else:
1884 n = orders[i]+1
1885 n1 = min(n//2, len(y1))
1886 n2 = min(n - n1, len(y2))
1887 n1 = min(n - n2, len(y2))
1888 if n1+n2 != n:
1889 mesg = ("Point %g has %d derivatives, point %g"
1890 " has %d derivatives, but order %d requested" % (
1891 xi[i], len(y1), xi[i+1], len(y2), orders[i]))
1892 raise ValueError(mesg)
1894 if not (n1 <= len(y1) and n2 <= len(y2)):
1895 raise ValueError("`order` input incompatible with"
1896 " length y1 or y2.")
1898 b = BPoly._construct_from_derivatives(xi[i], xi[i+1],
1899 y1[:n1], y2[:n2])
1900 if len(b) < k:
1901 b = BPoly._raise_degree(b, k - len(b))
1902 c.append(b)
1904 c = np.asarray(c)
1905 return cls(c.swapaxes(0, 1), xi, extrapolate)
1907 @staticmethod
1908 def _construct_from_derivatives(xa, xb, ya, yb):
1909 r"""Compute the coefficients of a polynomial in the Bernstein basis
1910 given the values and derivatives at the edges.
1912 Return the coefficients of a polynomial in the Bernstein basis
1913 defined on ``[xa, xb]`` and having the values and derivatives at the
1914 endpoints `xa` and `xb` as specified by `ya`` and `yb`.
1915 The polynomial constructed is of the minimal possible degree, i.e.,
1916 if the lengths of `ya` and `yb` are `na` and `nb`, the degree
1917 of the polynomial is ``na + nb - 1``.
1919 Parameters
1920 ----------
1921 xa : float
1922 Left-hand end point of the interval
1923 xb : float
1924 Right-hand end point of the interval
1925 ya : array_like
1926 Derivatives at `xa`. `ya[0]` is the value of the function, and
1927 `ya[i]` for ``i > 0`` is the value of the ``i``th derivative.
1928 yb : array_like
1929 Derivatives at `xb`.
1931 Returns
1932 -------
1933 array
1934 coefficient array of a polynomial having specified derivatives
1936 Notes
1937 -----
1938 This uses several facts from life of Bernstein basis functions.
1939 First of all,
1941 .. math:: b'_{a, n} = n (b_{a-1, n-1} - b_{a, n-1})
1943 If B(x) is a linear combination of the form
1945 .. math:: B(x) = \sum_{a=0}^{n} c_a b_{a, n},
1947 then :math: B'(x) = n \sum_{a=0}^{n-1} (c_{a+1} - c_{a}) b_{a, n-1}.
1948 Iterating the latter one, one finds for the q-th derivative
1950 .. math:: B^{q}(x) = n!/(n-q)! \sum_{a=0}^{n-q} Q_a b_{a, n-q},
1952 with
1954 .. math:: Q_a = \sum_{j=0}^{q} (-)^{j+q} comb(q, j) c_{j+a}
1956 This way, only `a=0` contributes to :math: `B^{q}(x = xa)`, and
1957 `c_q` are found one by one by iterating `q = 0, ..., na`.
1959 At ``x = xb`` it's the same with ``a = n - q``.
1961 """
1962 ya, yb = np.asarray(ya), np.asarray(yb)
1963 if ya.shape[1:] != yb.shape[1:]:
1964 raise ValueError('Shapes of ya {} and yb {} are incompatible'
1965 .format(ya.shape, yb.shape))
1967 dta, dtb = ya.dtype, yb.dtype
1968 if (np.issubdtype(dta, np.complexfloating) or
1969 np.issubdtype(dtb, np.complexfloating)):
1970 dt = np.complex_
1971 else:
1972 dt = np.float_
1974 na, nb = len(ya), len(yb)
1975 n = na + nb
1977 c = np.empty((na+nb,) + ya.shape[1:], dtype=dt)
1979 # compute coefficients of a polynomial degree na+nb-1
1980 # walk left-to-right
1981 for q in range(0, na):
1982 c[q] = ya[q] / spec.poch(n - q, q) * (xb - xa)**q
1983 for j in range(0, q):
1984 c[q] -= (-1)**(j+q) * comb(q, j) * c[j]
1986 # now walk right-to-left
1987 for q in range(0, nb):
1988 c[-q-1] = yb[q] / spec.poch(n - q, q) * (-1)**q * (xb - xa)**q
1989 for j in range(0, q):
1990 c[-q-1] -= (-1)**(j+1) * comb(q, j+1) * c[-q+j]
1992 return c
1994 @staticmethod
1995 def _raise_degree(c, d):
1996 r"""Raise a degree of a polynomial in the Bernstein basis.
1998 Given the coefficients of a polynomial degree `k`, return (the
1999 coefficients of) the equivalent polynomial of degree `k+d`.
2001 Parameters
2002 ----------
2003 c : array_like
2004 coefficient array, 1-D
2005 d : integer
2007 Returns
2008 -------
2009 array
2010 coefficient array, 1-D array of length `c.shape[0] + d`
2012 Notes
2013 -----
2014 This uses the fact that a Bernstein polynomial `b_{a, k}` can be
2015 identically represented as a linear combination of polynomials of
2016 a higher degree `k+d`:
2018 .. math:: b_{a, k} = comb(k, a) \sum_{j=0}^{d} b_{a+j, k+d} \
2019 comb(d, j) / comb(k+d, a+j)
2021 """
2022 if d == 0:
2023 return c
2025 k = c.shape[0] - 1
2026 out = np.zeros((c.shape[0] + d,) + c.shape[1:], dtype=c.dtype)
2028 for a in range(c.shape[0]):
2029 f = c[a] * comb(k, a)
2030 for j in range(d+1):
2031 out[a+j] += f * comb(d, j) / comb(k+d, a+j)
2032 return out
2035class NdPPoly:
2036 """
2037 Piecewise tensor product polynomial
2039 The value at point ``xp = (x', y', z', ...)`` is evaluated by first
2040 computing the interval indices `i` such that::
2042 x[0][i[0]] <= x' < x[0][i[0]+1]
2043 x[1][i[1]] <= y' < x[1][i[1]+1]
2044 ...
2046 and then computing::
2048 S = sum(c[k0-m0-1,...,kn-mn-1,i[0],...,i[n]]
2049 * (xp[0] - x[0][i[0]])**m0
2050 * ...
2051 * (xp[n] - x[n][i[n]])**mn
2052 for m0 in range(k[0]+1)
2053 ...
2054 for mn in range(k[n]+1))
2056 where ``k[j]`` is the degree of the polynomial in dimension j. This
2057 representation is the piecewise multivariate power basis.
2059 Parameters
2060 ----------
2061 c : ndarray, shape (k0, ..., kn, m0, ..., mn, ...)
2062 Polynomial coefficients, with polynomial order `kj` and
2063 `mj+1` intervals for each dimension `j`.
2064 x : ndim-tuple of ndarrays, shapes (mj+1,)
2065 Polynomial breakpoints for each dimension. These must be
2066 sorted in increasing order.
2067 extrapolate : bool, optional
2068 Whether to extrapolate to out-of-bounds points based on first
2069 and last intervals, or to return NaNs. Default: True.
2071 Attributes
2072 ----------
2073 x : tuple of ndarrays
2074 Breakpoints.
2075 c : ndarray
2076 Coefficients of the polynomials.
2078 Methods
2079 -------
2080 __call__
2081 derivative
2082 antiderivative
2083 integrate
2084 integrate_1d
2085 construct_fast
2087 See also
2088 --------
2089 PPoly : piecewise polynomials in 1D
2091 Notes
2092 -----
2093 High-order polynomials in the power basis can be numerically
2094 unstable.
2096 """
2098 def __init__(self, c, x, extrapolate=None):
2099 self.x = tuple(np.ascontiguousarray(v, dtype=np.float64) for v in x)
2100 self.c = np.asarray(c)
2101 if extrapolate is None:
2102 extrapolate = True
2103 self.extrapolate = bool(extrapolate)
2105 ndim = len(self.x)
2106 if any(v.ndim != 1 for v in self.x):
2107 raise ValueError("x arrays must all be 1-dimensional")
2108 if any(v.size < 2 for v in self.x):
2109 raise ValueError("x arrays must all contain at least 2 points")
2110 if c.ndim < 2*ndim:
2111 raise ValueError("c must have at least 2*len(x) dimensions")
2112 if any(np.any(v[1:] - v[:-1] < 0) for v in self.x):
2113 raise ValueError("x-coordinates are not in increasing order")
2114 if any(a != b.size - 1 for a, b in zip(c.shape[ndim:2*ndim], self.x)):
2115 raise ValueError("x and c do not agree on the number of intervals")
2117 dtype = self._get_dtype(self.c.dtype)
2118 self.c = np.ascontiguousarray(self.c, dtype=dtype)
2120 @classmethod
2121 def construct_fast(cls, c, x, extrapolate=None):
2122 """
2123 Construct the piecewise polynomial without making checks.
2125 Takes the same parameters as the constructor. Input arguments
2126 ``c`` and ``x`` must be arrays of the correct shape and type. The
2127 ``c`` array can only be of dtypes float and complex, and ``x``
2128 array must have dtype float.
2130 """
2131 self = object.__new__(cls)
2132 self.c = c
2133 self.x = x
2134 if extrapolate is None:
2135 extrapolate = True
2136 self.extrapolate = extrapolate
2137 return self
2139 def _get_dtype(self, dtype):
2140 if np.issubdtype(dtype, np.complexfloating) \
2141 or np.issubdtype(self.c.dtype, np.complexfloating):
2142 return np.complex_
2143 else:
2144 return np.float_
2146 def _ensure_c_contiguous(self):
2147 if not self.c.flags.c_contiguous:
2148 self.c = self.c.copy()
2149 if not isinstance(self.x, tuple):
2150 self.x = tuple(self.x)
2152 def __call__(self, x, nu=None, extrapolate=None):
2153 """
2154 Evaluate the piecewise polynomial or its derivative
2156 Parameters
2157 ----------
2158 x : array-like
2159 Points to evaluate the interpolant at.
2160 nu : tuple, optional
2161 Orders of derivatives to evaluate. Each must be non-negative.
2162 extrapolate : bool, optional
2163 Whether to extrapolate to out-of-bounds points based on first
2164 and last intervals, or to return NaNs.
2166 Returns
2167 -------
2168 y : array-like
2169 Interpolated values. Shape is determined by replacing
2170 the interpolation axis in the original array with the shape of x.
2172 Notes
2173 -----
2174 Derivatives are evaluated piecewise for each polynomial
2175 segment, even if the polynomial is not differentiable at the
2176 breakpoints. The polynomial intervals are considered half-open,
2177 ``[a, b)``, except for the last interval which is closed
2178 ``[a, b]``.
2180 """
2181 if extrapolate is None:
2182 extrapolate = self.extrapolate
2183 else:
2184 extrapolate = bool(extrapolate)
2186 ndim = len(self.x)
2188 x = _ndim_coords_from_arrays(x)
2189 x_shape = x.shape
2190 x = np.ascontiguousarray(x.reshape(-1, x.shape[-1]), dtype=np.float_)
2192 if nu is None:
2193 nu = np.zeros((ndim,), dtype=np.intc)
2194 else:
2195 nu = np.asarray(nu, dtype=np.intc)
2196 if nu.ndim != 1 or nu.shape[0] != ndim:
2197 raise ValueError("invalid number of derivative orders nu")
2199 dim1 = prod(self.c.shape[:ndim])
2200 dim2 = prod(self.c.shape[ndim:2*ndim])
2201 dim3 = prod(self.c.shape[2*ndim:])
2202 ks = np.array(self.c.shape[:ndim], dtype=np.intc)
2204 out = np.empty((x.shape[0], dim3), dtype=self.c.dtype)
2205 self._ensure_c_contiguous()
2207 _ppoly.evaluate_nd(self.c.reshape(dim1, dim2, dim3),
2208 self.x,
2209 ks,
2210 x,
2211 nu,
2212 bool(extrapolate),
2213 out)
2215 return out.reshape(x_shape[:-1] + self.c.shape[2*ndim:])
2217 def _derivative_inplace(self, nu, axis):
2218 """
2219 Compute 1-D derivative along a selected dimension in-place
2220 May result to non-contiguous c array.
2221 """
2222 if nu < 0:
2223 return self._antiderivative_inplace(-nu, axis)
2225 ndim = len(self.x)
2226 axis = axis % ndim
2228 # reduce order
2229 if nu == 0:
2230 # noop
2231 return
2232 else:
2233 sl = [slice(None)]*ndim
2234 sl[axis] = slice(None, -nu, None)
2235 c2 = self.c[tuple(sl)]
2237 if c2.shape[axis] == 0:
2238 # derivative of order 0 is zero
2239 shp = list(c2.shape)
2240 shp[axis] = 1
2241 c2 = np.zeros(shp, dtype=c2.dtype)
2243 # multiply by the correct rising factorials
2244 factor = spec.poch(np.arange(c2.shape[axis], 0, -1), nu)
2245 sl = [None]*c2.ndim
2246 sl[axis] = slice(None)
2247 c2 *= factor[tuple(sl)]
2249 self.c = c2
2251 def _antiderivative_inplace(self, nu, axis):
2252 """
2253 Compute 1-D antiderivative along a selected dimension
2254 May result to non-contiguous c array.
2255 """
2256 if nu <= 0:
2257 return self._derivative_inplace(-nu, axis)
2259 ndim = len(self.x)
2260 axis = axis % ndim
2262 perm = list(range(ndim))
2263 perm[0], perm[axis] = perm[axis], perm[0]
2264 perm = perm + list(range(ndim, self.c.ndim))
2266 c = self.c.transpose(perm)
2268 c2 = np.zeros((c.shape[0] + nu,) + c.shape[1:],
2269 dtype=c.dtype)
2270 c2[:-nu] = c
2272 # divide by the correct rising factorials
2273 factor = spec.poch(np.arange(c.shape[0], 0, -1), nu)
2274 c2[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]
2276 # fix continuity of added degrees of freedom
2277 perm2 = list(range(c2.ndim))
2278 perm2[1], perm2[ndim+axis] = perm2[ndim+axis], perm2[1]
2280 c2 = c2.transpose(perm2)
2281 c2 = c2.copy()
2282 _ppoly.fix_continuity(c2.reshape(c2.shape[0], c2.shape[1], -1),
2283 self.x[axis], nu-1)
2285 c2 = c2.transpose(perm2)
2286 c2 = c2.transpose(perm)
2288 # Done
2289 self.c = c2
2291 def derivative(self, nu):
2292 """
2293 Construct a new piecewise polynomial representing the derivative.
2295 Parameters
2296 ----------
2297 nu : ndim-tuple of int
2298 Order of derivatives to evaluate for each dimension.
2299 If negative, the antiderivative is returned.
2301 Returns
2302 -------
2303 pp : NdPPoly
2304 Piecewise polynomial of orders (k[0] - nu[0], ..., k[n] - nu[n])
2305 representing the derivative of this polynomial.
2307 Notes
2308 -----
2309 Derivatives are evaluated piecewise for each polynomial
2310 segment, even if the polynomial is not differentiable at the
2311 breakpoints. The polynomial intervals in each dimension are
2312 considered half-open, ``[a, b)``, except for the last interval
2313 which is closed ``[a, b]``.
2315 """
2316 p = self.construct_fast(self.c.copy(), self.x, self.extrapolate)
2318 for axis, n in enumerate(nu):
2319 p._derivative_inplace(n, axis)
2321 p._ensure_c_contiguous()
2322 return p
2324 def antiderivative(self, nu):
2325 """
2326 Construct a new piecewise polynomial representing the antiderivative.
2328 Antiderivative is also the indefinite integral of the function,
2329 and derivative is its inverse operation.
2331 Parameters
2332 ----------
2333 nu : ndim-tuple of int
2334 Order of derivatives to evaluate for each dimension.
2335 If negative, the derivative is returned.
2337 Returns
2338 -------
2339 pp : PPoly
2340 Piecewise polynomial of order k2 = k + n representing
2341 the antiderivative of this polynomial.
2343 Notes
2344 -----
2345 The antiderivative returned by this function is continuous and
2346 continuously differentiable to order n-1, up to floating point
2347 rounding error.
2349 """
2350 p = self.construct_fast(self.c.copy(), self.x, self.extrapolate)
2352 for axis, n in enumerate(nu):
2353 p._antiderivative_inplace(n, axis)
2355 p._ensure_c_contiguous()
2356 return p
2358 def integrate_1d(self, a, b, axis, extrapolate=None):
2359 r"""
2360 Compute NdPPoly representation for one dimensional definite integral
2362 The result is a piecewise polynomial representing the integral:
2364 .. math::
2366 p(y, z, ...) = \int_a^b dx\, p(x, y, z, ...)
2368 where the dimension integrated over is specified with the
2369 `axis` parameter.
2371 Parameters
2372 ----------
2373 a, b : float
2374 Lower and upper bound for integration.
2375 axis : int
2376 Dimension over which to compute the 1-D integrals
2377 extrapolate : bool, optional
2378 Whether to extrapolate to out-of-bounds points based on first
2379 and last intervals, or to return NaNs.
2381 Returns
2382 -------
2383 ig : NdPPoly or array-like
2384 Definite integral of the piecewise polynomial over [a, b].
2385 If the polynomial was 1D, an array is returned,
2386 otherwise, an NdPPoly object.
2388 """
2389 if extrapolate is None:
2390 extrapolate = self.extrapolate
2391 else:
2392 extrapolate = bool(extrapolate)
2394 ndim = len(self.x)
2395 axis = int(axis) % ndim
2397 # reuse 1-D integration routines
2398 c = self.c
2399 swap = list(range(c.ndim))
2400 swap.insert(0, swap[axis])
2401 del swap[axis + 1]
2402 swap.insert(1, swap[ndim + axis])
2403 del swap[ndim + axis + 1]
2405 c = c.transpose(swap)
2406 p = PPoly.construct_fast(c.reshape(c.shape[0], c.shape[1], -1),
2407 self.x[axis],
2408 extrapolate=extrapolate)
2409 out = p.integrate(a, b, extrapolate=extrapolate)
2411 # Construct result
2412 if ndim == 1:
2413 return out.reshape(c.shape[2:])
2414 else:
2415 c = out.reshape(c.shape[2:])
2416 x = self.x[:axis] + self.x[axis+1:]
2417 return self.construct_fast(c, x, extrapolate=extrapolate)
2419 def integrate(self, ranges, extrapolate=None):
2420 """
2421 Compute a definite integral over a piecewise polynomial.
2423 Parameters
2424 ----------
2425 ranges : ndim-tuple of 2-tuples float
2426 Sequence of lower and upper bounds for each dimension,
2427 ``[(a[0], b[0]), ..., (a[ndim-1], b[ndim-1])]``
2428 extrapolate : bool, optional
2429 Whether to extrapolate to out-of-bounds points based on first
2430 and last intervals, or to return NaNs.
2432 Returns
2433 -------
2434 ig : array_like
2435 Definite integral of the piecewise polynomial over
2436 [a[0], b[0]] x ... x [a[ndim-1], b[ndim-1]]
2438 """
2440 ndim = len(self.x)
2442 if extrapolate is None:
2443 extrapolate = self.extrapolate
2444 else:
2445 extrapolate = bool(extrapolate)
2447 if not hasattr(ranges, '__len__') or len(ranges) != ndim:
2448 raise ValueError("Range not a sequence of correct length")
2450 self._ensure_c_contiguous()
2452 # Reuse 1D integration routine
2453 c = self.c
2454 for n, (a, b) in enumerate(ranges):
2455 swap = list(range(c.ndim))
2456 swap.insert(1, swap[ndim - n])
2457 del swap[ndim - n + 1]
2459 c = c.transpose(swap)
2461 p = PPoly.construct_fast(c, self.x[n], extrapolate=extrapolate)
2462 out = p.integrate(a, b, extrapolate=extrapolate)
2463 c = out.reshape(c.shape[2:])
2465 return c