Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_interpolate.py: 11%
786 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1__all__ = ['interp1d', 'interp2d', 'lagrange', 'PPoly', 'BPoly', 'NdPPoly']
4import numpy as np
5from numpy import (array, transpose, searchsorted, atleast_1d, atleast_2d,
6 ravel, poly1d, asarray, intp)
8import scipy.special as spec
9from scipy.special import comb
10from scipy._lib._util import prod
12from . import _fitpack_py
13from . import dfitpack
14from . import _fitpack
15from ._polyint import _Interpolator1D
16from . import _ppoly
17from ._fitpack2 import RectBivariateSpline
18from .interpnd import _ndim_coords_from_arrays
19from ._bsplines import make_interp_spline, BSpline
22def lagrange(x, w):
23 r"""
24 Return a Lagrange interpolating polynomial.
26 Given two 1-D arrays `x` and `w,` returns the Lagrange interpolating
27 polynomial through the points ``(x, w)``.
29 Warning: This implementation is numerically unstable. Do not expect to
30 be able to use more than about 20 points even if they are chosen optimally.
32 Parameters
33 ----------
34 x : array_like
35 `x` represents the x-coordinates of a set of datapoints.
36 w : array_like
37 `w` represents the y-coordinates of a set of datapoints, i.e., f(`x`).
39 Returns
40 -------
41 lagrange : `numpy.poly1d` instance
42 The Lagrange interpolating polynomial.
44 Examples
45 --------
46 Interpolate :math:`f(x) = x^3` by 3 points.
48 >>> import numpy as np
49 >>> from scipy.interpolate import lagrange
50 >>> x = np.array([0, 1, 2])
51 >>> y = x**3
52 >>> poly = lagrange(x, y)
54 Since there are only 3 points, Lagrange polynomial has degree 2. Explicitly,
55 it is given by
57 .. math::
59 \begin{aligned}
60 L(x) &= 1\times \frac{x (x - 2)}{-1} + 8\times \frac{x (x-1)}{2} \\
61 &= x (-2 + 3x)
62 \end{aligned}
64 >>> from numpy.polynomial.polynomial import Polynomial
65 >>> Polynomial(poly.coef[::-1]).coef
66 array([ 0., -2., 3.])
68 >>> import matplotlib.pyplot as plt
69 >>> x_new = np.arange(0, 2.1, 0.1)
70 >>> plt.scatter(x, y, label='data')
71 >>> plt.plot(x_new, Polynomial(poly.coef[::-1])(x_new), label='Polynomial')
72 >>> plt.plot(x_new, 3*x_new**2 - 2*x_new + 0*x_new,
73 ... label=r"$3 x^2 - 2 x$", linestyle='-.')
74 >>> plt.legend()
75 >>> plt.show()
77 """
79 M = len(x)
80 p = poly1d(0.0)
81 for j in range(M):
82 pt = poly1d(w[j])
83 for k in range(M):
84 if k == j:
85 continue
86 fac = x[j]-x[k]
87 pt *= poly1d([1.0, -x[k]])/fac
88 p += pt
89 return p
92# !! Need to find argument for keeping initialize. If it isn't
93# !! found, get rid of it!
96dep_mesg = """\
97`interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy 1.12.0.
99For legacy code, nearly bug-for-bug compatible replacements are
100`RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
101scattered 2D data.
103In new code, for regular grids use `RegularGridInterpolator` instead.
104For scattered data, prefer `LinearNDInterpolator` or
105`CloughTocher2DInterpolator`.
107For more details see
108`https://gist.github.com/ev-br/8544371b40f414b7eaf3fe6217209bff`
109"""
111class interp2d:
112 """
113 interp2d(x, y, z, kind='linear', copy=True, bounds_error=False,
114 fill_value=None)
116 .. deprecated:: 1.10.0
118 `interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy
119 1.12.0.
121 For legacy code, nearly bug-for-bug compatible replacements are
122 `RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
123 scattered 2D data.
125 In new code, for regular grids use `RegularGridInterpolator` instead.
126 For scattered data, prefer `LinearNDInterpolator` or
127 `CloughTocher2DInterpolator`.
129 For more details see
130 `https://gist.github.com/ev-br/8544371b40f414b7eaf3fe6217209bff`
133 Interpolate over a 2-D grid.
135 `x`, `y` and `z` are arrays of values used to approximate some function
136 f: ``z = f(x, y)`` which returns a scalar value `z`. This class returns a
137 function whose call method uses spline interpolation to find the value
138 of new points.
140 If `x` and `y` represent a regular grid, consider using
141 `RectBivariateSpline`.
143 If `z` is a vector value, consider using `interpn`.
145 Note that calling `interp2d` with NaNs present in input values, or with
146 decreasing values in `x` an `y` results in undefined behaviour.
148 Methods
149 -------
150 __call__
152 Parameters
153 ----------
154 x, y : array_like
155 Arrays defining the data point coordinates.
156 The data point coordinates need to be sorted by increasing order.
158 If the points lie on a regular grid, `x` can specify the column
159 coordinates and `y` the row coordinates, for example::
161 >>> x = [0,1,2]; y = [0,3]; z = [[1,2,3], [4,5,6]]
163 Otherwise, `x` and `y` must specify the full coordinates for each
164 point, for example::
166 >>> x = [0,1,2,0,1,2]; y = [0,0,0,3,3,3]; z = [1,4,2,5,3,6]
168 If `x` and `y` are multidimensional, they are flattened before use.
169 z : array_like
170 The values of the function to interpolate at the data points. If
171 `z` is a multidimensional array, it is flattened before use assuming
172 Fortran-ordering (order='F'). The length of a flattened `z` array
173 is either len(`x`)*len(`y`) if `x` and `y` specify the column and
174 row coordinates or ``len(z) == len(x) == len(y)`` if `x` and `y`
175 specify coordinates for each point.
176 kind : {'linear', 'cubic', 'quintic'}, optional
177 The kind of spline interpolation to use. Default is 'linear'.
178 copy : bool, optional
179 If True, the class makes internal copies of x, y and z.
180 If False, references may be used. The default is to copy.
181 bounds_error : bool, optional
182 If True, when interpolated values are requested outside of the
183 domain of the input data (x,y), a ValueError is raised.
184 If False, then `fill_value` is used.
185 fill_value : number, optional
186 If provided, the value to use for points outside of the
187 interpolation domain. If omitted (None), values outside
188 the domain are extrapolated via nearest-neighbor extrapolation.
190 See Also
191 --------
192 RectBivariateSpline :
193 Much faster 2-D interpolation if your input data is on a grid
194 bisplrep, bisplev :
195 Spline interpolation based on FITPACK
196 BivariateSpline : a more recent wrapper of the FITPACK routines
197 interp1d : 1-D version of this function
198 RegularGridInterpolator : interpolation on a regular or rectilinear grid
199 in arbitrary dimensions.
200 interpn : Multidimensional interpolation on regular grids (wraps
201 `RegularGridInterpolator` and `RectBivariateSpline`).
203 Notes
204 -----
205 The minimum number of data points required along the interpolation
206 axis is ``(k+1)**2``, with k=1 for linear, k=3 for cubic and k=5 for
207 quintic interpolation.
209 The interpolator is constructed by `bisplrep`, with a smoothing factor
210 of 0. If more control over smoothing is needed, `bisplrep` should be
211 used directly.
213 The coordinates of the data points to interpolate `xnew` and `ynew`
214 have to be sorted by ascending order.
215 `interp2d` is legacy and is not
216 recommended for use in new code. New code should use
217 `RegularGridInterpolator` instead.
219 Examples
220 --------
221 Construct a 2-D grid and interpolate on it:
223 >>> import numpy as np
224 >>> from scipy import interpolate
225 >>> x = np.arange(-5.01, 5.01, 0.25)
226 >>> y = np.arange(-5.01, 5.01, 0.25)
227 >>> xx, yy = np.meshgrid(x, y)
228 >>> z = np.sin(xx**2+yy**2)
229 >>> f = interpolate.interp2d(x, y, z, kind='cubic')
231 Now use the obtained interpolation function and plot the result:
233 >>> import matplotlib.pyplot as plt
234 >>> xnew = np.arange(-5.01, 5.01, 1e-2)
235 >>> ynew = np.arange(-5.01, 5.01, 1e-2)
236 >>> znew = f(xnew, ynew)
237 >>> plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
238 >>> plt.show()
239 """
241 @np.deprecate(old_name='interp2d', message=dep_mesg)
242 def __init__(self, x, y, z, kind='linear', copy=True, bounds_error=False,
243 fill_value=None):
244 x = ravel(x)
245 y = ravel(y)
246 z = asarray(z)
248 rectangular_grid = (z.size == len(x) * len(y))
249 if rectangular_grid:
250 if z.ndim == 2:
251 if z.shape != (len(y), len(x)):
252 raise ValueError("When on a regular grid with x.size = m "
253 "and y.size = n, if z.ndim == 2, then z "
254 "must have shape (n, m)")
255 if not np.all(x[1:] >= x[:-1]):
256 j = np.argsort(x)
257 x = x[j]
258 z = z[:, j]
259 if not np.all(y[1:] >= y[:-1]):
260 j = np.argsort(y)
261 y = y[j]
262 z = z[j, :]
263 z = ravel(z.T)
264 else:
265 z = ravel(z)
266 if len(x) != len(y):
267 raise ValueError(
268 "x and y must have equal lengths for non rectangular grid")
269 if len(z) != len(x):
270 raise ValueError(
271 "Invalid length for input z for non rectangular grid")
273 interpolation_types = {'linear': 1, 'cubic': 3, 'quintic': 5}
274 try:
275 kx = ky = interpolation_types[kind]
276 except KeyError as e:
277 raise ValueError(
278 f"Unsupported interpolation type {repr(kind)}, must be "
279 f"either of {', '.join(map(repr, interpolation_types))}."
280 ) from e
282 if not rectangular_grid:
283 # TODO: surfit is really not meant for interpolation!
284 self.tck = _fitpack_py.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
285 else:
286 nx, tx, ny, ty, c, fp, ier = dfitpack.regrid_smth(
287 x, y, z, None, None, None, None,
288 kx=kx, ky=ky, s=0.0)
289 self.tck = (tx[:nx], ty[:ny], c[:(nx - kx - 1) * (ny - ky - 1)],
290 kx, ky)
292 self.bounds_error = bounds_error
293 self.fill_value = fill_value
294 self.x, self.y, self.z = [array(a, copy=copy) for a in (x, y, z)]
296 self.x_min, self.x_max = np.amin(x), np.amax(x)
297 self.y_min, self.y_max = np.amin(y), np.amax(y)
299 @np.deprecate(old_name='interp2d', message=dep_mesg)
300 def __call__(self, x, y, dx=0, dy=0, assume_sorted=False):
301 """Interpolate the function.
303 Parameters
304 ----------
305 x : 1-D array
306 x-coordinates of the mesh on which to interpolate.
307 y : 1-D array
308 y-coordinates of the mesh on which to interpolate.
309 dx : int >= 0, < kx
310 Order of partial derivatives in x.
311 dy : int >= 0, < ky
312 Order of partial derivatives in y.
313 assume_sorted : bool, optional
314 If False, values of `x` and `y` can be in any order and they are
315 sorted first.
316 If True, `x` and `y` have to be arrays of monotonically
317 increasing values.
319 Returns
320 -------
321 z : 2-D array with shape (len(y), len(x))
322 The interpolated values.
323 """
325 x = atleast_1d(x)
326 y = atleast_1d(y)
328 if x.ndim != 1 or y.ndim != 1:
329 raise ValueError("x and y should both be 1-D arrays")
331 if not assume_sorted:
332 x = np.sort(x, kind="mergesort")
333 y = np.sort(y, kind="mergesort")
335 if self.bounds_error or self.fill_value is not None:
336 out_of_bounds_x = (x < self.x_min) | (x > self.x_max)
337 out_of_bounds_y = (y < self.y_min) | (y > self.y_max)
339 any_out_of_bounds_x = np.any(out_of_bounds_x)
340 any_out_of_bounds_y = np.any(out_of_bounds_y)
342 if self.bounds_error and (any_out_of_bounds_x or any_out_of_bounds_y):
343 raise ValueError("Values out of range; x must be in %r, y in %r"
344 % ((self.x_min, self.x_max),
345 (self.y_min, self.y_max)))
347 z = _fitpack_py.bisplev(x, y, self.tck, dx, dy)
348 z = atleast_2d(z)
349 z = transpose(z)
351 if self.fill_value is not None:
352 if any_out_of_bounds_x:
353 z[:, out_of_bounds_x] = self.fill_value
354 if any_out_of_bounds_y:
355 z[out_of_bounds_y, :] = self.fill_value
357 if len(z) == 1:
358 z = z[0]
359 return array(z)
362def _check_broadcast_up_to(arr_from, shape_to, name):
363 """Helper to check that arr_from broadcasts up to shape_to"""
364 shape_from = arr_from.shape
365 if len(shape_to) >= len(shape_from):
366 for t, f in zip(shape_to[::-1], shape_from[::-1]):
367 if f != 1 and f != t:
368 break
369 else: # all checks pass, do the upcasting that we need later
370 if arr_from.size != 1 and arr_from.shape != shape_to:
371 arr_from = np.ones(shape_to, arr_from.dtype) * arr_from
372 return arr_from.ravel()
373 # at least one check failed
374 raise ValueError('%s argument must be able to broadcast up '
375 'to shape %s but had shape %s'
376 % (name, shape_to, shape_from))
379def _do_extrapolate(fill_value):
380 """Helper to check if fill_value == "extrapolate" without warnings"""
381 return (isinstance(fill_value, str) and
382 fill_value == 'extrapolate')
385class interp1d(_Interpolator1D):
386 """
387 Interpolate a 1-D function.
389 `x` and `y` are arrays of values used to approximate some function f:
390 ``y = f(x)``. This class returns a function whose call method uses
391 interpolation to find the value of new points.
393 Parameters
394 ----------
395 x : (N,) array_like
396 A 1-D array of real values.
397 y : (...,N,...) array_like
398 A N-D array of real values. The length of `y` along the interpolation
399 axis must be equal to the length of `x`.
400 kind : str or int, optional
401 Specifies the kind of interpolation as a string or as an integer
402 specifying the order of the spline interpolator to use.
403 The string has to be one of 'linear', 'nearest', 'nearest-up', 'zero',
404 'slinear', 'quadratic', 'cubic', 'previous', or 'next'. 'zero',
405 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation of
406 zeroth, first, second or third order; 'previous' and 'next' simply
407 return the previous or next value of the point; 'nearest-up' and
408 'nearest' differ when interpolating half-integers (e.g. 0.5, 1.5)
409 in that 'nearest-up' rounds up and 'nearest' rounds down. Default
410 is 'linear'.
411 axis : int, optional
412 Specifies the axis of `y` along which to interpolate.
413 Interpolation defaults to the last axis of `y`.
414 copy : bool, optional
415 If True, the class makes internal copies of x and y.
416 If False, references to `x` and `y` are used. The default is to copy.
417 bounds_error : bool, optional
418 If True, a ValueError is raised any time interpolation is attempted on
419 a value outside of the range of x (where extrapolation is
420 necessary). If False, out of bounds values are assigned `fill_value`.
421 By default, an error is raised unless ``fill_value="extrapolate"``.
422 fill_value : array-like or (array-like, array_like) or "extrapolate", optional
423 - if a ndarray (or float), this value will be used to fill in for
424 requested points outside of the data range. If not provided, then
425 the default is NaN. The array-like must broadcast properly to the
426 dimensions of the non-interpolation axes.
427 - If a two-element tuple, then the first element is used as a
428 fill value for ``x_new < x[0]`` and the second element is used for
429 ``x_new > x[-1]``. Anything that is not a 2-element tuple (e.g.,
430 list or ndarray, regardless of shape) is taken to be a single
431 array-like argument meant to be used for both bounds as
432 ``below, above = fill_value, fill_value``. Using a two-element tuple
433 or ndarray requires ``bounds_error=False``.
435 .. versionadded:: 0.17.0
436 - If "extrapolate", then points outside the data range will be
437 extrapolated.
439 .. versionadded:: 0.17.0
440 assume_sorted : bool, optional
441 If False, values of `x` can be in any order and they are sorted first.
442 If True, `x` has to be an array of monotonically increasing values.
444 Attributes
445 ----------
446 fill_value
448 Methods
449 -------
450 __call__
452 See Also
453 --------
454 splrep, splev
455 Spline interpolation/smoothing based on FITPACK.
456 UnivariateSpline : An object-oriented wrapper of the FITPACK routines.
457 interp2d : 2-D interpolation
459 Notes
460 -----
461 Calling `interp1d` with NaNs present in input values results in
462 undefined behaviour.
464 Input values `x` and `y` must be convertible to `float` values like
465 `int` or `float`.
467 If the values in `x` are not unique, the resulting behavior is
468 undefined and specific to the choice of `kind`, i.e., changing
469 `kind` will change the behavior for duplicates.
472 Examples
473 --------
474 >>> import numpy as np
475 >>> import matplotlib.pyplot as plt
476 >>> from scipy import interpolate
477 >>> x = np.arange(0, 10)
478 >>> y = np.exp(-x/3.0)
479 >>> f = interpolate.interp1d(x, y)
481 >>> xnew = np.arange(0, 9, 0.1)
482 >>> ynew = f(xnew) # use interpolation function returned by `interp1d`
483 >>> plt.plot(x, y, 'o', xnew, ynew, '-')
484 >>> plt.show()
485 """
487 def __init__(self, x, y, kind='linear', axis=-1,
488 copy=True, bounds_error=None, fill_value=np.nan,
489 assume_sorted=False):
490 """ Initialize a 1-D linear interpolation class."""
491 _Interpolator1D.__init__(self, x, y, axis=axis)
493 self.bounds_error = bounds_error # used by fill_value setter
494 self.copy = copy
496 if kind in ['zero', 'slinear', 'quadratic', 'cubic']:
497 order = {'zero': 0, 'slinear': 1,
498 'quadratic': 2, 'cubic': 3}[kind]
499 kind = 'spline'
500 elif isinstance(kind, int):
501 order = kind
502 kind = 'spline'
503 elif kind not in ('linear', 'nearest', 'nearest-up', 'previous',
504 'next'):
505 raise NotImplementedError("%s is unsupported: Use fitpack "
506 "routines for other types." % kind)
507 x = array(x, copy=self.copy)
508 y = array(y, copy=self.copy)
510 if not assume_sorted:
511 ind = np.argsort(x, kind="mergesort")
512 x = x[ind]
513 y = np.take(y, ind, axis=axis)
515 if x.ndim != 1:
516 raise ValueError("the x array must have exactly one dimension.")
517 if y.ndim == 0:
518 raise ValueError("the y array must have at least one dimension.")
520 # Force-cast y to a floating-point type, if it's not yet one
521 if not issubclass(y.dtype.type, np.inexact):
522 y = y.astype(np.float_)
524 # Backward compatibility
525 self.axis = axis % y.ndim
527 # Interpolation goes internally along the first axis
528 self.y = y
529 self._y = self._reshape_yi(self.y)
530 self.x = x
531 del y, x # clean up namespace to prevent misuse; use attributes
532 self._kind = kind
534 # Adjust to interpolation kind; store reference to *unbound*
535 # interpolation methods, in order to avoid circular references to self
536 # stored in the bound instance methods, and therefore delayed garbage
537 # collection. See: https://docs.python.org/reference/datamodel.html
538 if kind in ('linear', 'nearest', 'nearest-up', 'previous', 'next'):
539 # Make a "view" of the y array that is rotated to the interpolation
540 # axis.
541 minval = 1
542 if kind == 'nearest':
543 # Do division before addition to prevent possible integer
544 # overflow
545 self._side = 'left'
546 self.x_bds = self.x / 2.0
547 self.x_bds = self.x_bds[1:] + self.x_bds[:-1]
549 self._call = self.__class__._call_nearest
550 elif kind == 'nearest-up':
551 # Do division before addition to prevent possible integer
552 # overflow
553 self._side = 'right'
554 self.x_bds = self.x / 2.0
555 self.x_bds = self.x_bds[1:] + self.x_bds[:-1]
557 self._call = self.__class__._call_nearest
558 elif kind == 'previous':
559 # Side for np.searchsorted and index for clipping
560 self._side = 'left'
561 self._ind = 0
562 # Move x by one floating point value to the left
563 self._x_shift = np.nextafter(self.x, -np.inf)
564 self._call = self.__class__._call_previousnext
565 if _do_extrapolate(fill_value):
566 self._check_and_update_bounds_error_for_extrapolation()
567 # assume y is sorted by x ascending order here.
568 fill_value = (np.nan, np.take(self.y, -1, axis))
569 elif kind == 'next':
570 self._side = 'right'
571 self._ind = 1
572 # Move x by one floating point value to the right
573 self._x_shift = np.nextafter(self.x, np.inf)
574 self._call = self.__class__._call_previousnext
575 if _do_extrapolate(fill_value):
576 self._check_and_update_bounds_error_for_extrapolation()
577 # assume y is sorted by x ascending order here.
578 fill_value = (np.take(self.y, 0, axis), np.nan)
579 else:
580 # Check if we can delegate to numpy.interp (2x-10x faster).
581 np_types = (np.float_, np.int_)
582 cond = self.x.dtype in np_types and self.y.dtype in np_types
583 cond = cond and self.y.ndim == 1
584 cond = cond and not _do_extrapolate(fill_value)
586 if cond:
587 self._call = self.__class__._call_linear_np
588 else:
589 self._call = self.__class__._call_linear
590 else:
591 minval = order + 1
593 rewrite_nan = False
594 xx, yy = self.x, self._y
595 if order > 1:
596 # Quadratic or cubic spline. If input contains even a single
597 # nan, then the output is all nans. We cannot just feed data
598 # with nans to make_interp_spline because it calls LAPACK.
599 # So, we make up a bogus x and y with no nans and use it
600 # to get the correct shape of the output, which we then fill
601 # with nans.
602 # For slinear or zero order spline, we just pass nans through.
603 mask = np.isnan(self.x)
604 if mask.any():
605 sx = self.x[~mask]
606 if sx.size == 0:
607 raise ValueError("`x` array is all-nan")
608 xx = np.linspace(np.nanmin(self.x),
609 np.nanmax(self.x),
610 len(self.x))
611 rewrite_nan = True
612 if np.isnan(self._y).any():
613 yy = np.ones_like(self._y)
614 rewrite_nan = True
616 self._spline = make_interp_spline(xx, yy, k=order,
617 check_finite=False)
618 if rewrite_nan:
619 self._call = self.__class__._call_nan_spline
620 else:
621 self._call = self.__class__._call_spline
623 if len(self.x) < minval:
624 raise ValueError("x and y arrays must have at "
625 "least %d entries" % minval)
627 self.fill_value = fill_value # calls the setter, can modify bounds_err
629 @property
630 def fill_value(self):
631 """The fill value."""
632 # backwards compat: mimic a public attribute
633 return self._fill_value_orig
635 @fill_value.setter
636 def fill_value(self, fill_value):
637 # extrapolation only works for nearest neighbor and linear methods
638 if _do_extrapolate(fill_value):
639 self._check_and_update_bounds_error_for_extrapolation()
640 self._extrapolate = True
641 else:
642 broadcast_shape = (self.y.shape[:self.axis] +
643 self.y.shape[self.axis + 1:])
644 if len(broadcast_shape) == 0:
645 broadcast_shape = (1,)
646 # it's either a pair (_below_range, _above_range) or a single value
647 # for both above and below range
648 if isinstance(fill_value, tuple) and len(fill_value) == 2:
649 below_above = [np.asarray(fill_value[0]),
650 np.asarray(fill_value[1])]
651 names = ('fill_value (below)', 'fill_value (above)')
652 for ii in range(2):
653 below_above[ii] = _check_broadcast_up_to(
654 below_above[ii], broadcast_shape, names[ii])
655 else:
656 fill_value = np.asarray(fill_value)
657 below_above = [_check_broadcast_up_to(
658 fill_value, broadcast_shape, 'fill_value')] * 2
659 self._fill_value_below, self._fill_value_above = below_above
660 self._extrapolate = False
661 if self.bounds_error is None:
662 self.bounds_error = True
663 # backwards compat: fill_value was a public attr; make it writeable
664 self._fill_value_orig = fill_value
666 def _check_and_update_bounds_error_for_extrapolation(self):
667 if self.bounds_error:
668 raise ValueError("Cannot extrapolate and raise "
669 "at the same time.")
670 self.bounds_error = False
672 def _call_linear_np(self, x_new):
673 # Note that out-of-bounds values are taken care of in self._evaluate
674 return np.interp(x_new, self.x, self.y)
676 def _call_linear(self, x_new):
677 # 2. Find where in the original data, the values to interpolate
678 # would be inserted.
679 # Note: If x_new[n] == x[m], then m is returned by searchsorted.
680 x_new_indices = searchsorted(self.x, x_new)
682 # 3. Clip x_new_indices so that they are within the range of
683 # self.x indices and at least 1. Removes mis-interpolation
684 # of x_new[n] = x[0]
685 x_new_indices = x_new_indices.clip(1, len(self.x)-1).astype(int)
687 # 4. Calculate the slope of regions that each x_new value falls in.
688 lo = x_new_indices - 1
689 hi = x_new_indices
691 x_lo = self.x[lo]
692 x_hi = self.x[hi]
693 y_lo = self._y[lo]
694 y_hi = self._y[hi]
696 # Note that the following two expressions rely on the specifics of the
697 # broadcasting semantics.
698 slope = (y_hi - y_lo) / (x_hi - x_lo)[:, None]
700 # 5. Calculate the actual value for each entry in x_new.
701 y_new = slope*(x_new - x_lo)[:, None] + y_lo
703 return y_new
705 def _call_nearest(self, x_new):
706 """ Find nearest neighbor interpolated y_new = f(x_new)."""
708 # 2. Find where in the averaged data the values to interpolate
709 # would be inserted.
710 # Note: use side='left' (right) to searchsorted() to define the
711 # halfway point to be nearest to the left (right) neighbor
712 x_new_indices = searchsorted(self.x_bds, x_new, side=self._side)
714 # 3. Clip x_new_indices so that they are within the range of x indices.
715 x_new_indices = x_new_indices.clip(0, len(self.x)-1).astype(intp)
717 # 4. Calculate the actual value for each entry in x_new.
718 y_new = self._y[x_new_indices]
720 return y_new
722 def _call_previousnext(self, x_new):
723 """Use previous/next neighbor of x_new, y_new = f(x_new)."""
725 # 1. Get index of left/right value
726 x_new_indices = searchsorted(self._x_shift, x_new, side=self._side)
728 # 2. Clip x_new_indices so that they are within the range of x indices.
729 x_new_indices = x_new_indices.clip(1-self._ind,
730 len(self.x)-self._ind).astype(intp)
732 # 3. Calculate the actual value for each entry in x_new.
733 y_new = self._y[x_new_indices+self._ind-1]
735 return y_new
737 def _call_spline(self, x_new):
738 return self._spline(x_new)
740 def _call_nan_spline(self, x_new):
741 out = self._spline(x_new)
742 out[...] = np.nan
743 return out
745 def _evaluate(self, x_new):
746 # 1. Handle values in x_new that are outside of x. Throw error,
747 # or return a list of mask array indicating the outofbounds values.
748 # The behavior is set by the bounds_error variable.
749 x_new = asarray(x_new)
750 y_new = self._call(self, x_new)
751 if not self._extrapolate:
752 below_bounds, above_bounds = self._check_bounds(x_new)
753 if len(y_new) > 0:
754 # Note fill_value must be broadcast up to the proper size
755 # and flattened to work here
756 y_new[below_bounds] = self._fill_value_below
757 y_new[above_bounds] = self._fill_value_above
758 return y_new
760 def _check_bounds(self, x_new):
761 """Check the inputs for being in the bounds of the interpolated data.
763 Parameters
764 ----------
765 x_new : array
767 Returns
768 -------
769 out_of_bounds : bool array
770 The mask on x_new of values that are out of the bounds.
771 """
773 # If self.bounds_error is True, we raise an error if any x_new values
774 # fall outside the range of x. Otherwise, we return an array indicating
775 # which values are outside the boundary region.
776 below_bounds = x_new < self.x[0]
777 above_bounds = x_new > self.x[-1]
779 if self.bounds_error and below_bounds.any():
780 below_bounds_value = x_new[np.argmax(below_bounds)]
781 raise ValueError("A value ({}) in x_new is below "
782 "the interpolation range's minimum value ({})."
783 .format(below_bounds_value, self.x[0]))
784 if self.bounds_error and above_bounds.any():
785 above_bounds_value = x_new[np.argmax(above_bounds)]
786 raise ValueError("A value ({}) in x_new is above "
787 "the interpolation range's maximum value ({})."
788 .format(above_bounds_value, self.x[-1]))
790 # !! Should we emit a warning if some values are out of bounds?
791 # !! matlab does not.
792 return below_bounds, above_bounds
795class _PPolyBase:
796 """Base class for piecewise polynomials."""
797 __slots__ = ('c', 'x', 'extrapolate', 'axis')
799 def __init__(self, c, x, extrapolate=None, axis=0):
800 self.c = np.asarray(c)
801 self.x = np.ascontiguousarray(x, dtype=np.float64)
803 if extrapolate is None:
804 extrapolate = True
805 elif extrapolate != 'periodic':
806 extrapolate = bool(extrapolate)
807 self.extrapolate = extrapolate
809 if self.c.ndim < 2:
810 raise ValueError("Coefficients array must be at least "
811 "2-dimensional.")
813 if not (0 <= axis < self.c.ndim - 1):
814 raise ValueError("axis=%s must be between 0 and %s" %
815 (axis, self.c.ndim-1))
817 self.axis = axis
818 if axis != 0:
819 # move the interpolation axis to be the first one in self.c
820 # More specifically, the target shape for self.c is (k, m, ...),
821 # and axis !=0 means that we have c.shape (..., k, m, ...)
822 # ^
823 # axis
824 # So we roll two of them.
825 self.c = np.moveaxis(self.c, axis+1, 0)
826 self.c = np.moveaxis(self.c, axis+1, 0)
828 if self.x.ndim != 1:
829 raise ValueError("x must be 1-dimensional")
830 if self.x.size < 2:
831 raise ValueError("at least 2 breakpoints are needed")
832 if self.c.ndim < 2:
833 raise ValueError("c must have at least 2 dimensions")
834 if self.c.shape[0] == 0:
835 raise ValueError("polynomial must be at least of order 0")
836 if self.c.shape[1] != self.x.size-1:
837 raise ValueError("number of coefficients != len(x)-1")
838 dx = np.diff(self.x)
839 if not (np.all(dx >= 0) or np.all(dx <= 0)):
840 raise ValueError("`x` must be strictly increasing or decreasing.")
842 dtype = self._get_dtype(self.c.dtype)
843 self.c = np.ascontiguousarray(self.c, dtype=dtype)
845 def _get_dtype(self, dtype):
846 if np.issubdtype(dtype, np.complexfloating) \
847 or np.issubdtype(self.c.dtype, np.complexfloating):
848 return np.complex_
849 else:
850 return np.float_
852 @classmethod
853 def construct_fast(cls, c, x, extrapolate=None, axis=0):
854 """
855 Construct the piecewise polynomial without making checks.
857 Takes the same parameters as the constructor. Input arguments
858 ``c`` and ``x`` must be arrays of the correct shape and type. The
859 ``c`` array can only be of dtypes float and complex, and ``x``
860 array must have dtype float.
861 """
862 self = object.__new__(cls)
863 self.c = c
864 self.x = x
865 self.axis = axis
866 if extrapolate is None:
867 extrapolate = True
868 self.extrapolate = extrapolate
869 return self
871 def _ensure_c_contiguous(self):
872 """
873 c and x may be modified by the user. The Cython code expects
874 that they are C contiguous.
875 """
876 if not self.x.flags.c_contiguous:
877 self.x = self.x.copy()
878 if not self.c.flags.c_contiguous:
879 self.c = self.c.copy()
881 def extend(self, c, x):
882 """
883 Add additional breakpoints and coefficients to the polynomial.
885 Parameters
886 ----------
887 c : ndarray, size (k, m, ...)
888 Additional coefficients for polynomials in intervals. Note that
889 the first additional interval will be formed using one of the
890 ``self.x`` end points.
891 x : ndarray, size (m,)
892 Additional breakpoints. Must be sorted in the same order as
893 ``self.x`` and either to the right or to the left of the current
894 breakpoints.
895 """
897 c = np.asarray(c)
898 x = np.asarray(x)
900 if c.ndim < 2:
901 raise ValueError("invalid dimensions for c")
902 if x.ndim != 1:
903 raise ValueError("invalid dimensions for x")
904 if x.shape[0] != c.shape[1]:
905 raise ValueError("Shapes of x {} and c {} are incompatible"
906 .format(x.shape, c.shape))
907 if c.shape[2:] != self.c.shape[2:] or c.ndim != self.c.ndim:
908 raise ValueError("Shapes of c {} and self.c {} are incompatible"
909 .format(c.shape, self.c.shape))
911 if c.size == 0:
912 return
914 dx = np.diff(x)
915 if not (np.all(dx >= 0) or np.all(dx <= 0)):
916 raise ValueError("`x` is not sorted.")
918 if self.x[-1] >= self.x[0]:
919 if not x[-1] >= x[0]:
920 raise ValueError("`x` is in the different order "
921 "than `self.x`.")
923 if x[0] >= self.x[-1]:
924 action = 'append'
925 elif x[-1] <= self.x[0]:
926 action = 'prepend'
927 else:
928 raise ValueError("`x` is neither on the left or on the right "
929 "from `self.x`.")
930 else:
931 if not x[-1] <= x[0]:
932 raise ValueError("`x` is in the different order "
933 "than `self.x`.")
935 if x[0] <= self.x[-1]:
936 action = 'append'
937 elif x[-1] >= self.x[0]:
938 action = 'prepend'
939 else:
940 raise ValueError("`x` is neither on the left or on the right "
941 "from `self.x`.")
943 dtype = self._get_dtype(c.dtype)
945 k2 = max(c.shape[0], self.c.shape[0])
946 c2 = np.zeros((k2, self.c.shape[1] + c.shape[1]) + self.c.shape[2:],
947 dtype=dtype)
949 if action == 'append':
950 c2[k2-self.c.shape[0]:, :self.c.shape[1]] = self.c
951 c2[k2-c.shape[0]:, self.c.shape[1]:] = c
952 self.x = np.r_[self.x, x]
953 elif action == 'prepend':
954 c2[k2-self.c.shape[0]:, :c.shape[1]] = c
955 c2[k2-c.shape[0]:, c.shape[1]:] = self.c
956 self.x = np.r_[x, self.x]
958 self.c = c2
960 def __call__(self, x, nu=0, extrapolate=None):
961 """
962 Evaluate the piecewise polynomial or its derivative.
964 Parameters
965 ----------
966 x : array_like
967 Points to evaluate the interpolant at.
968 nu : int, optional
969 Order of derivative to evaluate. Must be non-negative.
970 extrapolate : {bool, 'periodic', None}, optional
971 If bool, determines whether to extrapolate to out-of-bounds points
972 based on first and last intervals, or to return NaNs.
973 If 'periodic', periodic extrapolation is used.
974 If None (default), use `self.extrapolate`.
976 Returns
977 -------
978 y : array_like
979 Interpolated values. Shape is determined by replacing
980 the interpolation axis in the original array with the shape of x.
982 Notes
983 -----
984 Derivatives are evaluated piecewise for each polynomial
985 segment, even if the polynomial is not differentiable at the
986 breakpoints. The polynomial intervals are considered half-open,
987 ``[a, b)``, except for the last interval which is closed
988 ``[a, b]``.
989 """
990 if extrapolate is None:
991 extrapolate = self.extrapolate
992 x = np.asarray(x)
993 x_shape, x_ndim = x.shape, x.ndim
994 x = np.ascontiguousarray(x.ravel(), dtype=np.float_)
996 # With periodic extrapolation we map x to the segment
997 # [self.x[0], self.x[-1]].
998 if extrapolate == 'periodic':
999 x = self.x[0] + (x - self.x[0]) % (self.x[-1] - self.x[0])
1000 extrapolate = False
1002 out = np.empty((len(x), prod(self.c.shape[2:])), dtype=self.c.dtype)
1003 self._ensure_c_contiguous()
1004 self._evaluate(x, nu, extrapolate, out)
1005 out = out.reshape(x_shape + self.c.shape[2:])
1006 if self.axis != 0:
1007 # transpose to move the calculated values to the interpolation axis
1008 l = list(range(out.ndim))
1009 l = l[x_ndim:x_ndim+self.axis] + l[:x_ndim] + l[x_ndim+self.axis:]
1010 out = out.transpose(l)
1011 return out
1014class PPoly(_PPolyBase):
1015 """
1016 Piecewise polynomial in terms of coefficients and breakpoints
1018 The polynomial between ``x[i]`` and ``x[i + 1]`` is written in the
1019 local power basis::
1021 S = sum(c[m, i] * (xp - x[i])**(k-m) for m in range(k+1))
1023 where ``k`` is the degree of the polynomial.
1025 Parameters
1026 ----------
1027 c : ndarray, shape (k, m, ...)
1028 Polynomial coefficients, order `k` and `m` intervals.
1029 x : ndarray, shape (m+1,)
1030 Polynomial breakpoints. Must be sorted in either increasing or
1031 decreasing order.
1032 extrapolate : bool or 'periodic', optional
1033 If bool, determines whether to extrapolate to out-of-bounds points
1034 based on first and last intervals, or to return NaNs. If 'periodic',
1035 periodic extrapolation is used. Default is True.
1036 axis : int, optional
1037 Interpolation axis. Default is zero.
1039 Attributes
1040 ----------
1041 x : ndarray
1042 Breakpoints.
1043 c : ndarray
1044 Coefficients of the polynomials. They are reshaped
1045 to a 3-D array with the last dimension representing
1046 the trailing dimensions of the original coefficient array.
1047 axis : int
1048 Interpolation axis.
1050 Methods
1051 -------
1052 __call__
1053 derivative
1054 antiderivative
1055 integrate
1056 solve
1057 roots
1058 extend
1059 from_spline
1060 from_bernstein_basis
1061 construct_fast
1063 See also
1064 --------
1065 BPoly : piecewise polynomials in the Bernstein basis
1067 Notes
1068 -----
1069 High-order polynomials in the power basis can be numerically
1070 unstable. Precision problems can start to appear for orders
1071 larger than 20-30.
1072 """
1074 def _evaluate(self, x, nu, extrapolate, out):
1075 _ppoly.evaluate(self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1076 self.x, x, nu, bool(extrapolate), out)
1078 def derivative(self, nu=1):
1079 """
1080 Construct a new piecewise polynomial representing the derivative.
1082 Parameters
1083 ----------
1084 nu : int, optional
1085 Order of derivative to evaluate. Default is 1, i.e., compute the
1086 first derivative. If negative, the antiderivative is returned.
1088 Returns
1089 -------
1090 pp : PPoly
1091 Piecewise polynomial of order k2 = k - n representing the derivative
1092 of this polynomial.
1094 Notes
1095 -----
1096 Derivatives are evaluated piecewise for each polynomial
1097 segment, even if the polynomial is not differentiable at the
1098 breakpoints. The polynomial intervals are considered half-open,
1099 ``[a, b)``, except for the last interval which is closed
1100 ``[a, b]``.
1101 """
1102 if nu < 0:
1103 return self.antiderivative(-nu)
1105 # reduce order
1106 if nu == 0:
1107 c2 = self.c.copy()
1108 else:
1109 c2 = self.c[:-nu, :].copy()
1111 if c2.shape[0] == 0:
1112 # derivative of order 0 is zero
1113 c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)
1115 # multiply by the correct rising factorials
1116 factor = spec.poch(np.arange(c2.shape[0], 0, -1), nu)
1117 c2 *= factor[(slice(None),) + (None,)*(c2.ndim-1)]
1119 # construct a compatible polynomial
1120 return self.construct_fast(c2, self.x, self.extrapolate, self.axis)
1122 def antiderivative(self, nu=1):
1123 """
1124 Construct a new piecewise polynomial representing the antiderivative.
1126 Antiderivative is also the indefinite integral of the function,
1127 and derivative is its inverse operation.
1129 Parameters
1130 ----------
1131 nu : int, optional
1132 Order of antiderivative to evaluate. Default is 1, i.e., compute
1133 the first integral. If negative, the derivative is returned.
1135 Returns
1136 -------
1137 pp : PPoly
1138 Piecewise polynomial of order k2 = k + n representing
1139 the antiderivative of this polynomial.
1141 Notes
1142 -----
1143 The antiderivative returned by this function is continuous and
1144 continuously differentiable to order n-1, up to floating point
1145 rounding error.
1147 If antiderivative is computed and ``self.extrapolate='periodic'``,
1148 it will be set to False for the returned instance. This is done because
1149 the antiderivative is no longer periodic and its correct evaluation
1150 outside of the initially given x interval is difficult.
1151 """
1152 if nu <= 0:
1153 return self.derivative(-nu)
1155 c = np.zeros((self.c.shape[0] + nu, self.c.shape[1]) + self.c.shape[2:],
1156 dtype=self.c.dtype)
1157 c[:-nu] = self.c
1159 # divide by the correct rising factorials
1160 factor = spec.poch(np.arange(self.c.shape[0], 0, -1), nu)
1161 c[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]
1163 # fix continuity of added degrees of freedom
1164 self._ensure_c_contiguous()
1165 _ppoly.fix_continuity(c.reshape(c.shape[0], c.shape[1], -1),
1166 self.x, nu - 1)
1168 if self.extrapolate == 'periodic':
1169 extrapolate = False
1170 else:
1171 extrapolate = self.extrapolate
1173 # construct a compatible polynomial
1174 return self.construct_fast(c, self.x, extrapolate, self.axis)
1176 def integrate(self, a, b, extrapolate=None):
1177 """
1178 Compute a definite integral over a piecewise polynomial.
1180 Parameters
1181 ----------
1182 a : float
1183 Lower integration bound
1184 b : float
1185 Upper integration bound
1186 extrapolate : {bool, 'periodic', None}, optional
1187 If bool, determines whether to extrapolate to out-of-bounds points
1188 based on first and last intervals, or to return NaNs.
1189 If 'periodic', periodic extrapolation is used.
1190 If None (default), use `self.extrapolate`.
1192 Returns
1193 -------
1194 ig : array_like
1195 Definite integral of the piecewise polynomial over [a, b]
1196 """
1197 if extrapolate is None:
1198 extrapolate = self.extrapolate
1200 # Swap integration bounds if needed
1201 sign = 1
1202 if b < a:
1203 a, b = b, a
1204 sign = -1
1206 range_int = np.empty((prod(self.c.shape[2:]),), dtype=self.c.dtype)
1207 self._ensure_c_contiguous()
1209 # Compute the integral.
1210 if extrapolate == 'periodic':
1211 # Split the integral into the part over period (can be several
1212 # of them) and the remaining part.
1214 xs, xe = self.x[0], self.x[-1]
1215 period = xe - xs
1216 interval = b - a
1217 n_periods, left = divmod(interval, period)
1219 if n_periods > 0:
1220 _ppoly.integrate(
1221 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1222 self.x, xs, xe, False, out=range_int)
1223 range_int *= n_periods
1224 else:
1225 range_int.fill(0)
1227 # Map a to [xs, xe], b is always a + left.
1228 a = xs + (a - xs) % period
1229 b = a + left
1231 # If b <= xe then we need to integrate over [a, b], otherwise
1232 # over [a, xe] and from xs to what is remained.
1233 remainder_int = np.empty_like(range_int)
1234 if b <= xe:
1235 _ppoly.integrate(
1236 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1237 self.x, a, b, False, out=remainder_int)
1238 range_int += remainder_int
1239 else:
1240 _ppoly.integrate(
1241 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1242 self.x, a, xe, False, out=remainder_int)
1243 range_int += remainder_int
1245 _ppoly.integrate(
1246 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1247 self.x, xs, xs + left + a - xe, False, out=remainder_int)
1248 range_int += remainder_int
1249 else:
1250 _ppoly.integrate(
1251 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1252 self.x, a, b, bool(extrapolate), out=range_int)
1254 # Return
1255 range_int *= sign
1256 return range_int.reshape(self.c.shape[2:])
1258 def solve(self, y=0., discontinuity=True, extrapolate=None):
1259 """
1260 Find real solutions of the equation ``pp(x) == y``.
1262 Parameters
1263 ----------
1264 y : float, optional
1265 Right-hand side. Default is zero.
1266 discontinuity : bool, optional
1267 Whether to report sign changes across discontinuities at
1268 breakpoints as roots.
1269 extrapolate : {bool, 'periodic', None}, optional
1270 If bool, determines whether to return roots from the polynomial
1271 extrapolated based on first and last intervals, 'periodic' works
1272 the same as False. If None (default), use `self.extrapolate`.
1274 Returns
1275 -------
1276 roots : ndarray
1277 Roots of the polynomial(s).
1279 If the PPoly object describes multiple polynomials, the
1280 return value is an object array whose each element is an
1281 ndarray containing the roots.
1283 Notes
1284 -----
1285 This routine works only on real-valued polynomials.
1287 If the piecewise polynomial contains sections that are
1288 identically zero, the root list will contain the start point
1289 of the corresponding interval, followed by a ``nan`` value.
1291 If the polynomial is discontinuous across a breakpoint, and
1292 there is a sign change across the breakpoint, this is reported
1293 if the `discont` parameter is True.
1295 Examples
1296 --------
1298 Finding roots of ``[x**2 - 1, (x - 1)**2]`` defined on intervals
1299 ``[-2, 1], [1, 2]``:
1301 >>> import numpy as np
1302 >>> from scipy.interpolate import PPoly
1303 >>> pp = PPoly(np.array([[1, -4, 3], [1, 0, 0]]).T, [-2, 1, 2])
1304 >>> pp.solve()
1305 array([-1., 1.])
1306 """
1307 if extrapolate is None:
1308 extrapolate = self.extrapolate
1310 self._ensure_c_contiguous()
1312 if np.issubdtype(self.c.dtype, np.complexfloating):
1313 raise ValueError("Root finding is only for "
1314 "real-valued polynomials")
1316 y = float(y)
1317 r = _ppoly.real_roots(self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1318 self.x, y, bool(discontinuity),
1319 bool(extrapolate))
1320 if self.c.ndim == 2:
1321 return r[0]
1322 else:
1323 r2 = np.empty(prod(self.c.shape[2:]), dtype=object)
1324 # this for-loop is equivalent to ``r2[...] = r``, but that's broken
1325 # in NumPy 1.6.0
1326 for ii, root in enumerate(r):
1327 r2[ii] = root
1329 return r2.reshape(self.c.shape[2:])
1331 def roots(self, discontinuity=True, extrapolate=None):
1332 """
1333 Find real roots of the piecewise polynomial.
1335 Parameters
1336 ----------
1337 discontinuity : bool, optional
1338 Whether to report sign changes across discontinuities at
1339 breakpoints as roots.
1340 extrapolate : {bool, 'periodic', None}, optional
1341 If bool, determines whether to return roots from the polynomial
1342 extrapolated based on first and last intervals, 'periodic' works
1343 the same as False. If None (default), use `self.extrapolate`.
1345 Returns
1346 -------
1347 roots : ndarray
1348 Roots of the polynomial(s).
1350 If the PPoly object describes multiple polynomials, the
1351 return value is an object array whose each element is an
1352 ndarray containing the roots.
1354 See Also
1355 --------
1356 PPoly.solve
1357 """
1358 return self.solve(0, discontinuity, extrapolate)
1360 @classmethod
1361 def from_spline(cls, tck, extrapolate=None):
1362 """
1363 Construct a piecewise polynomial from a spline
1365 Parameters
1366 ----------
1367 tck
1368 A spline, as returned by `splrep` or a BSpline object.
1369 extrapolate : bool or 'periodic', optional
1370 If bool, determines whether to extrapolate to out-of-bounds points
1371 based on first and last intervals, or to return NaNs.
1372 If 'periodic', periodic extrapolation is used. Default is True.
1374 Examples
1375 --------
1376 Construct an interpolating spline and convert it to a `PPoly` instance
1378 >>> import numpy as np
1379 >>> from scipy.interpolate import splrep, PPoly
1380 >>> x = np.linspace(0, 1, 11)
1381 >>> y = np.sin(2*np.pi*x)
1382 >>> tck = splrep(x, y, s=0)
1383 >>> p = PPoly.from_spline(tck)
1384 >>> isinstance(p, PPoly)
1385 True
1387 Note that this function only supports 1D splines out of the box.
1389 If the ``tck`` object represents a parametric spline (e.g. constructed
1390 by `splprep` or a `BSpline` with ``c.ndim > 1``), you will need to loop
1391 over the dimensions manually.
1393 >>> from scipy.interpolate import splprep, splev
1394 >>> t = np.linspace(0, 1, 11)
1395 >>> x = np.sin(2*np.pi*t)
1396 >>> y = np.cos(2*np.pi*t)
1397 >>> (t, c, k), u = splprep([x, y], s=0)
1399 Note that ``c`` is a list of two arrays of length 11.
1401 >>> unew = np.arange(0, 1.01, 0.01)
1402 >>> out = splev(unew, (t, c, k))
1404 To convert this spline to the power basis, we convert each
1405 component of the list of b-spline coefficients, ``c``, into the
1406 corresponding cubic polynomial.
1408 >>> polys = [PPoly.from_spline((t, cj, k)) for cj in c]
1409 >>> polys[0].c.shape
1410 (4, 14)
1412 Note that the coefficients of the polynomials `polys` are in the
1413 power basis and their dimensions reflect just that: here 4 is the order
1414 (degree+1), and 14 is the number of intervals---which is nothing but
1415 the length of the knot array of the original `tck` minus one.
1417 Optionally, we can stack the components into a single `PPoly` along
1418 the third dimension:
1420 >>> cc = np.dstack([p.c for p in polys]) # has shape = (4, 14, 2)
1421 >>> poly = PPoly(cc, polys[0].x)
1422 >>> np.allclose(poly(unew).T, # note the transpose to match `splev`
1423 ... out, atol=1e-15)
1424 True
1426 """
1427 if isinstance(tck, BSpline):
1428 t, c, k = tck.tck
1429 if extrapolate is None:
1430 extrapolate = tck.extrapolate
1431 else:
1432 t, c, k = tck
1434 cvals = np.empty((k + 1, len(t)-1), dtype=c.dtype)
1435 for m in range(k, -1, -1):
1436 y = _fitpack_py.splev(t[:-1], tck, der=m)
1437 cvals[k - m, :] = y/spec.gamma(m+1)
1439 return cls.construct_fast(cvals, t, extrapolate)
1441 @classmethod
1442 def from_bernstein_basis(cls, bp, extrapolate=None):
1443 """
1444 Construct a piecewise polynomial in the power basis
1445 from a polynomial in Bernstein basis.
1447 Parameters
1448 ----------
1449 bp : BPoly
1450 A Bernstein basis polynomial, as created by BPoly
1451 extrapolate : bool or 'periodic', optional
1452 If bool, determines whether to extrapolate to out-of-bounds points
1453 based on first and last intervals, or to return NaNs.
1454 If 'periodic', periodic extrapolation is used. Default is True.
1455 """
1456 if not isinstance(bp, BPoly):
1457 raise TypeError(".from_bernstein_basis only accepts BPoly instances. "
1458 "Got %s instead." % type(bp))
1460 dx = np.diff(bp.x)
1461 k = bp.c.shape[0] - 1 # polynomial order
1463 rest = (None,)*(bp.c.ndim-2)
1465 c = np.zeros_like(bp.c)
1466 for a in range(k+1):
1467 factor = (-1)**a * comb(k, a) * bp.c[a]
1468 for s in range(a, k+1):
1469 val = comb(k-a, s-a) * (-1)**s
1470 c[k-s] += factor * val / dx[(slice(None),)+rest]**s
1472 if extrapolate is None:
1473 extrapolate = bp.extrapolate
1475 return cls.construct_fast(c, bp.x, extrapolate, bp.axis)
1478class BPoly(_PPolyBase):
1479 """Piecewise polynomial in terms of coefficients and breakpoints.
1481 The polynomial between ``x[i]`` and ``x[i + 1]`` is written in the
1482 Bernstein polynomial basis::
1484 S = sum(c[a, i] * b(a, k; x) for a in range(k+1)),
1486 where ``k`` is the degree of the polynomial, and::
1488 b(a, k; x) = binom(k, a) * t**a * (1 - t)**(k - a),
1490 with ``t = (x - x[i]) / (x[i+1] - x[i])`` and ``binom`` is the binomial
1491 coefficient.
1493 Parameters
1494 ----------
1495 c : ndarray, shape (k, m, ...)
1496 Polynomial coefficients, order `k` and `m` intervals
1497 x : ndarray, shape (m+1,)
1498 Polynomial breakpoints. Must be sorted in either increasing or
1499 decreasing order.
1500 extrapolate : bool, optional
1501 If bool, determines whether to extrapolate to out-of-bounds points
1502 based on first and last intervals, or to return NaNs. If 'periodic',
1503 periodic extrapolation is used. Default is True.
1504 axis : int, optional
1505 Interpolation axis. Default is zero.
1507 Attributes
1508 ----------
1509 x : ndarray
1510 Breakpoints.
1511 c : ndarray
1512 Coefficients of the polynomials. They are reshaped
1513 to a 3-D array with the last dimension representing
1514 the trailing dimensions of the original coefficient array.
1515 axis : int
1516 Interpolation axis.
1518 Methods
1519 -------
1520 __call__
1521 extend
1522 derivative
1523 antiderivative
1524 integrate
1525 construct_fast
1526 from_power_basis
1527 from_derivatives
1529 See also
1530 --------
1531 PPoly : piecewise polynomials in the power basis
1533 Notes
1534 -----
1535 Properties of Bernstein polynomials are well documented in the literature,
1536 see for example [1]_ [2]_ [3]_.
1538 References
1539 ----------
1540 .. [1] https://en.wikipedia.org/wiki/Bernstein_polynomial
1542 .. [2] Kenneth I. Joy, Bernstein polynomials,
1543 http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf
1545 .. [3] E. H. Doha, A. H. Bhrawy, and M. A. Saker, Boundary Value Problems,
1546 vol 2011, article ID 829546, :doi:`10.1155/2011/829543`.
1548 Examples
1549 --------
1550 >>> from scipy.interpolate import BPoly
1551 >>> x = [0, 1]
1552 >>> c = [[1], [2], [3]]
1553 >>> bp = BPoly(c, x)
1555 This creates a 2nd order polynomial
1557 .. math::
1559 B(x) = 1 \\times b_{0, 2}(x) + 2 \\times b_{1, 2}(x) + 3 \\times b_{2, 2}(x) \\\\
1560 = 1 \\times (1-x)^2 + 2 \\times 2 x (1 - x) + 3 \\times x^2
1562 """
1564 def _evaluate(self, x, nu, extrapolate, out):
1565 _ppoly.evaluate_bernstein(
1566 self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
1567 self.x, x, nu, bool(extrapolate), out)
1569 def derivative(self, nu=1):
1570 """
1571 Construct a new piecewise polynomial representing the derivative.
1573 Parameters
1574 ----------
1575 nu : int, optional
1576 Order of derivative to evaluate. Default is 1, i.e., compute the
1577 first derivative. If negative, the antiderivative is returned.
1579 Returns
1580 -------
1581 bp : BPoly
1582 Piecewise polynomial of order k - nu representing the derivative of
1583 this polynomial.
1585 """
1586 if nu < 0:
1587 return self.antiderivative(-nu)
1589 if nu > 1:
1590 bp = self
1591 for k in range(nu):
1592 bp = bp.derivative()
1593 return bp
1595 # reduce order
1596 if nu == 0:
1597 c2 = self.c.copy()
1598 else:
1599 # For a polynomial
1600 # B(x) = \sum_{a=0}^{k} c_a b_{a, k}(x),
1601 # we use the fact that
1602 # b'_{a, k} = k ( b_{a-1, k-1} - b_{a, k-1} ),
1603 # which leads to
1604 # B'(x) = \sum_{a=0}^{k-1} (c_{a+1} - c_a) b_{a, k-1}
1605 #
1606 # finally, for an interval [y, y + dy] with dy != 1,
1607 # we need to correct for an extra power of dy
1609 rest = (None,)*(self.c.ndim-2)
1611 k = self.c.shape[0] - 1
1612 dx = np.diff(self.x)[(None, slice(None))+rest]
1613 c2 = k * np.diff(self.c, axis=0) / dx
1615 if c2.shape[0] == 0:
1616 # derivative of order 0 is zero
1617 c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)
1619 # construct a compatible polynomial
1620 return self.construct_fast(c2, self.x, self.extrapolate, self.axis)
1622 def antiderivative(self, nu=1):
1623 """
1624 Construct a new piecewise polynomial representing the antiderivative.
1626 Parameters
1627 ----------
1628 nu : int, optional
1629 Order of antiderivative to evaluate. Default is 1, i.e., compute
1630 the first integral. If negative, the derivative is returned.
1632 Returns
1633 -------
1634 bp : BPoly
1635 Piecewise polynomial of order k + nu representing the
1636 antiderivative of this polynomial.
1638 Notes
1639 -----
1640 If antiderivative is computed and ``self.extrapolate='periodic'``,
1641 it will be set to False for the returned instance. This is done because
1642 the antiderivative is no longer periodic and its correct evaluation
1643 outside of the initially given x interval is difficult.
1644 """
1645 if nu <= 0:
1646 return self.derivative(-nu)
1648 if nu > 1:
1649 bp = self
1650 for k in range(nu):
1651 bp = bp.antiderivative()
1652 return bp
1654 # Construct the indefinite integrals on individual intervals
1655 c, x = self.c, self.x
1656 k = c.shape[0]
1657 c2 = np.zeros((k+1,) + c.shape[1:], dtype=c.dtype)
1659 c2[1:, ...] = np.cumsum(c, axis=0) / k
1660 delta = x[1:] - x[:-1]
1661 c2 *= delta[(None, slice(None)) + (None,)*(c.ndim-2)]
1663 # Now fix continuity: on the very first interval, take the integration
1664 # constant to be zero; on an interval [x_j, x_{j+1}) with j>0,
1665 # the integration constant is then equal to the jump of the `bp` at x_j.
1666 # The latter is given by the coefficient of B_{n+1, n+1}
1667 # *on the previous interval* (other B. polynomials are zero at the
1668 # breakpoint). Finally, use the fact that BPs form a partition of unity.
1669 c2[:,1:] += np.cumsum(c2[k, :], axis=0)[:-1]
1671 if self.extrapolate == 'periodic':
1672 extrapolate = False
1673 else:
1674 extrapolate = self.extrapolate
1676 return self.construct_fast(c2, x, extrapolate, axis=self.axis)
1678 def integrate(self, a, b, extrapolate=None):
1679 """
1680 Compute a definite integral over a piecewise polynomial.
1682 Parameters
1683 ----------
1684 a : float
1685 Lower integration bound
1686 b : float
1687 Upper integration bound
1688 extrapolate : {bool, 'periodic', None}, optional
1689 Whether to extrapolate to out-of-bounds points based on first
1690 and last intervals, or to return NaNs. If 'periodic', periodic
1691 extrapolation is used. If None (default), use `self.extrapolate`.
1693 Returns
1694 -------
1695 array_like
1696 Definite integral of the piecewise polynomial over [a, b]
1698 """
1699 # XXX: can probably use instead the fact that
1700 # \int_0^{1} B_{j, n}(x) \dx = 1/(n+1)
1701 ib = self.antiderivative()
1702 if extrapolate is None:
1703 extrapolate = self.extrapolate
1705 # ib.extrapolate shouldn't be 'periodic', it is converted to
1706 # False for 'periodic. in antiderivative() call.
1707 if extrapolate != 'periodic':
1708 ib.extrapolate = extrapolate
1710 if extrapolate == 'periodic':
1711 # Split the integral into the part over period (can be several
1712 # of them) and the remaining part.
1714 # For simplicity and clarity convert to a <= b case.
1715 if a <= b:
1716 sign = 1
1717 else:
1718 a, b = b, a
1719 sign = -1
1721 xs, xe = self.x[0], self.x[-1]
1722 period = xe - xs
1723 interval = b - a
1724 n_periods, left = divmod(interval, period)
1725 res = n_periods * (ib(xe) - ib(xs))
1727 # Map a and b to [xs, xe].
1728 a = xs + (a - xs) % period
1729 b = a + left
1731 # If b <= xe then we need to integrate over [a, b], otherwise
1732 # over [a, xe] and from xs to what is remained.
1733 if b <= xe:
1734 res += ib(b) - ib(a)
1735 else:
1736 res += ib(xe) - ib(a) + ib(xs + left + a - xe) - ib(xs)
1738 return sign * res
1739 else:
1740 return ib(b) - ib(a)
1742 def extend(self, c, x):
1743 k = max(self.c.shape[0], c.shape[0])
1744 self.c = self._raise_degree(self.c, k - self.c.shape[0])
1745 c = self._raise_degree(c, k - c.shape[0])
1746 return _PPolyBase.extend(self, c, x)
1747 extend.__doc__ = _PPolyBase.extend.__doc__
1749 @classmethod
1750 def from_power_basis(cls, pp, extrapolate=None):
1751 """
1752 Construct a piecewise polynomial in Bernstein basis
1753 from a power basis polynomial.
1755 Parameters
1756 ----------
1757 pp : PPoly
1758 A piecewise polynomial in the power basis
1759 extrapolate : bool or 'periodic', optional
1760 If bool, determines whether to extrapolate to out-of-bounds points
1761 based on first and last intervals, or to return NaNs.
1762 If 'periodic', periodic extrapolation is used. Default is True.
1763 """
1764 if not isinstance(pp, PPoly):
1765 raise TypeError(".from_power_basis only accepts PPoly instances. "
1766 "Got %s instead." % type(pp))
1768 dx = np.diff(pp.x)
1769 k = pp.c.shape[0] - 1 # polynomial order
1771 rest = (None,)*(pp.c.ndim-2)
1773 c = np.zeros_like(pp.c)
1774 for a in range(k+1):
1775 factor = pp.c[a] / comb(k, k-a) * dx[(slice(None),)+rest]**(k-a)
1776 for j in range(k-a, k+1):
1777 c[j] += factor * comb(j, k-a)
1779 if extrapolate is None:
1780 extrapolate = pp.extrapolate
1782 return cls.construct_fast(c, pp.x, extrapolate, pp.axis)
1784 @classmethod
1785 def from_derivatives(cls, xi, yi, orders=None, extrapolate=None):
1786 """Construct a piecewise polynomial in the Bernstein basis,
1787 compatible with the specified values and derivatives at breakpoints.
1789 Parameters
1790 ----------
1791 xi : array_like
1792 sorted 1-D array of x-coordinates
1793 yi : array_like or list of array_likes
1794 ``yi[i][j]`` is the ``j``th derivative known at ``xi[i]``
1795 orders : None or int or array_like of ints. Default: None.
1796 Specifies the degree of local polynomials. If not None, some
1797 derivatives are ignored.
1798 extrapolate : bool or 'periodic', optional
1799 If bool, determines whether to extrapolate to out-of-bounds points
1800 based on first and last intervals, or to return NaNs.
1801 If 'periodic', periodic extrapolation is used. Default is True.
1803 Notes
1804 -----
1805 If ``k`` derivatives are specified at a breakpoint ``x``, the
1806 constructed polynomial is exactly ``k`` times continuously
1807 differentiable at ``x``, unless the ``order`` is provided explicitly.
1808 In the latter case, the smoothness of the polynomial at
1809 the breakpoint is controlled by the ``order``.
1811 Deduces the number of derivatives to match at each end
1812 from ``order`` and the number of derivatives available. If
1813 possible it uses the same number of derivatives from
1814 each end; if the number is odd it tries to take the
1815 extra one from y2. In any case if not enough derivatives
1816 are available at one end or another it draws enough to
1817 make up the total from the other end.
1819 If the order is too high and not enough derivatives are available,
1820 an exception is raised.
1822 Examples
1823 --------
1825 >>> from scipy.interpolate import BPoly
1826 >>> BPoly.from_derivatives([0, 1], [[1, 2], [3, 4]])
1828 Creates a polynomial `f(x)` of degree 3, defined on `[0, 1]`
1829 such that `f(0) = 1, df/dx(0) = 2, f(1) = 3, df/dx(1) = 4`
1831 >>> BPoly.from_derivatives([0, 1, 2], [[0, 1], [0], [2]])
1833 Creates a piecewise polynomial `f(x)`, such that
1834 `f(0) = f(1) = 0`, `f(2) = 2`, and `df/dx(0) = 1`.
1835 Based on the number of derivatives provided, the order of the
1836 local polynomials is 2 on `[0, 1]` and 1 on `[1, 2]`.
1837 Notice that no restriction is imposed on the derivatives at
1838 ``x = 1`` and ``x = 2``.
1840 Indeed, the explicit form of the polynomial is::
1842 f(x) = | x * (1 - x), 0 <= x < 1
1843 | 2 * (x - 1), 1 <= x <= 2
1845 So that f'(1-0) = -1 and f'(1+0) = 2
1847 """
1848 xi = np.asarray(xi)
1849 if len(xi) != len(yi):
1850 raise ValueError("xi and yi need to have the same length")
1851 if np.any(xi[1:] - xi[:1] <= 0):
1852 raise ValueError("x coordinates are not in increasing order")
1854 # number of intervals
1855 m = len(xi) - 1
1857 # global poly order is k-1, local orders are <=k and can vary
1858 try:
1859 k = max(len(yi[i]) + len(yi[i+1]) for i in range(m))
1860 except TypeError as e:
1861 raise ValueError(
1862 "Using a 1-D array for y? Please .reshape(-1, 1)."
1863 ) from e
1865 if orders is None:
1866 orders = [None] * m
1867 else:
1868 if isinstance(orders, (int, np.integer)):
1869 orders = [orders] * m
1870 k = max(k, max(orders))
1872 if any(o <= 0 for o in orders):
1873 raise ValueError("Orders must be positive.")
1875 c = []
1876 for i in range(m):
1877 y1, y2 = yi[i], yi[i+1]
1878 if orders[i] is None:
1879 n1, n2 = len(y1), len(y2)
1880 else:
1881 n = orders[i]+1
1882 n1 = min(n//2, len(y1))
1883 n2 = min(n - n1, len(y2))
1884 n1 = min(n - n2, len(y2))
1885 if n1+n2 != n:
1886 mesg = ("Point %g has %d derivatives, point %g"
1887 " has %d derivatives, but order %d requested" % (
1888 xi[i], len(y1), xi[i+1], len(y2), orders[i]))
1889 raise ValueError(mesg)
1891 if not (n1 <= len(y1) and n2 <= len(y2)):
1892 raise ValueError("`order` input incompatible with"
1893 " length y1 or y2.")
1895 b = BPoly._construct_from_derivatives(xi[i], xi[i+1],
1896 y1[:n1], y2[:n2])
1897 if len(b) < k:
1898 b = BPoly._raise_degree(b, k - len(b))
1899 c.append(b)
1901 c = np.asarray(c)
1902 return cls(c.swapaxes(0, 1), xi, extrapolate)
1904 @staticmethod
1905 def _construct_from_derivatives(xa, xb, ya, yb):
1906 r"""Compute the coefficients of a polynomial in the Bernstein basis
1907 given the values and derivatives at the edges.
1909 Return the coefficients of a polynomial in the Bernstein basis
1910 defined on ``[xa, xb]`` and having the values and derivatives at the
1911 endpoints `xa` and `xb` as specified by `ya`` and `yb`.
1912 The polynomial constructed is of the minimal possible degree, i.e.,
1913 if the lengths of `ya` and `yb` are `na` and `nb`, the degree
1914 of the polynomial is ``na + nb - 1``.
1916 Parameters
1917 ----------
1918 xa : float
1919 Left-hand end point of the interval
1920 xb : float
1921 Right-hand end point of the interval
1922 ya : array_like
1923 Derivatives at `xa`. `ya[0]` is the value of the function, and
1924 `ya[i]` for ``i > 0`` is the value of the ``i``th derivative.
1925 yb : array_like
1926 Derivatives at `xb`.
1928 Returns
1929 -------
1930 array
1931 coefficient array of a polynomial having specified derivatives
1933 Notes
1934 -----
1935 This uses several facts from life of Bernstein basis functions.
1936 First of all,
1938 .. math:: b'_{a, n} = n (b_{a-1, n-1} - b_{a, n-1})
1940 If B(x) is a linear combination of the form
1942 .. math:: B(x) = \sum_{a=0}^{n} c_a b_{a, n},
1944 then :math: B'(x) = n \sum_{a=0}^{n-1} (c_{a+1} - c_{a}) b_{a, n-1}.
1945 Iterating the latter one, one finds for the q-th derivative
1947 .. math:: B^{q}(x) = n!/(n-q)! \sum_{a=0}^{n-q} Q_a b_{a, n-q},
1949 with
1951 .. math:: Q_a = \sum_{j=0}^{q} (-)^{j+q} comb(q, j) c_{j+a}
1953 This way, only `a=0` contributes to :math: `B^{q}(x = xa)`, and
1954 `c_q` are found one by one by iterating `q = 0, ..., na`.
1956 At ``x = xb`` it's the same with ``a = n - q``.
1958 """
1959 ya, yb = np.asarray(ya), np.asarray(yb)
1960 if ya.shape[1:] != yb.shape[1:]:
1961 raise ValueError('Shapes of ya {} and yb {} are incompatible'
1962 .format(ya.shape, yb.shape))
1964 dta, dtb = ya.dtype, yb.dtype
1965 if (np.issubdtype(dta, np.complexfloating) or
1966 np.issubdtype(dtb, np.complexfloating)):
1967 dt = np.complex_
1968 else:
1969 dt = np.float_
1971 na, nb = len(ya), len(yb)
1972 n = na + nb
1974 c = np.empty((na+nb,) + ya.shape[1:], dtype=dt)
1976 # compute coefficients of a polynomial degree na+nb-1
1977 # walk left-to-right
1978 for q in range(0, na):
1979 c[q] = ya[q] / spec.poch(n - q, q) * (xb - xa)**q
1980 for j in range(0, q):
1981 c[q] -= (-1)**(j+q) * comb(q, j) * c[j]
1983 # now walk right-to-left
1984 for q in range(0, nb):
1985 c[-q-1] = yb[q] / spec.poch(n - q, q) * (-1)**q * (xb - xa)**q
1986 for j in range(0, q):
1987 c[-q-1] -= (-1)**(j+1) * comb(q, j+1) * c[-q+j]
1989 return c
1991 @staticmethod
1992 def _raise_degree(c, d):
1993 r"""Raise a degree of a polynomial in the Bernstein basis.
1995 Given the coefficients of a polynomial degree `k`, return (the
1996 coefficients of) the equivalent polynomial of degree `k+d`.
1998 Parameters
1999 ----------
2000 c : array_like
2001 coefficient array, 1-D
2002 d : integer
2004 Returns
2005 -------
2006 array
2007 coefficient array, 1-D array of length `c.shape[0] + d`
2009 Notes
2010 -----
2011 This uses the fact that a Bernstein polynomial `b_{a, k}` can be
2012 identically represented as a linear combination of polynomials of
2013 a higher degree `k+d`:
2015 .. math:: b_{a, k} = comb(k, a) \sum_{j=0}^{d} b_{a+j, k+d} \
2016 comb(d, j) / comb(k+d, a+j)
2018 """
2019 if d == 0:
2020 return c
2022 k = c.shape[0] - 1
2023 out = np.zeros((c.shape[0] + d,) + c.shape[1:], dtype=c.dtype)
2025 for a in range(c.shape[0]):
2026 f = c[a] * comb(k, a)
2027 for j in range(d+1):
2028 out[a+j] += f * comb(d, j) / comb(k+d, a+j)
2029 return out
2032class NdPPoly:
2033 """
2034 Piecewise tensor product polynomial
2036 The value at point ``xp = (x', y', z', ...)`` is evaluated by first
2037 computing the interval indices `i` such that::
2039 x[0][i[0]] <= x' < x[0][i[0]+1]
2040 x[1][i[1]] <= y' < x[1][i[1]+1]
2041 ...
2043 and then computing::
2045 S = sum(c[k0-m0-1,...,kn-mn-1,i[0],...,i[n]]
2046 * (xp[0] - x[0][i[0]])**m0
2047 * ...
2048 * (xp[n] - x[n][i[n]])**mn
2049 for m0 in range(k[0]+1)
2050 ...
2051 for mn in range(k[n]+1))
2053 where ``k[j]`` is the degree of the polynomial in dimension j. This
2054 representation is the piecewise multivariate power basis.
2056 Parameters
2057 ----------
2058 c : ndarray, shape (k0, ..., kn, m0, ..., mn, ...)
2059 Polynomial coefficients, with polynomial order `kj` and
2060 `mj+1` intervals for each dimension `j`.
2061 x : ndim-tuple of ndarrays, shapes (mj+1,)
2062 Polynomial breakpoints for each dimension. These must be
2063 sorted in increasing order.
2064 extrapolate : bool, optional
2065 Whether to extrapolate to out-of-bounds points based on first
2066 and last intervals, or to return NaNs. Default: True.
2068 Attributes
2069 ----------
2070 x : tuple of ndarrays
2071 Breakpoints.
2072 c : ndarray
2073 Coefficients of the polynomials.
2075 Methods
2076 -------
2077 __call__
2078 derivative
2079 antiderivative
2080 integrate
2081 integrate_1d
2082 construct_fast
2084 See also
2085 --------
2086 PPoly : piecewise polynomials in 1D
2088 Notes
2089 -----
2090 High-order polynomials in the power basis can be numerically
2091 unstable.
2093 """
2095 def __init__(self, c, x, extrapolate=None):
2096 self.x = tuple(np.ascontiguousarray(v, dtype=np.float64) for v in x)
2097 self.c = np.asarray(c)
2098 if extrapolate is None:
2099 extrapolate = True
2100 self.extrapolate = bool(extrapolate)
2102 ndim = len(self.x)
2103 if any(v.ndim != 1 for v in self.x):
2104 raise ValueError("x arrays must all be 1-dimensional")
2105 if any(v.size < 2 for v in self.x):
2106 raise ValueError("x arrays must all contain at least 2 points")
2107 if c.ndim < 2*ndim:
2108 raise ValueError("c must have at least 2*len(x) dimensions")
2109 if any(np.any(v[1:] - v[:-1] < 0) for v in self.x):
2110 raise ValueError("x-coordinates are not in increasing order")
2111 if any(a != b.size - 1 for a, b in zip(c.shape[ndim:2*ndim], self.x)):
2112 raise ValueError("x and c do not agree on the number of intervals")
2114 dtype = self._get_dtype(self.c.dtype)
2115 self.c = np.ascontiguousarray(self.c, dtype=dtype)
2117 @classmethod
2118 def construct_fast(cls, c, x, extrapolate=None):
2119 """
2120 Construct the piecewise polynomial without making checks.
2122 Takes the same parameters as the constructor. Input arguments
2123 ``c`` and ``x`` must be arrays of the correct shape and type. The
2124 ``c`` array can only be of dtypes float and complex, and ``x``
2125 array must have dtype float.
2127 """
2128 self = object.__new__(cls)
2129 self.c = c
2130 self.x = x
2131 if extrapolate is None:
2132 extrapolate = True
2133 self.extrapolate = extrapolate
2134 return self
2136 def _get_dtype(self, dtype):
2137 if np.issubdtype(dtype, np.complexfloating) \
2138 or np.issubdtype(self.c.dtype, np.complexfloating):
2139 return np.complex_
2140 else:
2141 return np.float_
2143 def _ensure_c_contiguous(self):
2144 if not self.c.flags.c_contiguous:
2145 self.c = self.c.copy()
2146 if not isinstance(self.x, tuple):
2147 self.x = tuple(self.x)
2149 def __call__(self, x, nu=None, extrapolate=None):
2150 """
2151 Evaluate the piecewise polynomial or its derivative
2153 Parameters
2154 ----------
2155 x : array-like
2156 Points to evaluate the interpolant at.
2157 nu : tuple, optional
2158 Orders of derivatives to evaluate. Each must be non-negative.
2159 extrapolate : bool, optional
2160 Whether to extrapolate to out-of-bounds points based on first
2161 and last intervals, or to return NaNs.
2163 Returns
2164 -------
2165 y : array-like
2166 Interpolated values. Shape is determined by replacing
2167 the interpolation axis in the original array with the shape of x.
2169 Notes
2170 -----
2171 Derivatives are evaluated piecewise for each polynomial
2172 segment, even if the polynomial is not differentiable at the
2173 breakpoints. The polynomial intervals are considered half-open,
2174 ``[a, b)``, except for the last interval which is closed
2175 ``[a, b]``.
2177 """
2178 if extrapolate is None:
2179 extrapolate = self.extrapolate
2180 else:
2181 extrapolate = bool(extrapolate)
2183 ndim = len(self.x)
2185 x = _ndim_coords_from_arrays(x)
2186 x_shape = x.shape
2187 x = np.ascontiguousarray(x.reshape(-1, x.shape[-1]), dtype=np.float_)
2189 if nu is None:
2190 nu = np.zeros((ndim,), dtype=np.intc)
2191 else:
2192 nu = np.asarray(nu, dtype=np.intc)
2193 if nu.ndim != 1 or nu.shape[0] != ndim:
2194 raise ValueError("invalid number of derivative orders nu")
2196 dim1 = prod(self.c.shape[:ndim])
2197 dim2 = prod(self.c.shape[ndim:2*ndim])
2198 dim3 = prod(self.c.shape[2*ndim:])
2199 ks = np.array(self.c.shape[:ndim], dtype=np.intc)
2201 out = np.empty((x.shape[0], dim3), dtype=self.c.dtype)
2202 self._ensure_c_contiguous()
2204 _ppoly.evaluate_nd(self.c.reshape(dim1, dim2, dim3),
2205 self.x,
2206 ks,
2207 x,
2208 nu,
2209 bool(extrapolate),
2210 out)
2212 return out.reshape(x_shape[:-1] + self.c.shape[2*ndim:])
2214 def _derivative_inplace(self, nu, axis):
2215 """
2216 Compute 1-D derivative along a selected dimension in-place
2217 May result to non-contiguous c array.
2218 """
2219 if nu < 0:
2220 return self._antiderivative_inplace(-nu, axis)
2222 ndim = len(self.x)
2223 axis = axis % ndim
2225 # reduce order
2226 if nu == 0:
2227 # noop
2228 return
2229 else:
2230 sl = [slice(None)]*ndim
2231 sl[axis] = slice(None, -nu, None)
2232 c2 = self.c[tuple(sl)]
2234 if c2.shape[axis] == 0:
2235 # derivative of order 0 is zero
2236 shp = list(c2.shape)
2237 shp[axis] = 1
2238 c2 = np.zeros(shp, dtype=c2.dtype)
2240 # multiply by the correct rising factorials
2241 factor = spec.poch(np.arange(c2.shape[axis], 0, -1), nu)
2242 sl = [None]*c2.ndim
2243 sl[axis] = slice(None)
2244 c2 *= factor[tuple(sl)]
2246 self.c = c2
2248 def _antiderivative_inplace(self, nu, axis):
2249 """
2250 Compute 1-D antiderivative along a selected dimension
2251 May result to non-contiguous c array.
2252 """
2253 if nu <= 0:
2254 return self._derivative_inplace(-nu, axis)
2256 ndim = len(self.x)
2257 axis = axis % ndim
2259 perm = list(range(ndim))
2260 perm[0], perm[axis] = perm[axis], perm[0]
2261 perm = perm + list(range(ndim, self.c.ndim))
2263 c = self.c.transpose(perm)
2265 c2 = np.zeros((c.shape[0] + nu,) + c.shape[1:],
2266 dtype=c.dtype)
2267 c2[:-nu] = c
2269 # divide by the correct rising factorials
2270 factor = spec.poch(np.arange(c.shape[0], 0, -1), nu)
2271 c2[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]
2273 # fix continuity of added degrees of freedom
2274 perm2 = list(range(c2.ndim))
2275 perm2[1], perm2[ndim+axis] = perm2[ndim+axis], perm2[1]
2277 c2 = c2.transpose(perm2)
2278 c2 = c2.copy()
2279 _ppoly.fix_continuity(c2.reshape(c2.shape[0], c2.shape[1], -1),
2280 self.x[axis], nu-1)
2282 c2 = c2.transpose(perm2)
2283 c2 = c2.transpose(perm)
2285 # Done
2286 self.c = c2
2288 def derivative(self, nu):
2289 """
2290 Construct a new piecewise polynomial representing the derivative.
2292 Parameters
2293 ----------
2294 nu : ndim-tuple of int
2295 Order of derivatives to evaluate for each dimension.
2296 If negative, the antiderivative is returned.
2298 Returns
2299 -------
2300 pp : NdPPoly
2301 Piecewise polynomial of orders (k[0] - nu[0], ..., k[n] - nu[n])
2302 representing the derivative of this polynomial.
2304 Notes
2305 -----
2306 Derivatives are evaluated piecewise for each polynomial
2307 segment, even if the polynomial is not differentiable at the
2308 breakpoints. The polynomial intervals in each dimension are
2309 considered half-open, ``[a, b)``, except for the last interval
2310 which is closed ``[a, b]``.
2312 """
2313 p = self.construct_fast(self.c.copy(), self.x, self.extrapolate)
2315 for axis, n in enumerate(nu):
2316 p._derivative_inplace(n, axis)
2318 p._ensure_c_contiguous()
2319 return p
2321 def antiderivative(self, nu):
2322 """
2323 Construct a new piecewise polynomial representing the antiderivative.
2325 Antiderivative is also the indefinite integral of the function,
2326 and derivative is its inverse operation.
2328 Parameters
2329 ----------
2330 nu : ndim-tuple of int
2331 Order of derivatives to evaluate for each dimension.
2332 If negative, the derivative is returned.
2334 Returns
2335 -------
2336 pp : PPoly
2337 Piecewise polynomial of order k2 = k + n representing
2338 the antiderivative of this polynomial.
2340 Notes
2341 -----
2342 The antiderivative returned by this function is continuous and
2343 continuously differentiable to order n-1, up to floating point
2344 rounding error.
2346 """
2347 p = self.construct_fast(self.c.copy(), self.x, self.extrapolate)
2349 for axis, n in enumerate(nu):
2350 p._antiderivative_inplace(n, axis)
2352 p._ensure_c_contiguous()
2353 return p
2355 def integrate_1d(self, a, b, axis, extrapolate=None):
2356 r"""
2357 Compute NdPPoly representation for one dimensional definite integral
2359 The result is a piecewise polynomial representing the integral:
2361 .. math::
2363 p(y, z, ...) = \int_a^b dx\, p(x, y, z, ...)
2365 where the dimension integrated over is specified with the
2366 `axis` parameter.
2368 Parameters
2369 ----------
2370 a, b : float
2371 Lower and upper bound for integration.
2372 axis : int
2373 Dimension over which to compute the 1-D integrals
2374 extrapolate : bool, optional
2375 Whether to extrapolate to out-of-bounds points based on first
2376 and last intervals, or to return NaNs.
2378 Returns
2379 -------
2380 ig : NdPPoly or array-like
2381 Definite integral of the piecewise polynomial over [a, b].
2382 If the polynomial was 1D, an array is returned,
2383 otherwise, an NdPPoly object.
2385 """
2386 if extrapolate is None:
2387 extrapolate = self.extrapolate
2388 else:
2389 extrapolate = bool(extrapolate)
2391 ndim = len(self.x)
2392 axis = int(axis) % ndim
2394 # reuse 1-D integration routines
2395 c = self.c
2396 swap = list(range(c.ndim))
2397 swap.insert(0, swap[axis])
2398 del swap[axis + 1]
2399 swap.insert(1, swap[ndim + axis])
2400 del swap[ndim + axis + 1]
2402 c = c.transpose(swap)
2403 p = PPoly.construct_fast(c.reshape(c.shape[0], c.shape[1], -1),
2404 self.x[axis],
2405 extrapolate=extrapolate)
2406 out = p.integrate(a, b, extrapolate=extrapolate)
2408 # Construct result
2409 if ndim == 1:
2410 return out.reshape(c.shape[2:])
2411 else:
2412 c = out.reshape(c.shape[2:])
2413 x = self.x[:axis] + self.x[axis+1:]
2414 return self.construct_fast(c, x, extrapolate=extrapolate)
2416 def integrate(self, ranges, extrapolate=None):
2417 """
2418 Compute a definite integral over a piecewise polynomial.
2420 Parameters
2421 ----------
2422 ranges : ndim-tuple of 2-tuples float
2423 Sequence of lower and upper bounds for each dimension,
2424 ``[(a[0], b[0]), ..., (a[ndim-1], b[ndim-1])]``
2425 extrapolate : bool, optional
2426 Whether to extrapolate to out-of-bounds points based on first
2427 and last intervals, or to return NaNs.
2429 Returns
2430 -------
2431 ig : array_like
2432 Definite integral of the piecewise polynomial over
2433 [a[0], b[0]] x ... x [a[ndim-1], b[ndim-1]]
2435 """
2437 ndim = len(self.x)
2439 if extrapolate is None:
2440 extrapolate = self.extrapolate
2441 else:
2442 extrapolate = bool(extrapolate)
2444 if not hasattr(ranges, '__len__') or len(ranges) != ndim:
2445 raise ValueError("Range not a sequence of correct length")
2447 self._ensure_c_contiguous()
2449 # Reuse 1D integration routine
2450 c = self.c
2451 for n, (a, b) in enumerate(ranges):
2452 swap = list(range(c.ndim))
2453 swap.insert(1, swap[ndim - n])
2454 del swap[ndim - n + 1]
2456 c = c.transpose(swap)
2458 p = PPoly.construct_fast(c, self.x[n], extrapolate=extrapolate)
2459 out = p.integrate(a, b, extrapolate=extrapolate)
2460 c = out.reshape(c.shape[2:])
2462 return c