use runmat_builtins::{StructValue, Tensor, Value};
use crate::builtins::common::random_args::keyword_of;
use crate::builtins::common::tensor;
use crate::{build_runtime_error, dispatcher, BuiltinResult, RuntimeError};
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum InterpMethod {
Linear,
Nearest,
Spline,
Pchip,
}
#[derive(Clone, Debug)]
pub enum Extrapolation {
Nan,
Extrapolate,
Value(f64),
}
impl Default for Extrapolation {
fn default() -> Self {
Self::Nan
}
}
#[derive(Clone, Debug)]
pub struct NumericSeries {
pub x: Vec<f64>,
pub y: Vec<f64>,
pub series: usize,
pub trailing_shape: Vec<usize>,
}
#[derive(Clone, Debug)]
pub struct QueryPoints {
pub values: Vec<f64>,
pub shape: Vec<usize>,
}
#[derive(Clone, Debug)]
pub struct PiecewisePolynomial {
pub breaks: Vec<f64>,
pub coefs: Vec<f64>,
pub pieces: usize,
pub order: usize,
pub dim: usize,
}
pub fn interp_error(name: &'static str, message: impl Into<String>) -> RuntimeError {
build_runtime_error(message).with_builtin(name).build()
}
pub fn parse_method(value: &Value, name: &'static str) -> BuiltinResult<Option<InterpMethod>> {
let Some(keyword) = keyword_of(value) else {
return Ok(None);
};
let method = match keyword.as_str() {
"linear" => InterpMethod::Linear,
"nearest" => InterpMethod::Nearest,
"spline" => InterpMethod::Spline,
"pchip" | "cubic" => InterpMethod::Pchip,
_ => {
return Err(interp_error(
name,
format!("{name}: unsupported interpolation method '{keyword}'"),
))
}
};
Ok(Some(method))
}
pub async fn parse_extrapolation(
value: &Value,
name: &'static str,
) -> BuiltinResult<Option<Extrapolation>> {
if let Some(keyword) = keyword_of(value) {
return match keyword.as_str() {
"extrap" => Ok(Some(Extrapolation::Extrapolate)),
_ => Ok(None),
};
}
let gathered = dispatcher::gather_if_needed_async(value).await?;
let Some(scalar) = tensor::scalar_f64_from_value_async(&gathered)
.await
.map_err(|err| interp_error(name, format!("{name}: {err}")))?
else {
return Ok(None);
};
Ok(Some(Extrapolation::Value(scalar)))
}
pub async fn query_points(value: Value, name: &'static str) -> BuiltinResult<QueryPoints> {
let gathered = dispatcher::gather_if_needed_async(&value).await?;
let tensor = tensor::value_into_tensor_for(name, gathered)
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
let shape = canonical_shape(&tensor.shape, tensor.data.len());
Ok(QueryPoints {
values: tensor.data,
shape,
})
}
pub async fn vector_from_value(
value: Value,
label: &str,
name: &'static str,
) -> BuiltinResult<Vec<f64>> {
let gathered = dispatcher::gather_if_needed_async(&value).await?;
let tensor = tensor::value_into_tensor_for(name, gathered)
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
if !is_vector_shape(&tensor.shape) && tensor.data.len() > 1 {
return Err(interp_error(
name,
format!("{name}: {label} must be a vector"),
));
}
Ok(tensor.data)
}
pub async fn series_from_values(
x_value: Value,
y_value: Value,
name: &'static str,
) -> BuiltinResult<NumericSeries> {
let x = vector_from_value(x_value, "X", name).await?;
let y_host = dispatcher::gather_if_needed_async(&y_value).await?;
let y_tensor = tensor::value_into_tensor_for(name, y_host)
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
series_from_tensor(x, y_tensor, name)
}
pub async fn implicit_series_from_values(
y_value: Value,
name: &'static str,
) -> BuiltinResult<NumericSeries> {
let y_host = dispatcher::gather_if_needed_async(&y_value).await?;
let y_tensor = tensor::value_into_tensor_for(name, y_host)
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
let n = first_non_singleton_len(&y_tensor).unwrap_or(y_tensor.data.len());
let x: Vec<f64> = (1..=n).map(|value| value as f64).collect();
series_from_tensor(x, y_tensor, name)
}
fn series_from_tensor(
x: Vec<f64>,
y_tensor: Tensor,
name: &'static str,
) -> BuiltinResult<NumericSeries> {
validate_breaks(&x, name)?;
let n = x.len();
if y_tensor.data.len() != n && !y_tensor.data.len().is_multiple_of(n) {
return Err(interp_error(
name,
format!("{name}: Y length must match X or have X as its first dimension"),
));
}
let y_shape = canonical_shape(&y_tensor.shape, y_tensor.data.len());
let (series, trailing_shape) = if y_tensor.data.len() == n && is_vector_shape(&y_shape) {
(1, Vec::new())
} else if y_shape.first().copied() == Some(n) {
let series = y_tensor.data.len() / n;
let trailing = if y_shape.len() > 1 {
y_shape[1..].to_vec()
} else {
vec![series]
};
(series, trailing)
} else if y_tensor.data.len() == n {
(1, Vec::new())
} else {
return Err(interp_error(
name,
format!("{name}: size of Y must be compatible with X"),
));
};
Ok(NumericSeries {
x,
y: y_tensor.data,
series,
trailing_shape,
})
}
pub fn validate_breaks(x: &[f64], name: &'static str) -> BuiltinResult<()> {
if x.len() < 2 {
return Err(interp_error(
name,
format!("{name}: X must contain at least two points"),
));
}
if x.iter().any(|v| !v.is_finite()) {
return Err(interp_error(name, format!("{name}: X must be finite")));
}
for pair in x.windows(2) {
if pair[1] <= pair[0] {
return Err(interp_error(
name,
format!("{name}: X must be strictly increasing"),
));
}
}
Ok(())
}
pub fn evaluate_linear_or_nearest(
series: &NumericSeries,
query: &QueryPoints,
method: InterpMethod,
extrap: &Extrapolation,
name: &'static str,
) -> BuiltinResult<Value> {
let mut out = vec![0.0; query.values.len() * series.series];
for s in 0..series.series {
let y = &series.y[s * series.x.len()..(s + 1) * series.x.len()];
for (q_index, &xq) in query.values.iter().enumerate() {
out[q_index + s * query.values.len()] = match method {
InterpMethod::Linear => eval_linear(&series.x, y, xq, extrap),
InterpMethod::Nearest => eval_nearest(&series.x, y, xq, extrap),
_ => unreachable!("only direct methods are accepted here"),
};
}
}
tensor_from_query_output(out, query, series, name)
}
pub fn build_spline_pp(
series: &NumericSeries,
name: &'static str,
) -> BuiltinResult<PiecewisePolynomial> {
let n = series.x.len();
let pieces = n - 1;
let order = 4;
let rows = pieces * series.series;
let mut coefs = vec![0.0; rows * order];
for s in 0..series.series {
let y = &series.y[s * n..(s + 1) * n];
let local = spline_coefs_for_series(&series.x, y)?;
for piece in 0..pieces {
let row = s * pieces + piece;
for col in 0..order {
coefs[row + col * rows] = local[piece][col];
}
}
}
validate_finite_coefs(&coefs, name)?;
Ok(PiecewisePolynomial {
breaks: series.x.clone(),
coefs,
pieces,
order,
dim: series.series,
})
}
pub fn build_pchip_pp(
series: &NumericSeries,
name: &'static str,
) -> BuiltinResult<PiecewisePolynomial> {
let n = series.x.len();
let pieces = n - 1;
let order = 4;
let rows = pieces * series.series;
let mut coefs = vec![0.0; rows * order];
for s in 0..series.series {
let y = &series.y[s * n..(s + 1) * n];
let local = pchip_coefs_for_series(&series.x, y);
for piece in 0..pieces {
let row = s * pieces + piece;
for col in 0..order {
coefs[row + col * rows] = local[piece][col];
}
}
}
validate_finite_coefs(&coefs, name)?;
Ok(PiecewisePolynomial {
breaks: series.x.clone(),
coefs,
pieces,
order,
dim: series.series,
})
}
pub fn pp_to_value(pp: PiecewisePolynomial, name: &'static str) -> BuiltinResult<Value> {
let breaks = Tensor::new(pp.breaks.clone(), vec![1, pp.breaks.len()])
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
let coefs = Tensor::new(pp.coefs.clone(), vec![pp.pieces * pp.dim, pp.order])
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
let mut st = StructValue::new();
st.insert("form", Value::String("pp".to_string()));
st.insert("breaks", Value::Tensor(breaks));
st.insert("coefs", Value::Tensor(coefs));
st.insert("pieces", Value::Num(pp.pieces as f64));
st.insert("order", Value::Num(pp.order as f64));
st.insert("dim", Value::Num(pp.dim as f64));
Ok(Value::Struct(st))
}
pub async fn pp_from_value(value: Value, name: &'static str) -> BuiltinResult<PiecewisePolynomial> {
let gathered = dispatcher::gather_if_needed_async(&value).await?;
let Value::Struct(st) = gathered else {
return Err(interp_error(
name,
format!("{name}: first argument must be a pp struct"),
));
};
let form = st
.fields
.get("form")
.and_then(keyword_of)
.ok_or_else(|| interp_error(name, format!("{name}: pp struct is missing form")))?;
if form != "pp" {
return Err(interp_error(
name,
format!("{name}: struct form must be 'pp'"),
));
}
let breaks_value = st
.fields
.get("breaks")
.ok_or_else(|| interp_error(name, format!("{name}: pp struct is missing breaks")))?
.clone();
let breaks = vector_from_value(breaks_value, "breaks", name).await?;
validate_breaks(&breaks, name)?;
let coefs_value = st
.fields
.get("coefs")
.ok_or_else(|| interp_error(name, format!("{name}: pp struct is missing coefs")))?
.clone();
let coefs_host = dispatcher::gather_if_needed_async(&coefs_value).await?;
let coefs_tensor = tensor::value_into_tensor_for(name, coefs_host)
.map_err(|err| interp_error(name, format!("{name}: {err}")))?;
let pieces = scalar_field(&st, "pieces", name).await? as usize;
let order = scalar_field(&st, "order", name).await? as usize;
let dim = scalar_field(&st, "dim", name).await? as usize;
if pieces + 1 != breaks.len() {
return Err(interp_error(name, format!("{name}: malformed pp breaks")));
}
if order == 0 || dim == 0 {
return Err(interp_error(
name,
format!("{name}: malformed pp dimensions"),
));
}
if coefs_tensor.data.len() != pieces * dim * order {
return Err(interp_error(
name,
format!("{name}: malformed pp coefficients"),
));
}
Ok(PiecewisePolynomial {
breaks,
coefs: coefs_tensor.data,
pieces,
order,
dim,
})
}
async fn scalar_field(st: &StructValue, field: &str, name: &'static str) -> BuiltinResult<f64> {
let value = st
.fields
.get(field)
.ok_or_else(|| interp_error(name, format!("{name}: pp struct is missing {field}")))?;
let Some(raw) = tensor::scalar_f64_from_value_async(value)
.await
.map_err(|err| interp_error(name, format!("{name}: {err}")))?
else {
return Err(interp_error(
name,
format!("{name}: pp field {field} must be scalar"),
));
};
if !raw.is_finite() || raw < 1.0 || (raw.round() - raw).abs() > 1e-9 {
return Err(interp_error(
name,
format!("{name}: pp field {field} must be a positive integer"),
));
}
Ok(raw.round())
}
pub fn evaluate_pp(
pp: &PiecewisePolynomial,
query: &QueryPoints,
extrap: &Extrapolation,
name: &'static str,
) -> BuiltinResult<Value> {
let mut out = vec![0.0; query.values.len() * pp.dim];
for d in 0..pp.dim {
for (q_index, &xq) in query.values.iter().enumerate() {
out[q_index + d * query.values.len()] = eval_pp_scalar(pp, d, xq, extrap);
}
}
tensor_from_pp_output(out, query, pp.dim, name)
}
pub fn tensor_from_query_output(
data: Vec<f64>,
query: &QueryPoints,
series: &NumericSeries,
name: &'static str,
) -> BuiltinResult<Value> {
let shape = output_shape(&query.shape, series.series, &series.trailing_shape);
value_from_data(data, shape, name)
}
fn tensor_from_pp_output(
data: Vec<f64>,
query: &QueryPoints,
dim: usize,
name: &'static str,
) -> BuiltinResult<Value> {
let shape = if dim == 1 {
query.shape.clone()
} else {
let mut shape = query.shape.clone();
shape.push(dim);
shape
};
value_from_data(data, shape, name)
}
fn value_from_data(data: Vec<f64>, shape: Vec<usize>, name: &'static str) -> BuiltinResult<Value> {
if data.len() == 1 {
return Ok(Value::Num(data[0]));
}
let tensor =
Tensor::new(data, shape).map_err(|err| interp_error(name, format!("{name}: {err}")))?;
Ok(Value::Tensor(tensor))
}
fn output_shape(query_shape: &[usize], series: usize, trailing_shape: &[usize]) -> Vec<usize> {
if series == 1 {
return query_shape.to_vec();
}
let mut shape = query_shape.to_vec();
if trailing_shape.is_empty() {
shape.push(series);
} else {
shape.extend(trailing_shape.iter().copied());
}
shape
}
fn canonical_shape(shape: &[usize], len: usize) -> Vec<usize> {
tensor::default_shape_for(shape, len)
}
fn first_non_singleton_len(tensor: &Tensor) -> Option<usize> {
canonical_shape(&tensor.shape, tensor.data.len())
.into_iter()
.find(|&dim| dim > 1)
}
pub(super) fn is_vector_shape(shape: &[usize]) -> bool {
match shape {
[] => true,
[_] => true,
[rows, cols] => *rows == 1 || *cols == 1,
dims => dims.iter().filter(|&&dim| dim > 1).count() <= 1,
}
}
fn eval_linear(x: &[f64], y: &[f64], xq: f64, extrap: &Extrapolation) -> f64 {
if !xq.is_finite() {
return f64::NAN;
}
let Some(piece) = interval_index(x, xq, matches!(extrap, Extrapolation::Extrapolate)) else {
return out_of_range_value(extrap);
};
let h = x[piece + 1] - x[piece];
let t = (xq - x[piece]) / h;
y[piece] + t * (y[piece + 1] - y[piece])
}
fn eval_nearest(x: &[f64], y: &[f64], xq: f64, extrap: &Extrapolation) -> f64 {
if !xq.is_finite() {
return f64::NAN;
}
if xq < x[0] {
return match extrap {
Extrapolation::Extrapolate => y[0],
_ => out_of_range_value(extrap),
};
}
if xq > x[x.len() - 1] {
return match extrap {
Extrapolation::Extrapolate => y[y.len() - 1],
_ => out_of_range_value(extrap),
};
}
match x.binary_search_by(|probe| probe.partial_cmp(&xq).unwrap()) {
Ok(index) => y[index],
Err(index) => {
let left = index.saturating_sub(1);
let right = index.min(x.len() - 1);
if (xq - x[left]).abs() <= (x[right] - xq).abs() {
y[left]
} else {
y[right]
}
}
}
}
fn eval_pp_scalar(pp: &PiecewisePolynomial, series: usize, xq: f64, extrap: &Extrapolation) -> f64 {
if !xq.is_finite() {
return f64::NAN;
}
let Some(piece) = interval_index(&pp.breaks, xq, matches!(extrap, Extrapolation::Extrapolate))
else {
return out_of_range_value(extrap);
};
let t = xq - pp.breaks[piece];
let rows = pp.pieces * pp.dim;
let row = series * pp.pieces + piece;
let mut acc = 0.0;
for col in 0..pp.order {
acc = acc * t + pp.coefs[row + col * rows];
}
acc
}
pub(super) fn interval_index(x: &[f64], xq: f64, allow_extrapolation: bool) -> Option<usize> {
if xq < x[0] {
return allow_extrapolation.then_some(0);
}
let last = x.len() - 1;
if xq > x[last] {
return allow_extrapolation.then_some(last - 1);
}
if xq == x[last] {
return Some(last - 1);
}
match x.binary_search_by(|probe| probe.partial_cmp(&xq).unwrap()) {
Ok(index) => Some(index.min(last - 1)),
Err(index) if index > 0 && index < x.len() => Some(index - 1),
_ => None,
}
}
pub(super) fn out_of_range_value(extrap: &Extrapolation) -> f64 {
match extrap {
Extrapolation::Nan => f64::NAN,
Extrapolation::Extrapolate => f64::NAN,
Extrapolation::Value(value) => *value,
}
}
fn spline_coefs_for_series(x: &[f64], y: &[f64]) -> BuiltinResult<Vec<[f64; 4]>> {
let n = x.len();
if n == 2 {
let h = x[1] - x[0];
return Ok(vec![[0.0, 0.0, (y[1] - y[0]) / h, y[0]]]);
}
if n == 3 {
return Ok(quadratic_as_piecewise_cubic(x, y));
}
let h: Vec<f64> = x.windows(2).map(|pair| pair[1] - pair[0]).collect();
let mut a = vec![vec![0.0; n]; n];
let mut rhs = vec![0.0; n];
a[0][0] = -h[1];
a[0][1] = h[0] + h[1];
a[0][2] = -h[0];
for i in 1..(n - 1) {
a[i][i - 1] = h[i - 1];
a[i][i] = 2.0 * (h[i - 1] + h[i]);
a[i][i + 1] = h[i];
rhs[i] = 6.0 * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]);
}
a[n - 1][n - 3] = -h[n - 2];
a[n - 1][n - 2] = h[n - 3] + h[n - 2];
a[n - 1][n - 1] = -h[n - 3];
let second = solve_dense(a, rhs)?;
Ok(second_derivatives_to_coefs(x, y, &second))
}
fn quadratic_as_piecewise_cubic(x: &[f64], y: &[f64]) -> Vec<[f64; 4]> {
let x0 = x[0];
let x1 = x[1];
let x2 = x[2];
let y0 = y[0];
let y1 = y[1];
let y2 = y[2];
let c2 =
y0 / ((x0 - x1) * (x0 - x2)) + y1 / ((x1 - x0) * (x1 - x2)) + y2 / ((x2 - x0) * (x2 - x1));
let c1 = -y0 * (x1 + x2) / ((x0 - x1) * (x0 - x2))
- y1 * (x0 + x2) / ((x1 - x0) * (x1 - x2))
- y2 * (x0 + x1) / ((x2 - x0) * (x2 - x1));
let c0 = y0 * x1 * x2 / ((x0 - x1) * (x0 - x2))
+ y1 * x0 * x2 / ((x1 - x0) * (x1 - x2))
+ y2 * x0 * x1 / ((x2 - x0) * (x2 - x1));
(0..2)
.map(|i| {
let xi = x[i];
let a = c2 * xi * xi + c1 * xi + c0;
let b = 2.0 * c2 * xi + c1;
[0.0, c2, b, a]
})
.collect()
}
fn second_derivatives_to_coefs(x: &[f64], y: &[f64], second: &[f64]) -> Vec<[f64; 4]> {
let mut coefs = Vec::with_capacity(x.len() - 1);
for i in 0..(x.len() - 1) {
let h = x[i + 1] - x[i];
let a = y[i];
let b = (y[i + 1] - y[i]) / h - h * (2.0 * second[i] + second[i + 1]) / 6.0;
let c = second[i] / 2.0;
let d = (second[i + 1] - second[i]) / (6.0 * h);
coefs.push([d, c, b, a]);
}
coefs
}
fn pchip_coefs_for_series(x: &[f64], y: &[f64]) -> Vec<[f64; 4]> {
let n = x.len();
let h: Vec<f64> = x.windows(2).map(|pair| pair[1] - pair[0]).collect();
let delta: Vec<f64> = y
.windows(2)
.zip(h.iter())
.map(|(pair, &width)| (pair[1] - pair[0]) / width)
.collect();
let mut slopes = vec![0.0; n];
if n == 2 {
slopes[0] = delta[0];
slopes[1] = delta[0];
} else {
slopes[0] = endpoint_slope(h[0], h[1], delta[0], delta[1]);
slopes[n - 1] = endpoint_slope(h[n - 2], h[n - 3], delta[n - 2], delta[n - 3]);
for i in 1..(n - 1) {
if delta[i - 1] == 0.0 || delta[i] == 0.0 || delta[i - 1].signum() != delta[i].signum()
{
slopes[i] = 0.0;
} else {
let w1 = 2.0 * h[i] + h[i - 1];
let w2 = h[i] + 2.0 * h[i - 1];
slopes[i] = (w1 + w2) / (w1 / delta[i - 1] + w2 / delta[i]);
}
}
}
let mut coefs = Vec::with_capacity(n - 1);
for i in 0..(n - 1) {
let c0 = y[i];
let c1 = slopes[i];
let c2 = (3.0 * delta[i] - 2.0 * slopes[i] - slopes[i + 1]) / h[i];
let c3 = (slopes[i] + slopes[i + 1] - 2.0 * delta[i]) / (h[i] * h[i]);
coefs.push([c3, c2, c1, c0]);
}
coefs
}
fn endpoint_slope(h0: f64, h1: f64, del0: f64, del1: f64) -> f64 {
let mut slope = ((2.0 * h0 + h1) * del0 - h0 * del1) / (h0 + h1);
if slope.signum() != del0.signum() {
slope = 0.0;
} else if del0.signum() != del1.signum() && slope.abs() > (3.0 * del0).abs() {
slope = 3.0 * del0;
}
slope
}
fn solve_dense(mut a: Vec<Vec<f64>>, mut rhs: Vec<f64>) -> BuiltinResult<Vec<f64>> {
let n = rhs.len();
for col in 0..n {
let mut pivot = col;
let mut pivot_abs = a[col][col].abs();
for row in (col + 1)..n {
let candidate = a[row][col].abs();
if candidate > pivot_abs {
pivot = row;
pivot_abs = candidate;
}
}
if pivot_abs <= 1.0e-14 {
return Err(interp_error(
"spline",
"spline: singular interpolation system",
));
}
if pivot != col {
a.swap(pivot, col);
rhs.swap(pivot, col);
}
for row in (col + 1)..n {
let factor = a[row][col] / a[col][col];
a[row][col] = 0.0;
for k in (col + 1)..n {
a[row][k] -= factor * a[col][k];
}
rhs[row] -= factor * rhs[col];
}
}
let mut solution = vec![0.0; n];
for row in (0..n).rev() {
let mut acc = rhs[row];
for col in (row + 1)..n {
acc -= a[row][col] * solution[col];
}
solution[row] = acc / a[row][row];
}
Ok(solution)
}
fn validate_finite_coefs(coefs: &[f64], name: &'static str) -> BuiltinResult<()> {
if coefs.iter().any(|value| !value.is_finite()) {
return Err(interp_error(
name,
format!("{name}: interpolation coefficients must be finite"),
));
}
Ok(())
}
#[cfg(test)]
mod tests {
use super::*;
fn scalar_series(y: &[f64]) -> NumericSeries {
NumericSeries {
x: (1..=y.len()).map(|v| v as f64).collect(),
y: y.to_vec(),
series: 1,
trailing_shape: Vec::new(),
}
}
#[test]
fn pchip_preserves_monotone_samples() {
let series = scalar_series(&[0.0, 1.0, 1.5, 1.75]);
let pp = build_pchip_pp(&series, "pchip").expect("pchip");
let query = QueryPoints {
values: vec![1.25, 1.5, 2.5, 3.5],
shape: vec![1, 4],
};
let value = evaluate_pp(&pp, &query, &Extrapolation::Extrapolate, "pchip").expect("ppval");
let Value::Tensor(tensor) = value else {
panic!("expected tensor");
};
assert!(tensor.data.windows(2).all(|pair| pair[1] >= pair[0]));
}
#[test]
fn spline_matches_quadratic_for_three_points() {
let series = scalar_series(&[1.0, 4.0, 9.0]);
let pp = build_spline_pp(&series, "spline").expect("spline");
let query = QueryPoints {
values: vec![1.5, 2.5],
shape: vec![1, 2],
};
let value = evaluate_pp(&pp, &query, &Extrapolation::Extrapolate, "spline").expect("ppval");
let Value::Tensor(tensor) = value else {
panic!("expected tensor");
};
assert!((tensor.data[0] - 2.25).abs() < 1e-10);
assert!((tensor.data[1] - 6.25).abs() < 1e-10);
}
#[test]
fn ppval_multiseries_appends_series_dimension() {
let series = NumericSeries {
x: vec![1.0, 2.0, 3.0],
y: vec![1.0, 4.0, 9.0, 100.0, 400.0, 900.0],
series: 2,
trailing_shape: vec![2],
};
let pp = build_spline_pp(&series, "spline").expect("spline");
let query = QueryPoints {
values: vec![1.5, 2.5],
shape: vec![1, 2],
};
let value = evaluate_pp(&pp, &query, &Extrapolation::Extrapolate, "spline").expect("ppval");
let Value::Tensor(tensor) = value else {
panic!("expected tensor");
};
assert_eq!(tensor.shape, vec![1, 2, 2]);
assert_eq!(tensor.data.len(), 4);
assert!((tensor.data[0] - 2.25).abs() < 1e-10);
assert!((tensor.data[1] - 6.25).abs() < 1e-10);
assert!((tensor.data[2] - 225.0).abs() < 1e-10);
assert!((tensor.data[3] - 625.0).abs() < 1e-10);
}
}