scijit.interpolate.PchipInterpolator

scijit.interpolate.PchipInterpolator(x, y, axis=0, extrapolate=None)

Monotone piecewise-cubic Hermite interpolator (Fritsch-Carlson slopes).

Preserves monotonicity of the data and does not overshoot, at the price of being only C1 (the second derivative jumps at the nodes) where CubicSpline is C2.

Parameters:

x1-D float64 ndarray, length n >= 2

Breakpoints, strictly increasing. A non-increasing pair raises ValueError.

yfloat64 ndarray of rank 1, 2 or 3

Values at the breakpoints. Its length along axis must be n. Every position on the other axes is an independent series over the same x. Rank 4 and above raise ValueError.

axisint, optional

Which axis of y the interpolation runs along, default 0. Negative values count from the end. Ignored for a 1-D y. ev(xs) returns the shape of y with this axis replaced by len(xs); ev_one returns it with this axis removed.

extrapolateNone, bool or str, optional

What a query outside [x[0], x[-1]] gives:

True           the end segment polynomial, continued
False          NaN
'periodic'     the query folded into ``[x[0], x[-1]]`` first
None           the default, which is True

A string other than 'periodic' raises ValueError. 0 and 1 are accepted for False and True. May be a runtime value, not only a literal.

Attributes:

x1-D float64 ndarray

The breakpoints, as float64.

c(4, n-1) float64 ndarray, or (4, n-1, T) for an N-D y

PPoly coefficients: c[k, i] multiplies (x - x[i]) ** (3 - k) on [x[i], x[i+1]]. For an N-D y the trailing axes are flattened into one; get_coeffs unflattens them.

nint

Number of breakpoints.

extrapolateint

The resolved extrapolation mode: 0 NaN outside the data, 1 the end segment continued, 2 the query folded into the base period.

Methods

ev(xs), ev_one(x), __getitem__(xs), derivative_ev(xs, nu),
get_knots(), get_coeffs()

Returns:

_PchipInterpolator

A callable jitclass instance: f(xs) runs ev and f(x) runs ev_one. An N-D y gives _PchipInterpolatorND1 or _PchipInterpolatorND2, which carry the same methods.

See also

scipy.interpolate.PchipInterpolator

The scipy routine this mirrors.

Notes

Deviations from scipy. extrapolate reads back as the resolved int code rather than the caller’s value, and an unrecognised string raises where scipy extrapolates. y is capped at rank 3 and is float64, where scipy takes any rank; a jitclass field has a fixed rank and a fixed dtype, so each supported rank is a separate class. A rank-2 x raises TypingError where scipy raises ValueError('`x` must be 1-dimensional.'); both refuse and only the class differs. A complex y raises ValueError naming np.real, which is scipy’s own refusal, from Python and from inside @njit. Evaluation is ev(xs) and derivative_ev(xs, nu), where scipy spells both as f(xs, nu); a two-argument call raises on arity here. There is no .derivative() or .antiderivative() object-returning API, no .roots(), and no integral method, which CubicSpline here does carry.

PchipInterpolator is an @njit factory over the jitclass _PchipInterpolator, which is the shape every other evaluable name in this subpackage uses, so extrapolate may be omitted from Python AND from inside @njit: a jitclass constructor’s defaults are Python-only and a plain @njit function’s are not.

Measured against scipy 1.18 on 5 nodes spanning [0, 1.4], queried at 2.0 and -1.0: extrapolate=None gives 7.07678571 and 10.1547619 on both sides, extrapolate=False NaN on both, extrapolate='periodic' 0.7145361 and 0.8756787.

Accuracy vs scipy.interpolate.PchipInterpolator on 15 random nodes over 211 points: values 4.4e-16, first derivative 7.1e-15, and the PPoly coefficient array exactly 0.0.

prange-safe: yes – pure numpy, no state.

Examples

>>> import numpy as np
>>> from numba import njit
>>> from scijit.interpolate import PchipInterpolator
>>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
>>> y = np.array([0.0, 1.0, 1.0, 3.0, 3.0])
>>> f = PchipInterpolator(x, y)
>>> np.round(f(np.array([0.5, 2.5])), 6)
array([0.6875, 2.    ])

Inside compiled code, scanning the domain for the interpolant’s peak:

>>> @njit
... def peak(f):
...     xs = np.linspace(0.0, 4.0, 201)
...     return np.max(f(xs))
>>> float(np.round(peak(f), 6))
3.0