scirs2_interpolate/interp1d/
pchip.rs

1//! PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) interpolation
2//!
3//! This module provides an implementation of the PCHIP algorithm, which produces
4//! shape-preserving interpolants that maintain monotonicity in the data.
5
6use crate::error::{InterpolateError, InterpolateResult};
7use scirs2_core::ndarray::{Array1, ArrayView1};
8use scirs2_core::numeric::{Float, FromPrimitive};
9use std::fmt::Debug;
10
11/// Strategy for PCHIP extrapolation outside the data range.
12#[derive(Debug, Clone, Copy, PartialEq, Default)]
13pub enum PchipExtrapolateMode {
14    /// Linear extension using the endpoint derivative (bounded growth, stable
15    /// for far extrapolation). This is the default.
16    #[default]
17    Linear,
18    /// Hermite polynomial continuation of the boundary segment. Matches
19    /// `scipy.interpolate.PchipInterpolator(extrapolate=True)` but can grow
20    /// cubically for points far from the data range.
21    Polynomial,
22}
23
24/// PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) interpolator
25///
26/// This interpolator preserves monotonicity in the interpolation data and does
27/// not overshoot if the data is not smooth. The first derivatives are guaranteed
28/// to be continuous, but the second derivatives may jump at the knot points.
29///
30/// The algorithm determines the derivatives at points x_k by using the PCHIP algorithm
31/// from Fritsch and Butland (1984).
32///
33/// # References
34///
35/// 1. F. N. Fritsch and J. Butland, "A method for constructing local monotone
36///    piecewise cubic interpolants", SIAM J. Sci. Comput., 5(2), 300-304 (1984).
37/// 2. C. Moler, "Numerical Computing with Matlab", 2004.
38#[derive(Debug, Clone)]
39pub struct PchipInterpolator<F: Float> {
40    /// X coordinates (must be sorted)
41    x: Array1<F>,
42    /// Y coordinates
43    y: Array1<F>,
44    /// Derivatives at points
45    derivatives: Array1<F>,
46    /// Extrapolation mode
47    extrapolate: bool,
48    /// Extrapolation strategy (linear vs polynomial continuation)
49    extrapolate_mode: PchipExtrapolateMode,
50}
51
52impl<F: Float + FromPrimitive + Debug> PchipInterpolator<F> {
53    /// Create a new PCHIP interpolator
54    ///
55    /// # Arguments
56    ///
57    /// * `x` - The x coordinates (must be sorted in ascending order)
58    /// * `y` - The y coordinates (must have the same length as x)
59    /// * `extrapolate` - Whether to extrapolate beyond the data range
60    ///
61    /// # Returns
62    ///
63    /// A new `PchipInterpolator` object
64    ///
65    /// # Errors
66    ///
67    /// Returns an error if:
68    /// - `x` and `y` have different lengths
69    /// - `x` is not sorted in ascending order
70    /// - There are fewer than 2 points
71    /// - `y` contains complex values (not supported)
72    ///
73    /// # Examples
74    ///
75    /// ```
76    /// use scirs2_core::ndarray::array;
77    /// use scirs2_interpolate::interp1d::pchip::PchipInterpolator;
78    ///
79    /// let x = array![0.0f64, 1.0, 2.0, 3.0];
80    /// let y = array![0.0f64, 1.0, 4.0, 9.0];
81    ///
82    /// // Create a PCHIP interpolator
83    /// let interp = PchipInterpolator::new(&x.view(), &y.view(), true).expect("Operation failed");
84    ///
85    /// // Interpolate at x = 1.5
86    /// let y_interp = interp.evaluate(1.5);
87    /// ```
88    pub fn new(x: &ArrayView1<F>, y: &ArrayView1<F>, extrapolate: bool) -> InterpolateResult<Self> {
89        // Check inputs
90        if x.len() != y.len() {
91            return Err(InterpolateError::invalid_input(
92                "x and y arrays must have the same length".to_string(),
93            ));
94        }
95
96        if x.len() < 2 {
97            return Err(InterpolateError::insufficient_points(
98                2,
99                x.len(),
100                "PCHIP interpolation",
101            ));
102        }
103
104        // Check that x is sorted
105        for i in 1..x.len() {
106            if x[i] <= x[i - 1] {
107                return Err(InterpolateError::invalid_input(
108                    "x values must be sorted in ascending order".to_string(),
109                ));
110            }
111        }
112
113        // Clone input arrays
114        let x_arr = x.to_owned();
115        let y_arr = y.to_owned();
116
117        // Compute derivatives at points
118        let derivatives = Self::find_derivatives(&x_arr, &y_arr)?;
119
120        Ok(PchipInterpolator {
121            x: x_arr,
122            y: y_arr,
123            derivatives,
124            extrapolate,
125            extrapolate_mode: PchipExtrapolateMode::Linear,
126        })
127    }
128
129    /// Set the extrapolation mode (builder pattern).
130    ///
131    /// # Arguments
132    ///
133    /// * `mode` - The extrapolation strategy to use
134    pub fn with_extrapolate_mode(mut self, mode: PchipExtrapolateMode) -> Self {
135        self.extrapolate_mode = mode;
136        self
137    }
138
139    /// Evaluate the interpolation at the given point
140    ///
141    /// # Arguments
142    ///
143    /// * `xnew` - The x coordinate at which to evaluate the interpolation
144    ///
145    /// # Returns
146    ///
147    /// The interpolated y value at `xnew`
148    ///
149    /// # Errors
150    ///
151    /// Returns an error if `xnew` is outside the interpolation range and
152    /// extrapolation is disabled.
153    pub fn evaluate(&self, xnew: F) -> InterpolateResult<F> {
154        let n = self.x.len();
155
156        // Check if we're extrapolating
157        let is_extrapolating = xnew < self.x[0] || xnew > self.x[n - 1];
158        if is_extrapolating && !self.extrapolate {
159            return Err(InterpolateError::OutOfBounds(
160                "xnew is outside the interpolation range".to_string(),
161            ));
162        }
163
164        // Handle linear extrapolation mode
165        if is_extrapolating && self.extrapolate_mode == PchipExtrapolateMode::Linear {
166            if xnew < self.x[0] {
167                let dx = xnew - self.x[0];
168                return Ok(self.y[0] + self.derivatives[0] * dx);
169            } else {
170                let dx = xnew - self.x[n - 1];
171                return Ok(self.y[n - 1] + self.derivatives[n - 1] * dx);
172            }
173        }
174
175        // Special case: xnew is exactly the last point (non-extrapolating)
176        if !is_extrapolating && xnew == self.x[n - 1] {
177            return Ok(self.y[n - 1]);
178        }
179
180        // Find the segment index
181        let idx = if is_extrapolating {
182            // Polynomial extrapolation: use boundary segment
183            if xnew < self.x[0] {
184                0
185            } else {
186                n - 2
187            }
188        } else {
189            // In-bounds: binary/linear search
190            let mut seg = 0;
191            for i in 0..n - 1 {
192                if xnew >= self.x[i] && xnew <= self.x[i + 1] {
193                    seg = i;
194                    break;
195                }
196            }
197            seg
198        };
199
200        // Get coordinates and derivatives for the segment
201        let x1 = self.x[idx];
202        let x2 = self.x[idx + 1];
203        let y1 = self.y[idx];
204        let y2 = self.y[idx + 1];
205        let d1 = self.derivatives[idx];
206        let d2 = self.derivatives[idx + 1];
207
208        // Normalized position within the interval [x1, x2]
209        // For polynomial extrapolation t may be outside [0,1]
210        let h = x2 - x1;
211        let t = (xnew - x1) / h;
212
213        // Compute Hermite basis functions
214        let h00 = Self::h00(t);
215        let h10 = Self::h10(t);
216        let h01 = Self::h01(t);
217        let h11 = Self::h11(t);
218
219        // Evaluate cubic Hermite polynomial
220        let result = h00 * y1 + h10 * h * d1 + h01 * y2 + h11 * h * d2;
221
222        Ok(result)
223    }
224
225    /// Evaluate the interpolation at multiple points
226    ///
227    /// # Arguments
228    ///
229    /// * `xnew` - The x coordinates at which to evaluate the interpolation
230    ///
231    /// # Returns
232    ///
233    /// The interpolated y values at `xnew`
234    ///
235    /// # Errors
236    ///
237    /// Returns an error if any point in `xnew` is outside the interpolation range
238    /// and extrapolation is disabled.
239    pub fn evaluate_array(&self, xnew: &ArrayView1<F>) -> InterpolateResult<Array1<F>> {
240        let mut result = Array1::zeros(xnew.len());
241        for (i, &x) in xnew.iter().enumerate() {
242            result[i] = self.evaluate(x)?;
243        }
244        Ok(result)
245    }
246
247    /// Hermite basis function h₀₀(t)
248    fn h00(t: F) -> F {
249        let two = F::from_f64(2.0).expect("Operation failed");
250        let three = F::from_f64(3.0).expect("Operation failed");
251        (two * t * t * t) - (three * t * t) + F::one()
252    }
253
254    /// Hermite basis function h₁₀(t)
255    fn h10(t: F) -> F {
256        let two = F::from_f64(2.0).expect("Operation failed");
257        // three is not used in this function
258        (t * t * t) - (two * t * t) + t
259    }
260
261    /// Hermite basis function h₀₁(t)
262    fn h01(t: F) -> F {
263        let two = F::from_f64(2.0).expect("Operation failed");
264        let three = F::from_f64(3.0).expect("Operation failed");
265        -(two * t * t * t) + (three * t * t)
266    }
267
268    /// Hermite basis function h₁₁(t)
269    fn h11(t: F) -> F {
270        // No constants needed for this function
271        (t * t * t) - (t * t)
272    }
273
274    /// Helper for edge case derivative estimation
275    ///
276    /// # Arguments
277    ///
278    /// * `h0` - Spacing: x[1] - x[0] for the first point, or x[n-1] - x[n-2] for the last point
279    /// * `h1` - Spacing: x[2] - x[1] for the first point, or x[n] - x[n-1] for the last point
280    /// * `m0` - Slope of the first segment
281    /// * `m1` - Slope of the second segment
282    ///
283    /// # Returns
284    ///
285    /// The estimated derivative at the endpoint
286    fn edge_case(h0: F, h1: F, m0: F, m1: F) -> F {
287        // One-sided three-point estimate for the derivative
288        let two = F::from_f64(2.0).expect("Operation failed");
289        let three = F::from_f64(3.0).expect("Operation failed");
290
291        let d = ((two * h0 + h1) * m0 - h0 * m1) / (h0 + h1);
292
293        // Try to preserve shape
294        let sign_d = if d >= F::zero() { F::one() } else { -F::one() };
295        let sign_m0 = if m0 >= F::zero() { F::one() } else { -F::one() };
296        let sign_m1 = if m1 >= F::zero() { F::one() } else { -F::one() };
297
298        // If the signs are different or abs(d) > 3*abs(m0), adjust d
299        if sign_d != sign_m0 {
300            F::zero()
301        } else if (sign_m0 != sign_m1) && (d.abs() > three * m0.abs()) {
302            three * m0
303        } else {
304            d
305        }
306    }
307
308    /// Compute derivatives at points for PCHIP interpolation
309    ///
310    /// # Arguments
311    ///
312    /// * `x` - X coordinates
313    /// * `y` - Y coordinates
314    ///
315    /// # Returns
316    ///
317    /// Array of derivatives, one for each point in x
318    fn find_derivatives(x: &Array1<F>, y: &Array1<F>) -> InterpolateResult<Array1<F>> {
319        let n = x.len();
320        let mut derivatives = Array1::zeros(n);
321
322        // Handle special case: only two points, use linear interpolation
323        if n == 2 {
324            let slope = (y[1] - y[0]) / (x[1] - x[0]);
325            derivatives[0] = slope;
326            derivatives[1] = slope;
327            return Ok(derivatives);
328        }
329
330        // Calculate slopes between segments: m_k = (y_{k+1} - y_k) / (x_{k+1} - x_k)
331        let mut slopes = Array1::zeros(n - 1);
332        for i in 0..n - 1 {
333            slopes[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
334        }
335
336        // Calculate spacings: h_k = x_{k+1} - x_k
337        let mut h = Array1::zeros(n - 1);
338        for i in 0..n - 1 {
339            h[i] = x[i + 1] - x[i];
340        }
341
342        // For interior points, use PCHIP formula
343        let two = F::from_f64(2.0).expect("Operation failed");
344
345        for i in 1..n - 1 {
346            // Determine if slopes have different signs or if either is zero
347            let prev_slope = slopes[i - 1];
348            let curr_slope = slopes[i];
349
350            let sign_prev = if prev_slope > F::zero() {
351                F::one()
352            } else if prev_slope < F::zero() {
353                -F::one()
354            } else {
355                F::zero()
356            };
357
358            let sign_curr = if curr_slope > F::zero() {
359                F::one()
360            } else if curr_slope < F::zero() {
361                -F::one()
362            } else {
363                F::zero()
364            };
365
366            // If signs are different or either slope is zero, set derivative to zero
367            if sign_prev * sign_curr <= F::zero() {
368                derivatives[i] = F::zero();
369            } else {
370                // Use weighted harmonic mean
371                let w1 = two * h[i] + h[i - 1];
372                let w2 = h[i] + two * h[i - 1];
373
374                // Compute harmonic mean
375                if prev_slope.abs() < F::epsilon() || curr_slope.abs() < F::epsilon() {
376                    derivatives[i] = F::zero();
377                } else {
378                    let whmean_inv = (w1 / prev_slope + w2 / curr_slope) / (w1 + w2);
379                    derivatives[i] = F::one() / whmean_inv;
380                }
381            }
382        }
383
384        // Special case: endpoints
385        // For the first point
386        derivatives[0] = Self::edge_case(h[0], h[1], slopes[0], slopes[1]);
387
388        // For the last point
389        derivatives[n - 1] = Self::edge_case(h[n - 2], h[n - 3], slopes[n - 2], slopes[n - 3]);
390
391        Ok(derivatives)
392    }
393}
394
395/// PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) interpolation convenience function
396///
397/// # Arguments
398///
399/// * `x` - The x coordinates (must be sorted in ascending order)
400/// * `y` - The y coordinates
401/// * `xnew` - The points at which to interpolate
402/// * `extrapolate` - Whether to extrapolate beyond the data range
403///
404/// # Returns
405///
406/// The interpolated values
407///
408/// # Examples
409///
410/// ```
411/// use scirs2_core::ndarray::array;
412/// use scirs2_interpolate::pchip_interpolate;
413///
414/// let x = array![0.0f64, 1.0, 2.0, 3.0];
415/// let y = array![0.0f64, 1.0, 4.0, 9.0];
416/// let xnew = array![0.5f64, 1.5, 2.5];
417///
418/// let y_interp = pchip_interpolate(&x.view(), &y.view(), &xnew.view(), true).expect("Operation failed");
419/// ```
420#[allow(dead_code)]
421pub fn pchip_interpolate<F: crate::traits::InterpolationFloat>(
422    x: &ArrayView1<F>,
423    y: &ArrayView1<F>,
424    xnew: &ArrayView1<F>,
425    extrapolate: bool,
426) -> InterpolateResult<Array1<F>> {
427    let interp = PchipInterpolator::new(x, y, extrapolate)?;
428    interp.evaluate_array(xnew)
429}
430
431#[cfg(test)]
432mod tests {
433    use super::*;
434    use approx::assert_relative_eq;
435    use scirs2_core::ndarray::array;
436
437    #[test]
438    fn test_pchip_interpolation_basic() {
439        let x = array![0.0, 1.0, 2.0, 3.0];
440        let y = array![0.0, 1.0, 4.0, 9.0];
441
442        let interp = PchipInterpolator::new(&x.view(), &y.view(), false).expect("Operation failed");
443
444        // Test points exactly at data points
445        assert_relative_eq!(interp.evaluate(0.0).expect("Operation failed"), 0.0);
446        assert_relative_eq!(interp.evaluate(1.0).expect("Operation failed"), 1.0);
447        assert_relative_eq!(interp.evaluate(2.0).expect("Operation failed"), 4.0);
448        assert_relative_eq!(interp.evaluate(3.0).expect("Operation failed"), 9.0);
449
450        // Test points between data points - PCHIP should preserve shape
451        let y_interp_0_5 = interp.evaluate(0.5).expect("Operation failed");
452        let y_interp_1_5 = interp.evaluate(1.5).expect("Operation failed");
453        let y_interp_2_5 = interp.evaluate(2.5).expect("Operation failed");
454
455        // For this monotonically increasing data, PCHIP should preserve monotonicity
456        assert!(y_interp_0_5 > 0.0 && y_interp_0_5 < 1.0);
457        assert!(y_interp_1_5 > 1.0 && y_interp_1_5 < 4.0);
458        assert!(y_interp_2_5 > 4.0 && y_interp_2_5 < 9.0);
459    }
460
461    #[test]
462    fn test_pchip_monotonicity_preservation() {
463        // Test with data that has monotonic segments
464        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
465        let y = array![0.0, 1.0, 0.5, 0.0, 0.5, 2.0];
466
467        let interp = PchipInterpolator::new(&x.view(), &y.view(), false).expect("Operation failed");
468
469        // Check monotonicity preservation in the first segment (increasing)
470        let y_0_25 = interp.evaluate(0.25).expect("Operation failed");
471        let y_0_50 = interp.evaluate(0.50).expect("Operation failed");
472        let y_0_75 = interp.evaluate(0.75).expect("Operation failed");
473        assert!(y_0_25 <= y_0_50 && y_0_50 <= y_0_75);
474
475        // Check monotonicity preservation in the second segment (decreasing)
476        let y_1_25 = interp.evaluate(1.25).expect("Operation failed");
477        let y_1_50 = interp.evaluate(1.50).expect("Operation failed");
478        let y_1_75 = interp.evaluate(1.75).expect("Operation failed");
479        assert!(y_1_25 >= y_1_50 && y_1_50 >= y_1_75);
480
481        // Check monotonicity preservation in the last segment (increasing)
482        let y_4_25 = interp.evaluate(4.25).expect("Operation failed");
483        let y_4_50 = interp.evaluate(4.50).expect("Operation failed");
484        let y_4_75 = interp.evaluate(4.75).expect("Operation failed");
485        assert!(y_4_25 <= y_4_50 && y_4_50 <= y_4_75);
486    }
487
488    #[test]
489    fn test_pchip_extrapolation() {
490        let x = array![0.0, 1.0, 2.0, 3.0];
491        let y = array![0.0, 1.0, 4.0, 9.0];
492
493        // Test with extrapolation enabled
494        let interp_extrap =
495            PchipInterpolator::new(&x.view(), &y.view(), true).expect("Operation failed");
496        let y_minus_1 = interp_extrap.evaluate(-1.0).expect("Operation failed");
497        let y_plus_4 = interp_extrap.evaluate(4.0).expect("Operation failed");
498
499        // For this data, PCHIP computes derivatives that preserve shape
500        // The extrapolation uses linear extension based on endpoint derivatives
501        // Verify that extrapolation above the range increases beyond the last point
502        assert!(
503            y_plus_4 > 9.0,
504            "Extrapolation above should be greater than last point"
505        );
506
507        // Verify that extrapolation is finite (not NaN or infinite)
508        assert!(
509            y_minus_1.is_finite(),
510            "Extrapolation should produce finite values"
511        );
512        assert!(
513            y_plus_4.is_finite(),
514            "Extrapolation should produce finite values"
515        );
516
517        // Test with extrapolation disabled
518        let interp_no_extrap =
519            PchipInterpolator::new(&x.view(), &y.view(), false).expect("Operation failed");
520        assert!(interp_no_extrap.evaluate(-1.0).is_err());
521        assert!(interp_no_extrap.evaluate(4.0).is_err());
522    }
523
524    #[test]
525    fn test_pchip_extrapolation_far_beyond_range() {
526        // Regression test for issue #96
527        // PCHIP extrapolation should use linear extension, not cubic polynomial
528        let x = array![0.0, 1.0, 2.0, 3.0];
529        let y = array![0.0, 1.0, 4.0, 9.0];
530
531        let interp = PchipInterpolator::new(&x.view(), &y.view(), true).expect("Operation failed");
532
533        // Test far extrapolation (the bug case)
534        let y_50 = interp.evaluate(50.0).expect("Operation failed");
535
536        // THE KEY FIX: Before the fix, this returned -24008 (wrong!)
537        // After the fix with linear extrapolation, it should:
538        // 1. Be positive (not -24008 as in the bug)
539        // 2. Be greater than the last data point (9.0)
540        // 3. Be reasonable (linear extension, not explosive)
541        assert!(
542            y_50 > 9.0,
543            "Extrapolation at x=50 should be > 9.0, got {}",
544            y_50
545        );
546        assert!(
547            y_50 < 1000.0,
548            "Extrapolation should be reasonable (linear), got {}",
549            y_50
550        );
551        assert!(
552            y_50.is_finite(),
553            "Extrapolation should produce finite values"
554        );
555
556        // Test far extrapolation in the negative direction
557        // Should be finite and reasonable (not explosive)
558        let y_minus_50 = interp.evaluate(-50.0).expect("Operation failed");
559        assert!(
560            y_minus_50.is_finite(),
561            "Extrapolation should produce finite values"
562        );
563        assert!(
564            y_minus_50.abs() < 1000.0,
565            "Extrapolation should be reasonable (linear), got {}",
566            y_minus_50
567        );
568    }
569
570    #[test]
571    fn test_pchip_extrapolation_linear_behavior() {
572        // Verify that extrapolation uses linear extension (not cubic)
573        let x = array![0.0, 1.0, 2.0, 3.0];
574        let y = array![0.0, 1.0, 4.0, 9.0];
575
576        let interp = PchipInterpolator::new(&x.view(), &y.view(), true).expect("Operation failed");
577
578        // Get extrapolated values at two points beyond the range
579        let y_4 = interp.evaluate(4.0).expect("Operation failed");
580        let y_5 = interp.evaluate(5.0).expect("Operation failed");
581
582        // If using linear extrapolation, the slope should be constant
583        // slope = (y_5 - y_4) / (5.0 - 4.0) = y_5 - y_4
584        let extrap_slope = y_5 - y_4;
585
586        // The slope should be equal to the derivative at the last point
587        let last_derivative = interp.derivatives[interp.derivatives.len() - 1];
588
589        // Allow small numerical differences
590        assert_relative_eq!(extrap_slope, last_derivative, epsilon = 1e-10);
591    }
592
593    #[test]
594    fn test_pchip_interpolate_function() {
595        let x = array![0.0, 1.0, 2.0, 3.0];
596        let y = array![0.0, 1.0, 4.0, 9.0];
597        let xnew = array![0.5, 1.5, 2.5];
598
599        let y_interp =
600            pchip_interpolate(&x.view(), &y.view(), &xnew.view(), false).expect("Operation failed");
601
602        // Test that the function returns the expected number of points
603        assert_eq!(y_interp.len(), 3);
604
605        // For this monotonically increasing data, all interpolated values should be monotonically increasing
606        assert!(y_interp[0] > 0.0 && y_interp[0] < 1.0);
607        assert!(y_interp[1] > 1.0 && y_interp[1] < 4.0);
608        assert!(y_interp[2] > 4.0 && y_interp[2] < 9.0);
609    }
610
611    #[test]
612    fn test_pchip_error_conditions() {
613        // Test with different length arrays
614        let x = array![0.0, 1.0, 2.0, 3.0];
615        let y = array![0.0, 1.0, 4.0];
616        assert!(PchipInterpolator::new(&x.view(), &y.view(), false).is_err());
617
618        // Test with unsorted x
619        let x_unsorted = array![0.0, 2.0, 1.0, 3.0];
620        let y_valid = array![0.0, 1.0, 4.0, 9.0];
621        assert!(PchipInterpolator::new(&x_unsorted.view(), &y_valid.view(), false).is_err());
622
623        // Test with too few points
624        let x_short = array![0.0];
625        let y_short = array![0.0];
626        assert!(PchipInterpolator::new(&x_short.view(), &y_short.view(), false).is_err());
627    }
628}
scirs2_interpolate/interp1d/pchip.rs

scirs2_interpolate/interp1d/
pchip.rs
