Skip to main content

runmat_runtime/builtins/math/interpolation/
pchip.rs

1//! MATLAB-compatible `pchip` builtin for shape-preserving cubic interpolation.
2
3use runmat_builtins::{ResolveContext, Type, Value};
4use runmat_macros::runtime_builtin;
5
6use crate::builtins::common::spec::{
7    BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
8    ReductionNaN, ResidencyPolicy, ScalarType, ShapeRequirements,
9};
10
11use super::pp::{
12    build_pchip_pp, evaluate_pp, interp_error, pp_to_value, query_points, series_from_values,
13    Extrapolation,
14};
15
16const NAME: &str = "pchip";
17
18#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::interpolation::pchip")]
19pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
20    name: NAME,
21    op_kind: GpuOpKind::Custom("pchip"),
22    supported_precisions: &[ScalarType::F32, ScalarType::F64],
23    broadcast: BroadcastSemantics::Matlab,
24    provider_hooks: &[],
25    constant_strategy: ConstantStrategy::InlineLiteral,
26    residency: ResidencyPolicy::GatherImmediately,
27    nan_mode: ReductionNaN::Include,
28    two_pass_threshold: None,
29    workgroup_size: None,
30    accepts_nan_mode: false,
31    notes: "PCHIP coefficient construction and evaluation currently run on the CPU reference path after gathering GPU inputs.",
32};
33
34#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::interpolation::pchip")]
35pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
36    name: NAME,
37    shape: ShapeRequirements::Any,
38    constant_strategy: ConstantStrategy::InlineLiteral,
39    elementwise: None,
40    reduction: None,
41    emits_nan: false,
42    notes: "PCHIP builds a piecewise-polynomial representation and is not fused.",
43};
44
45fn pchip_type(args: &[Type], _ctx: &ResolveContext) -> Type {
46    match args.len() {
47        0 | 1 => Type::Unknown,
48        2 => Type::Struct { known_fields: None },
49        _ => match args.get(2) {
50            Some(Type::Num | Type::Int | Type::Bool) => Type::Num,
51            Some(Type::Tensor { shape }) | Some(Type::Logical { shape }) => Type::Tensor {
52                shape: shape.clone(),
53            },
54            _ => Type::tensor(),
55        },
56    }
57}
58
59#[runtime_builtin(
60    name = "pchip",
61    category = "math/interpolation",
62    summary = "Shape-preserving piecewise cubic Hermite interpolation.",
63    keywords = "pchip,shape preserving,cubic hermite,interpolation,ppval",
64    accel = "sink",
65    sink = true,
66    type_resolver(pchip_type),
67    builtin_path = "crate::builtins::math::interpolation::pchip"
68)]
69async fn pchip_builtin(x: Value, y: Value, rest: Vec<Value>) -> crate::BuiltinResult<Value> {
70    let series = series_from_values(x, y, NAME).await?;
71    let pp = build_pchip_pp(&series, NAME)?;
72    match rest.len() {
73        0 => pp_to_value(pp, NAME),
74        1 => {
75            let query = query_points(rest.into_iter().next().expect("query"), NAME).await?;
76            evaluate_pp(&pp, &query, &Extrapolation::Extrapolate, NAME)
77        }
78        _ => Err(interp_error(NAME, "pchip: too many input arguments")),
79    }
80}
81
82#[cfg(test)]
83mod tests {
84    use super::*;
85    use futures::executor::block_on;
86    use runmat_builtins::Tensor;
87
88    fn row(values: &[f64]) -> Value {
89        Value::Tensor(Tensor::new(values.to_vec(), vec![1, values.len()]).expect("tensor"))
90    }
91
92    #[test]
93    fn pchip_three_arg_evaluates_monotone_data() {
94        let value = block_on(pchip_builtin(
95            row(&[1.0, 2.0, 3.0, 4.0]),
96            row(&[0.0, 1.0, 1.5, 1.75]),
97            vec![row(&[1.5, 2.5, 3.5])],
98        ))
99        .expect("pchip");
100        let Value::Tensor(tensor) = value else {
101            panic!("expected tensor");
102        };
103        assert!(tensor.data.windows(2).all(|pair| pair[1] >= pair[0]));
104    }
105}