diff --git a/src/interpolator/enums.rs b/src/interpolator/enums.rs index 42249d4..905589b 100644 --- a/src/interpolator/enums.rs +++ b/src/interpolator/enums.rs @@ -386,7 +386,7 @@ where Interp1D(Interp1DBase>), Interp2D(Interp2DBase>), Interp3D(Interp3DBase>), - InterpND(InterpNDBase>), + InterpND(InterpNDBase), } /// Owned interpolator enum (see [`InterpolatorEnumBase`] for the generic form). pub type InterpolatorEnum = InterpolatorEnumBase>; @@ -493,14 +493,14 @@ where } } -impl From> for InterpolatorEnumBase +impl From> for InterpolatorEnumBase where D: Data + RawDataClone + Clone, D::Elem: Float + Debug, - S: Into> + Clone, + S: Into + Clone, { #[inline] - fn from(interpolator: InterpNDBase) -> Self { + fn from(interpolator: InterpNDBase) -> Self { InterpolatorEnumBase::InterpND(InterpNDBase { data: interpolator.data, strategy: interpolator.strategy.into(), @@ -579,7 +579,7 @@ where pub fn new_nd( grid: Vec>, values: ArrayBase, - strategy: impl Into>, + strategy: impl Into, extrapolate: Extrapolate, ) -> Result { Ok(Self::InterpND(InterpNDBase::new( @@ -759,7 +759,7 @@ mod tests { .unwrap(); let interp1 = InterpolatorEnumBase::from(interp0.clone()); - let interp2: InterpNDBase<_, strategy::enums::StrategyNDEnum> = InterpNDBase::new( + let interp2: InterpNDBase<_, _, strategy::enums::StrategyNDEnum> = InterpNDBase::new( vec![x.view(), y.view()], f_xy_dyn.view(), strategy::Nearest.into(), diff --git a/src/interpolator/n/mod.rs b/src/interpolator/n/mod.rs index c04e16e..71b3ae3 100644 --- a/src/interpolator/n/mod.rs +++ b/src/interpolator/n/mod.rs @@ -10,6 +10,12 @@ mod tests; /// Interpolator data for N-dimensional interpolators, where N can vary at runtime. /// +/// Split into a grid storage type `Dg` and a value storage type `Dv` (ninterp's +/// Tg/Tv split, prototyped here first, see issue #57). `InterpDataND` / +/// `InterpND` cover the common same-type case (`Dg::Elem == Dv::Elem == T`); +/// reach for `InterpDataNDBase, OwnedRepr>` directly when they +/// should genuinely differ. +/// /// See [`InterpDataBase`] and its dimension-specific aliases /// for concrete-dimensionality interpolator data structs. #[derive(Debug, Clone)] @@ -17,35 +23,43 @@ mod tests; #[cfg_attr( feature = "serde", serde(bound( - serialize = "D::Elem: Serialize", + serialize = "Dg::Elem: Serialize, Dv::Elem: Serialize", deserialize = " - D: DataOwned, - D::Elem: Deserialize<'de>, + Dg: DataOwned, + Dg::Elem: Deserialize<'de>, + Dv: DataOwned, + Dv::Elem: Deserialize<'de>, " )) )] -pub struct InterpDataNDBase +pub struct InterpDataNDBase where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, { - /// Coordinate grid: a vector of 1-dimensional [`ArrayBase`]. + /// Coordinate grid: a vector of 1-dimensional [`ArrayBase`]. #[cfg_attr(feature = "serde", serde(deserialize_with = "deserialize_grid_vec"))] - pub grid: Vec>, + pub grid: Vec>, /// Function values at coordinates: a single dynamic-dimensional [`ArrayBase`]. #[cfg_attr(feature = "serde", serde(deserialize_with = "deserialize_dyn"))] - pub values: ArrayBase, + pub values: ArrayBase, } -/// Owned data variant for N-D data (see [`InterpDataNDBase`] for the generic form). -pub type InterpDataND = InterpDataNDBase>; -/// Viewed data variant for N-D data (see [`InterpDataNDBase`] for the generic form). -pub type InterpDataNDView = InterpDataNDBase>; +/// Owned data variant for N-D data, same-type convenience (see [`InterpDataNDBase`] +/// for the general grid/value-split form). +pub type InterpDataND = InterpDataNDBase, OwnedRepr>; +/// Viewed data variant for N-D data, same-type convenience (see [`InterpDataNDBase`] +/// for the general grid/value-split form). +pub type InterpDataNDView = InterpDataNDBase, ViewRepr>; #[cfg(feature = "serde")] -impl SerializeNested for InterpDataNDBase +impl SerializeNested for InterpDataNDBase where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug + Serialize, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug + Serialize, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug + Serialize, { fn serialize_nested(&self, serializer: S) -> Result where @@ -58,29 +72,33 @@ where } } -impl PartialEq for InterpDataNDBase +impl PartialEq for InterpDataNDBase where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, - ArrayBase: PartialEq, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, + ArrayBase: PartialEq, { fn eq(&self, other: &Self) -> bool { self.grid == other.grid && self.values == other.values } } -impl InterpDataNDBase +impl InterpDataNDBase where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, { /// Construct and validate a new [`InterpDataND`]. pub fn new( - grid: Vec>, - values: ArrayBase, + grid: Vec>, + values: ArrayBase, ) -> Result where - D::Elem: PartialOrd, + Dg::Elem: PartialOrd, { let data = Self { grid, values }; data.validate()?; @@ -90,7 +108,7 @@ where /// Validate interpolator data. pub fn validate(&self) -> Result<(), ValidateError> where - D::Elem: PartialOrd, + Dg::Elem: PartialOrd, { let n = self.ndim(); if (self.grid.len() != n) && !(n == 0 && self.grid.iter().all(|g| g.is_empty())) { @@ -131,19 +149,20 @@ where } /// View interpolator data. - pub fn view(&self) -> InterpDataNDView<&D::Elem> { - InterpDataNDView { + pub fn view(&self) -> InterpDataNDBase, ViewRepr<&Dv::Elem>> { + InterpDataNDBase { grid: self.grid.iter().map(|g| g.view()).collect(), values: self.values.view(), } } /// Turn the data into an [`InterpDataND`], cloning the array elements if necessary. - pub fn into_owned(self) -> InterpDataND + pub fn into_owned(self) -> InterpDataNDBase, OwnedRepr> where - D::Elem: Clone, + Dg::Elem: Clone, + Dv::Elem: Clone, { - InterpDataND { + InterpDataNDBase { grid: self.grid.into_iter().map(|g| g.into_owned()).collect(), values: self.values.into_owned(), } @@ -151,55 +170,143 @@ where } /// N-D interpolator +/// +/// Split into a grid storage type `Dg` and a value storage type `Dv` (ninterp's +/// Tg/Tv split, prototyped here first, see issue #57). `InterpND` covers the +/// common same-type case (`Dg::Elem == Dv::Elem == T`), and is the only shape that +/// implements [`Interpolator`]: the general `Dg != Dv` case is reached through +/// the inherent methods below instead, since `Interpolator`'s single type +/// parameter can't express a split grid/value type. Reach for +/// `InterpNDBase, OwnedRepr, S>` directly when they should +/// genuinely differ. #[derive(Debug, Clone)] #[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] #[cfg_attr( feature = "serde", serde(bound( serialize = " - D::Elem: Serialize, + Dg::Elem: Serialize, + Dv::Elem: Serialize, S: Serialize, ", deserialize = " - D: DataOwned, - D::Elem: Deserialize<'de>, + Dg: DataOwned, + Dg::Elem: Deserialize<'de>, + Dv: DataOwned, + Dv::Elem: Deserialize<'de>, S: Deserialize<'de>, " )) )] -pub struct InterpNDBase +pub struct InterpNDBase where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, S: Clone, { /// Interpolator data. - pub data: InterpDataNDBase, + pub data: InterpDataNDBase, /// Interpolation strategy. pub strategy: S, /// Extrapolation setting. - pub extrapolate: Extrapolate, + pub extrapolate: Extrapolate, +} +/// Owned interpolator variant, same-type convenience (see [`InterpNDBase`] for the +/// general grid/value-split form). +pub type InterpND = InterpNDBase, OwnedRepr, S>; +/// Viewed interpolator variant, same-type convenience (see [`InterpNDBase`] for the +/// general grid/value-split form). +pub type InterpNDView = InterpNDBase, ViewRepr, S>; + +// The shared macros below (`partialeq_impl!`, `serialize_nested_impl!`, +// `extrapolate_impl!`, `set_strategy_box_impl!`, `set_strategy_enum_impl!`, +// `any_interpolator_impl!`) are single-`D` and shared with `Interp1D`/`2D`/`3D`, so +// reusing them here would force a signature change onto types this prototype isn't +// touching. Left commented out (not deleted): when `Interp1D`/`2D`/`3D` eventually +// follow the same Tg/Tv split, the macros themselves get updated to take `Dg`/`Dv`, +// and these call sites are what make that swap-back a diff instead of a rewrite. +// +// partialeq_impl!(InterpNDBase, InterpDataNDBase, StrategyND); +// #[cfg(feature = "serde")] +// serialize_nested_impl!(InterpNDBase, InterpDataNDBase, StrategyND); +// extrapolate_impl!(InterpNDBase, StrategyND); +// set_strategy_box_impl!(InterpNDBase, StrategyND); +// set_strategy_enum_impl!( +// InterpNDBase, +// strategy::enums::StrategyNDEnum, +// strategy::enums::StrategyNDEnum +// ); +// any_interpolator_impl!(InterpND, StrategyND); + +impl PartialEq for InterpNDBase +where + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, + S: Clone + PartialEq, + InterpDataNDBase: PartialEq, +{ + fn eq(&self, other: &Self) -> bool { + self.data == other.data + && self.strategy == other.strategy + && self.extrapolate == other.extrapolate + } } -/// Owned interpolator variant (see [`InterpNDBase`] for the generic form). -pub type InterpND = InterpNDBase, S>; -/// Viewed interpolator variant (see [`InterpNDBase`] for the generic form). -pub type InterpNDView = InterpNDBase, S>; -partialeq_impl!(InterpNDBase, InterpDataNDBase, StrategyND); #[cfg(feature = "serde")] -serialize_nested_impl!(InterpNDBase, InterpDataNDBase, StrategyND); +impl SerializeNested for InterpNDBase +where + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug + Serialize, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug + Serialize, + S: Clone + Serialize, + InterpDataNDBase: SerializeNested + Serialize, +{ + fn serialize_nested(&self, serializer: Ser) -> Result + where + Ser: Serializer, + { + let mut s = serializer.serialize_struct("InterpNDBase", 3)?; + s.serialize_field("data", &Nested(&self.data))?; + s.serialize_field("strategy", &self.strategy)?; + s.serialize_field("extrapolate", &self.extrapolate)?; + s.end() + } +} -impl InterpNDBase +impl InterpNDBase where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, - S: StrategyND + Clone, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, + S: StrategyND + Clone, { + /// Check that `extrapolate` is applicable to the current strategy. + /// + /// Only [`Extrapolate::Enable`] can be rejected, and only by a strategy whose + /// `allow_extrapolate` returns `false`. Takes the setting as an argument rather + /// than reading `self.extrapolate`, so `set_extrapolate` can vet a candidate + /// before storing it. + pub fn validate_extrapolate( + &self, + extrapolate: &Extrapolate, + ) -> Result<(), ValidateError> { + if matches!(extrapolate, Extrapolate::Enable) && !self.strategy.allow_extrapolate() { + return Err(ValidateError::ExtrapolateUnsupported); + } + Ok(()) + } + /// Re-run the strategy's [`StrategyND::validate`] against the current data. /// - /// `new`, `set_strategy`, and [`Interpolator::validate`] already call this - /// internally, so this is only needed after mutating the public `data`/`strategy` - /// fields directly. + /// `new`, `set_strategy`, and [`Self::validate`] already call this internally, + /// so this is only needed after mutating the public `data`/`strategy` fields + /// directly. pub fn validate_strategy(&self) -> Result<(), ValidateError> { self.strategy.validate(&self.data) } @@ -209,20 +316,77 @@ where /// `new` and `set_strategy` already call this internally, so this is only needed /// after bypassing them: mutating the public `data`/`strategy` fields directly, or /// deserializing an interpolator whose strategy skips its cached state from - /// serialization (e.g. via `#[serde(skip)]`, to avoid bloating the wire format with - /// a large derived array). `Deserialize` does not call `init`; if the cached state - /// is instead stored in ordinary serialized fields, it comes back as-is and this + /// serialization. `Deserialize` does not call `init`; if the cached state is + /// instead stored in ordinary serialized fields, it comes back as-is and this /// isn't needed. pub fn init_strategy(&mut self) -> Result<(), ValidateError> { self.strategy.init(&self.data) } } -impl InterpNDBase +impl InterpNDBase>> where - D: Data + RawDataClone + Clone, - D::Elem: PartialOrd + Debug, - S: StrategyND + Clone, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, +{ + /// Update strategy at runtime, calling [`StrategyND::init`] on the new strategy + /// against the current data. + /// + /// To swap in a strategy without re-running `init` (e.g. one whose state was + /// already established elsewhere), assign the `strategy` field directly instead. + pub fn set_strategy( + &mut self, + strategy: Box>, + ) -> Result<(), ValidateError> { + self.strategy = strategy; + self.validate_extrapolate(&self.extrapolate)?; + self.validate_strategy()?; + self.init_strategy() + } +} + +impl InterpNDBase +where + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug + NumCast + PartialOrd + Copy, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug + NumCast + Copy, +{ + /// Update strategy at runtime, calling [`strategy::enums::StrategyNDEnum::init`] on + /// the new strategy against the current data. + /// + /// To swap in a strategy without re-running `init` (e.g. one whose state was + /// already established elsewhere), assign the `strategy` field directly instead. + pub fn set_strategy( + &mut self, + strategy: impl Into, + ) -> Result<(), ValidateError> { + self.strategy = strategy.into(); + self.validate_extrapolate(&self.extrapolate)?; + self.validate_strategy()?; + self.init_strategy() + } +} + +impl AnyInterpolator for InterpND +where + T: Float + Euclid + Debug + Send + Sync + 'static, + S: StrategyND, OwnedRepr> + Clone + Send + Sync + 'static, +{ + fn as_any(&self) -> &dyn Any { + self + } +} + +impl InterpNDBase +where + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialOrd + Debug + NumCast + Copy, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug + Copy, + S: StrategyND + Clone, { /// Construct and validate an N-D (any dimensionality) interpolator. /// @@ -267,10 +431,10 @@ where /// )); /// ``` pub fn new( - grid: Vec>, - values: ArrayBase, + grid: Vec>, + values: ArrayBase, strategy: S, - extrapolate: Extrapolate, + extrapolate: Extrapolate, ) -> Result { let mut interpolator = Self { data: InterpDataNDBase::new(grid, values)?, @@ -283,52 +447,49 @@ where Ok(interpolator) } - /// Return an interpolator with viewed data. - pub fn view(&self) -> InterpNDView<&D::Elem, S> - where - S: for<'a> StrategyND>, - D::Elem: Clone, - { - InterpNDView { - data: self.data.view(), - strategy: self.strategy.clone(), - extrapolate: self.extrapolate.clone(), - } - } - - /// Turn the interpolator into an [`InterpND`], cloning the array elements if necessary. - pub fn into_owned(self) -> InterpND - where - S: StrategyND>, - D::Elem: Clone, - { - InterpND { - data: self.data.into_owned(), - strategy: self.strategy.clone(), - extrapolate: self.extrapolate.clone(), - } - } -} - -impl Interpolator for InterpNDBase -where - D: Data + RawDataClone + Clone, - D::Elem: Num + PartialOrd + Euclid + Copy + Debug, - S: StrategyND + Clone, -{ + /// Get data dimensionality. #[inline] - fn ndim(&self) -> usize { + pub fn ndim(&self) -> usize { self.data.ndim() } - fn validate(&self) -> Result<(), ValidateError> { + /// Validate interpolator data. + pub fn validate(&self) -> Result<(), ValidateError> { self.validate_extrapolate(&self.extrapolate)?; self.data.validate()?; self.validate_strategy()?; Ok(()) } - fn interpolate(&self, point: &[D::Elem]) -> Result { + /// Set [`Extrapolate`] variant, checking validity. + pub fn set_extrapolate( + &mut self, + extrapolate: Extrapolate, + ) -> Result<(), ValidateError> { + self.validate_extrapolate(&extrapolate)?; + self.extrapolate = extrapolate; + Ok(()) + } + + /// Casts `data.grid[dim]`'s bounds to `f64`, to compare against a `point` + /// (always `f64`-typed; see [`StrategyND`]'s docs for why). + fn grid_bounds_f64(&self, dim: usize) -> (f64, f64) { + let lo: f64 = num_traits::cast(*self.data.grid[dim].first().unwrap()) + .expect("grid element must cast to f64"); + let hi: f64 = num_traits::cast(*self.data.grid[dim].last().unwrap()) + .expect("grid element must cast to f64"); + (lo, hi) + } + + /// Interpolate at the supplied point. + /// + /// Named `interpolate_f64` (not `interpolate`) to avoid colliding with + /// [`Interpolator::interpolate`]'s `&[T]`-typed method of the same name on + /// `InterpND` (`Dg = Dv = T`): an inherent method always wins over a + /// trait method in dot-call resolution regardless of whether argument types + /// actually match, so a same-named `&[f64]`-typed inherent method here would + /// silently break `Interpolator::interpolate` for any `T != f64`. + pub fn interpolate_f64(&self, point: &[f64]) -> Result { let n = self.ndim(); if point.len() != n { return Err(InterpolateError::PointLength { @@ -341,22 +502,18 @@ where } let mut errors = Vec::new(); for dim in 0..n { - if !(self.data.grid[dim].first().unwrap()..=self.data.grid[dim].last().unwrap()) - .contains(&&point[dim]) - { + let (lo, hi) = self.grid_bounds_f64(dim); + if !(lo..=hi).contains(&point[dim]) { match &self.extrapolate { Extrapolate::Enable => {} Extrapolate::Fill(value) => return Ok(*value), Extrapolate::Clamp => { - let clamped_point: Vec<_> = point + let clamped_point: Vec = point .iter() .enumerate() - .map(|(dim, pt)| { - *clamp( - pt, - self.data.grid[dim].first().unwrap(), - self.data.grid[dim].last().unwrap(), - ) + .map(|(dim, &pt)| { + let (lo, hi) = self.grid_bounds_f64(dim); + clamp(pt, lo, hi) }) .collect(); return self.strategy.interpolate(&self.data, &clamped_point); @@ -376,13 +533,15 @@ where self.strategy.interpolate(&self.data, point) } - fn set_extrapolate(&mut self, extrapolate: Extrapolate) -> Result<(), ValidateError> { - self.validate_extrapolate(&extrapolate)?; - self.extrapolate = extrapolate; - Ok(()) - } - - fn interpolate_fast(&self, point: &[D::Elem]) -> D::Elem { + /// Interpolate without bounds/extrapolation checks, for use in hot loops where the + /// caller has already checked bounds and knows that extrapolation is not needed. + /// + /// Named `interpolate_f64_fast`, not `interpolate_fast`; see + /// [`Self::interpolate_f64`]'s doc for why. + /// + /// # Panics + /// Panics if `point.len()` doesn't match [`Self::ndim`]. + pub fn interpolate_f64_fast(&self, point: &[f64]) -> Dv::Elem { assert_eq!( point.len(), self.ndim(), @@ -391,10 +550,18 @@ where self.strategy.interpolate_fast(&self.data, point) } - fn batch_interpolate_into( + /// Interpolate at each of several points, writing results into `out` instead of + /// allocating. + /// + /// `self.extrapolate` is one setting for the whole call, not resolved per point: + /// every point still funnels into at most one call to the strategy. + /// + /// Named `batch_interpolate_f64_into`, not `batch_interpolate_into`; see + /// [`Self::interpolate_f64`]'s doc for why. + pub fn batch_interpolate_f64_into( &self, - points: &[&[D::Elem]], - out: &mut [D::Elem], + points: &[&[f64]], + out: &mut [Dv::Elem], ) -> Result<(), InterpolateError> { let n = self.ndim(); let failures: Vec = points @@ -425,53 +592,27 @@ where Extrapolate::Clamp => { // Clamping an in-bounds point is already identity, so every point // can be clamped unconditionally. - let clamped: Vec> = points + let clamped: Vec> = points .iter() .map(|&point| { point .iter() .enumerate() - .map(|(dim, pt)| { - *clamp( - pt, - self.data.grid[dim].first().unwrap(), - self.data.grid[dim].last().unwrap(), - ) + .map(|(dim, &pt)| { + let (lo, hi) = self.grid_bounds_f64(dim); + clamp(pt, lo, hi) }) .collect() }) .collect(); - let clamped: Vec<&[D::Elem]> = clamped.iter().map(Vec::as_slice).collect(); + let clamped: Vec<&[f64]> = clamped.iter().map(Vec::as_slice).collect(); self.strategy .batch_interpolate_into(&self.data, &clamped, out) } Extrapolate::Wrap => { - // Must go through `interpolate_wrapped`, not a raw-space `wrap()` - // here plus the checked `batch_interpolate_into`: a strategy whose - // working coordinate space differs from the grid's (e.g. - // `GridTransform`) needs to wrap in its own space, since wrapping - // doesn't commute with a nonlinear transform. Only out-of-bounds - // points go through it; `wrap()` isn't identity exactly at the - // boundary, so an in-bounds point must still take the plain - // `interpolate` path. - // - // Pre-scans the whole batch for domain violations before doing any - // actual interpolation work, so a nested `GridTransform`'s - // aggregation isn't lost to this per-point dispatch (it can't use - // `batch_interpolate_into` directly: which points wrap and which - // interpolate normally is decided per point). self.strategy.check_batch_domain(points)?; - // TODO: every point here calls `interpolate`/`interpolate_wrapped` - // one at a time rather than routing the in-bounds subset through - // `batch_interpolate_into`. No strategy shipped in this crate - // overrides `batch_interpolate_into` for real amortized work, so - // this costs nothing today; if one ever does, partition points into - // in-bounds/out-of-bounds first (like the `Fill` arm above) and - // route the in-bounds subset through `batch_interpolate_into` to - // pick it up. Not worth the added trait-surface complexity until - // there's a strategy that would actually benefit. for (o, &point) in out.iter_mut().zip(points) { - *o = if out_of_bounds(&self.data.grid, point) { + *o = if self.point_out_of_bounds(point) { self.strategy.interpolate_wrapped(&self.data, point)? } else { self.strategy.interpolate(&self.data, point)? @@ -486,15 +627,18 @@ where *o = *value; } let mut in_bounds_indices = Vec::new(); - let mut in_bounds_points: Vec<&[D::Elem]> = Vec::new(); + let mut in_bounds_points: Vec<&[f64]> = Vec::new(); for (i, &point) in points.iter().enumerate() { - if !out_of_bounds(&self.data.grid, point) { + if !self.point_out_of_bounds(point) { in_bounds_indices.push(i); in_bounds_points.push(point); } } if !in_bounds_indices.is_empty() { - let mut scratch = vec![D::Elem::zero(); in_bounds_indices.len()]; + let mut scratch: Vec = in_bounds_indices + .iter() + .map(|_| out[in_bounds_indices[0]]) + .collect(); self.strategy.batch_interpolate_into( &self.data, &in_bounds_points, @@ -508,12 +652,12 @@ where } Extrapolate::Error => { let mut errors = Vec::new(); - let mut in_bounds_points: Vec<&[D::Elem]> = Vec::new(); + let mut in_bounds_points: Vec<&[f64]> = Vec::new(); for (i, &point) in points.iter().enumerate() { let mut point_errors = Vec::new(); - for (dim, (axis, &coord)) in self.data.grid.iter().zip(point.iter()).enumerate() - { - if !(axis.first().unwrap()..=axis.last().unwrap()).contains(&&coord) { + for (dim, &pt) in point.iter().enumerate() { + let (lo, hi) = self.grid_bounds_f64(dim); + if !(lo..=hi).contains(&pt) { point_errors.push(OutOfBoundsAt { index: i, dim }); } } @@ -532,28 +676,166 @@ where } } - fn batch_interpolate_fast_into(&self, points: &[&[D::Elem]], out: &mut [D::Elem]) { + /// Is `point` out of `self.data.grid`'s bounds in any dimension? + fn point_out_of_bounds(&self, point: &[f64]) -> bool { + (0..point.len()).any(|dim| { + let (lo, hi) = self.grid_bounds_f64(dim); + !(lo..=hi).contains(&point[dim]) + }) + } + + /// Interpolate at each of several points, sharing one grid across all of them. + /// + /// Named `batch_interpolate_f64`, not `batch_interpolate`; see + /// [`Self::interpolate_f64`]'s doc for why. + pub fn batch_interpolate_f64( + &self, + points: &[&[f64]], + ) -> Result, InterpolateError> + where + Dv::Elem: Num, + { + let mut out = vec![Dv::Elem::zero(); points.len()]; + self.batch_interpolate_f64_into(points, &mut out)?; + Ok(out) + } + + /// Unchecked batched [`Self::interpolate_f64_fast`], assuming every point is valid. + /// + /// Named `batch_interpolate_f64_fast_into`, not `batch_interpolate_fast_into`; + /// see [`Self::interpolate_f64`]'s doc for why. + /// + /// # Panics + /// Panics if `out.len() != points.len()`. + pub fn batch_interpolate_f64_fast_into(&self, points: &[&[f64]], out: &mut [Dv::Elem]) { let n = self.ndim(); for point in points { assert_eq!( point.len(), n, - "batch_interpolate_fast_into: point length mismatch" + "batch_interpolate_f64_fast_into: point length mismatch" ); } assert_eq!( out.len(), points.len(), - "batch_interpolate_fast_into: length mismatch" + "batch_interpolate_f64_fast_into: length mismatch" ); self.strategy .batch_interpolate_fast_into(&self.data, points, out) } - fn batch_interpolate_fast(&self, points: &[&[D::Elem]]) -> Vec + /// Batched [`Self::interpolate_f64_fast`], for use in hot loops where the caller + /// has already checked bounds. + /// + /// Named `batch_interpolate_f64_fast`, not `batch_interpolate_fast`; see + /// [`Self::interpolate_f64`]'s doc for why. + pub fn batch_interpolate_f64_fast(&self, points: &[&[f64]]) -> Vec where - D::Elem: Num + Copy, + Dv::Elem: Num, { + let mut out = vec![Dv::Elem::zero(); points.len()]; + self.batch_interpolate_f64_fast_into(points, &mut out); + out + } + + /// Return an interpolator with viewed data. + pub fn view(&self) -> InterpNDBase, ViewRepr<&Dv::Elem>, S> + where + S: for<'a> StrategyND, ViewRepr<&'a Dv::Elem>>, + Dg::Elem: Clone, + Dv::Elem: Clone, + { + InterpNDBase { + data: self.data.view(), + strategy: self.strategy.clone(), + extrapolate: self.extrapolate, + } + } + + /// Turn the interpolator into an owned variant, cloning the array elements if necessary. + pub fn into_owned(self) -> InterpNDBase, OwnedRepr, S> + where + S: StrategyND, OwnedRepr>, + Dg::Elem: Clone, + Dv::Elem: Clone, + { + InterpNDBase { + data: self.data.clone().into_owned(), + strategy: self.strategy.clone(), + extrapolate: self.extrapolate, + } + } +} + +impl Interpolator for InterpNDBase +where + D: Data + RawDataClone + Clone, + D::Elem: Num + PartialOrd + Euclid + NumCast + Copy + Debug, + S: StrategyND + Clone, +{ + #[inline] + fn ndim(&self) -> usize { + self.data.ndim() + } + + fn validate(&self) -> Result<(), ValidateError> { + InterpNDBase::validate(self) + } + + fn interpolate(&self, point: &[D::Elem]) -> Result { + let point: Vec = point + .iter() + .map(|&x| num_traits::cast(x).expect("point element must cast to f64")) + .collect(); + InterpNDBase::interpolate_f64(self, &point) + } + + fn set_extrapolate(&mut self, extrapolate: Extrapolate) -> Result<(), ValidateError> { + InterpNDBase::set_extrapolate(self, extrapolate) + } + + fn interpolate_fast(&self, point: &[D::Elem]) -> D::Elem { + let point: Vec = point + .iter() + .map(|&x| num_traits::cast(x).expect("point element must cast to f64")) + .collect(); + InterpNDBase::interpolate_f64_fast(self, &point) + } + + fn batch_interpolate_into( + &self, + points: &[&[D::Elem]], + out: &mut [D::Elem], + ) -> Result<(), InterpolateError> { + let points: Vec> = points + .iter() + .map(|point| { + point + .iter() + .map(|&x| num_traits::cast(x).expect("point element must cast to f64")) + .collect() + }) + .collect(); + let points: Vec<&[f64]> = points.iter().map(Vec::as_slice).collect(); + InterpNDBase::batch_interpolate_f64_into(self, &points, out) + } + + fn batch_interpolate_fast_into(&self, points: &[&[D::Elem]], out: &mut [D::Elem]) { + let points: Vec> = points + .iter() + .map(|point| { + point + .iter() + .map(|&x| num_traits::cast(x).expect("point element must cast to f64")) + .collect() + }) + .collect(); + let points: Vec<&[f64]> = points.iter().map(Vec::as_slice).collect(); + InterpNDBase::batch_interpolate_f64_fast_into(self, &points, out) + } + + fn batch_interpolate_fast(&self, points: &[&[D::Elem]]) -> Vec { let n = self.ndim(); for point in points { assert_eq!( @@ -567,12 +849,3 @@ where out } } - -extrapolate_impl!(InterpNDBase, StrategyND); -set_strategy_box_impl!(InterpNDBase, StrategyND); -set_strategy_enum_impl!( - InterpNDBase, - strategy::enums::StrategyNDEnum, - strategy::enums::StrategyNDEnum -); -any_interpolator_impl!(InterpND, StrategyND); diff --git a/src/interpolator/n/strategies.rs b/src/interpolator/n/strategies.rs index 4db09f2..d72ada7 100644 --- a/src/interpolator/n/strategies.rs +++ b/src/interpolator/n/strategies.rs @@ -1,16 +1,48 @@ use super::*; use strategy::*; -impl StrategyND for Linear +/// Casts a single grid coordinate to `f64`, the shared blend/query precision every +/// ND strategy computes in (see [`StrategyND`]'s docs). +fn to_f64(x: T) -> f64 { + num_traits::cast(x).expect("grid element must cast to f64") +} + +/// Casts a whole grid axis to `f64`, for a strategy that needs to do real arithmetic +/// across the axis (`LinearUniform`'s uniformity check, `CubicC2`'s corner cache, +/// `GridTransform`'s forward transform) rather than a single per-probe comparison +/// (which `locate_lower_index_cast`/`exact_index_cast`/`locate_step_index_cast` +/// handle without this allocation). +fn grid_to_f64(grid: ArrayView1) -> Array1 { + grid.iter().map(|&x| to_f64(x)).collect() +} + +/// Casts a blended `f64` result back down to `Dv::Elem`, for the checked (`Result` +/// returning) interpolation path. +/// +/// Known limitation: this goes through `NumCast`, whose float-to-integer conversion +/// truncates toward zero (`as`'s own semantics) rather than rounding to nearest. So +/// for an integer `Dv`, a blend that lands on e.g. `15.5` casts down to `15`, not a +/// rounded `16`, silently biasing every non-exact blend downward. Fixing this +/// properly needs a type-aware cast (round for integer `Dv`, exact passthrough for a +/// float `Dv`, since rounding *that* would wrongly destroy real fractional output), +/// deferred for now; revisit alongside a real `Tp` type if this prototype proceeds. +fn from_f64_checked(x: f64) -> Result { + num_traits::cast(x) + .ok_or_else(|| InterpolateError::Other("blended value doesn't fit in value type".into())) +} + +impl StrategyND for Linear where - D: Data + RawDataClone + Clone, - D::Elem: Float + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: NumCast + PartialOrd + Copy + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: NumCast + Copy + Debug + PartialEq, { fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { // Dimensionality let mut n = data.values.ndim(); @@ -22,13 +54,11 @@ where let mut values_view = data.values.view(); for dim in (0..n).rev() { // Skip empty grid dimensions (e.g. the 0-D multilinear case uses an empty grid). - // The original iter().position() returned None on empty grids without touching point[dim]; - // the binary search path would panic on first().unwrap(), so we guard it here. if grid[dim].is_empty() { continue; } - let lower = locate_lower_index(grid[dim].view(), &point[dim]); - let pos = exact_index(grid[dim].view(), lower, &point[dim]); + let lower = locate_lower_index_cast(grid[dim], point[dim]); + let pos = exact_index_cast(grid[dim], lower, point[dim]); if let Some(pos) = pos { point.remove(dim); grid.remove(dim); @@ -37,6 +67,7 @@ where } if values_view.len() == 1 { // Supplied point is coincident with a grid point, so just return the value + // directly: no blending, so no precision lost to the f64 round trip below. return Ok(values_view.first().copied().unwrap()); } // Simplified dimensionality @@ -49,35 +80,35 @@ where for dim in 0..n { // Extrapolation is checked previously in Interpolator::interpolate, // meaning by now, point is within grid bounds or extrapolation is enabled - let lower_idx = locate_lower_index(grid[dim].view(), &point[dim]); - let interp_diff = (point[dim] - grid[dim][lower_idx]) - / (grid[dim][lower_idx + 1] - grid[dim][lower_idx]); + let lower_idx = locate_lower_index_cast(grid[dim], point[dim]); + let g_lower = to_f64(grid[dim][lower_idx]); + let g_upper = to_f64(grid[dim][lower_idx + 1]); + let interp_diff = (point[dim] - g_lower) / (g_upper - g_lower); lower_idxs.push(lower_idx); interp_diffs.push(interp_diff); } - // Fill all 2^n corner values into a flat array indexed by bitmask. - // Bit (n-1-d) of the mask = 1 selects the upper index in dimension d. - // This layout supports an in-place butterfly reduction with no coordinate permutation tables. + // Fill all 2^n corner values into a flat array indexed by bitmask, blending in + // f64 (the shared precision, independent of Dv) so an integer-ish Dv (e.g. an + // image's u8 pixel values) doesn't lose precision to repeated integer blends. let size = 1usize << n; - let mut vals = vec![D::Elem::zero(); size]; + let mut vals = vec![0f64; size]; let mut idx = vec![0usize; n]; for (mask, val) in vals.iter_mut().enumerate() { for d in 0..n { idx[d] = lower_idxs[d] + ((mask >> (n - 1 - d)) & 1); } - *val = values_view[idx.as_slice()]; + *val = to_f64(values_view[idx.as_slice()]); } // Butterfly reduction: one pass per dimension. - // After pass d, vals[0..2^(n-d-1)] holds the result with dimensions 0..=d blended. for (d, diff) in interp_diffs.iter().enumerate() { let half = 1 << (n - 1 - d); for i in 0..half { - vals[i] = vals[i] * (D::Elem::one() - *diff) + vals[i + half] * *diff; + vals[i] = vals[i] * (1.0 - *diff) + vals[i + half] * *diff; } } - Ok(vals[0]) + from_f64_checked(vals[0]) } /// Returns `true`. @@ -86,52 +117,57 @@ where } } -impl StrategyND for LinearUniform +impl StrategyND for LinearUniform where - D: Data + RawDataClone + Clone, - D::Elem: Float + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: NumCast + PartialOrd + Copy + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: NumCast + Copy + Debug + PartialEq, { /// Ensures grid uniformity in all dimensions - fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { for (dim, grid) in data.grid.iter().enumerate() { - validate_uniform_grid_epsilon(grid.view(), dim, None)?; + let grid_f64 = grid_to_f64(grid.view()); + validate_uniform_grid_epsilon(grid_f64.view(), dim, None)?; } Ok(()) } fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { let n = data.values.ndim(); let mut lower_idxs = Vec::with_capacity(n); let mut interp_diffs = Vec::with_capacity(n); for (grid_dim, &point_dim) in data.grid.iter().zip(point.iter()) { - let step = grid_dim[1] - grid_dim[0]; - let lower_idx = - locate_lower_index_uniform(grid_dim[0], step, grid_dim.len(), point_dim); - let diff = (point_dim - grid_dim[lower_idx]) / step; + let g0 = to_f64(grid_dim[0]); + let g1 = to_f64(grid_dim[1]); + let step = g1 - g0; + let lower_idx = locate_lower_index_uniform(g0, step, grid_dim.len(), point_dim); + let g_lower = to_f64(grid_dim[lower_idx]); + let diff = (point_dim - g_lower) / step; lower_idxs.push(lower_idx); interp_diffs.push(diff); } // Same bitmask/butterfly reduction as Linear ND let size = 1usize << n; - let mut vals = vec![D::Elem::zero(); size]; + let mut vals = vec![0f64; size]; let mut idx = vec![0usize; n]; for (mask, val) in vals.iter_mut().enumerate() { for d in 0..n { idx[d] = lower_idxs[d] + ((mask >> (n - 1 - d)) & 1); } - *val = data.values.view()[idx.as_slice()]; + *val = to_f64(data.values.view()[idx.as_slice()]); } for (d, diff) in interp_diffs.iter().enumerate() { let half = 1 << (n - 1 - d); for i in 0..half { - vals[i] = vals[i] * (D::Elem::one() - *diff) + vals[i + half] * *diff; + vals[i] = vals[i] * (1.0 - *diff) + vals[i + half] * *diff; } } - Ok(vals[0]) + from_f64_checked(vals[0]) } /// Returns `true`. @@ -140,26 +176,28 @@ where } } -impl StrategyND for Nearest +impl StrategyND for Nearest where - D: Data + RawDataClone + Clone, - D::Elem: Sub + PartialOrd + Copy + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: NumCast + PartialOrd + Copy + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: Copy + Debug + PartialEq, { fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { let n = data.values.ndim(); // Nearest-neighbor on a rectilinear grid factorizes: select the nearest index - // independently per dimension, then do a single lookup. No corner extraction or - // dimensionality reduction needed; the distance comparison handles exact matches correctly. + // independently per dimension, then do a single lookup: a direct, uncast + // `Dv::Elem` read, so no precision is lost regardless of `Dv`. let mut idx = vec![0usize; n]; for dim in 0..n { - let lower_idx = locate_lower_index(data.grid[dim].view(), &point[dim]); - idx[dim] = if point[dim] - data.grid[dim][lower_idx] - < data.grid[dim][lower_idx + 1] - point[dim] - { + let lower_idx = locate_lower_index_cast(data.grid[dim].view(), point[dim]); + let lo = to_f64(data.grid[dim][lower_idx]); + let hi = to_f64(data.grid[dim][lower_idx + 1]); + idx[dim] = if point[dim] - lo < hi - point[dim] { lower_idx } else { lower_idx + 1 @@ -174,26 +212,29 @@ where } } -impl StrategyND for Step +impl StrategyND for Step where - D: Data + RawDataClone + Clone, - D::Elem: PartialOrd + Copy + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: NumCast + PartialOrd + Copy + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: Copy + Debug + PartialEq, { /// Ensures the number of provided step directions matches the dimensionality of the interpolator - fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { self.directions .validate_len(data.values.ndim(), "Step", "directions") } fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { let n = data.values.ndim(); let mut idx = vec![0usize; n]; for dim in 0..n { - idx[dim] = locate_step_index(self.directions[dim], data.grid[dim].view(), &point[dim]); + idx[dim] = + locate_step_index_cast(self.directions[dim], data.grid[dim].view(), point[dim]); } Ok(data.values.view()[idx.as_slice()]) } @@ -204,12 +245,14 @@ where } } -impl StrategyND for CubicC2 +impl StrategyND for CubicC2 where - D: Data + RawDataClone + Clone, - D::Elem: Float + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: NumCast + PartialOrd + Copy + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: NumCast + Copy + Debug + PartialEq, { - fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { self.boundary_conditions .validate_len(data.ndim(), "CubicC2", "boundary conditions")?; for dim in 0..data.ndim() { @@ -218,33 +261,39 @@ where Ok(()) } - /// Precomputes the full corner-derivative tensor via `compute_corner_cache`. - fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + /// Precomputes the full corner-derivative tensor via `compute_corner_cache`, + /// casting the grid and values to `f64` first (the tensor is a blend + /// intermediate, like `Linear`'s corner blend, so it's `f64`-typed regardless of + /// `Dg`/`Dv`). A one-time cost per `init`, not per query. + fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { if data.ndim() == 0 { return Ok(()); } - let data_view = data.view(); + let grids_f64: Vec> = data.grid.iter().map(|g| grid_to_f64(g.view())).collect(); + let grid_views: Vec> = grids_f64.iter().map(|g| g.view()).collect(); + let values_f64: ArrayD = data.values.map(|&x| to_f64(x)); self.cache = - compute_corner_cache(&data_view.grid, data_view.values, &self.boundary_conditions); + compute_corner_cache(&grid_views, values_f64.view(), &self.boundary_conditions); Ok(()) } fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { if data.ndim() == 0 { return data.values.first().copied().ok_or_else(|| { InterpolateError::Other("internal: 0-D interpolation data has no value".into()) }); } - let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); - Ok(evaluate_spline_corner_cached( - &grids, - self.cache.view(), - point, - )) + // Re-cast the grid to `f64` per query: `CubicC2`'s own struct (shared with + // `Interp1D`/`2D`/`3D`) isn't touched by this prototype, so there's no spare + // field here to cache it in the way `GridTransform::grid_cache` does. + let grids_f64: Vec> = data.grid.iter().map(|g| grid_to_f64(g.view())).collect(); + let grid_views: Vec> = grids_f64.iter().map(|g| g.view()).collect(); + let result = evaluate_spline_corner_cached(&grid_views, self.cache.view(), point); + from_f64_checked(result) } /// Returns `true`: the boundary cubic polynomials extend naturally. @@ -253,48 +302,48 @@ where } } -impl StrategyND for GridTransform +impl StrategyND for GridTransform where - D: Data + RawDataClone + Clone, - D::Elem: Float + Debug, - S: for<'a> StrategyND> + Clone + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: NumCast + PartialOrd + Copy + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, + // Both sides are view-repr'd here: `data.values.view()` produces a view + // regardless of what `Dv` itself is (owned or already a view), matching how the + // pre-split code bounded this the same way on its single shared type param. + S: for<'a> StrategyND, ViewRepr<&'a Dv::Elem>> + Clone + Debug, { - /// Checks the axis count and that every raw grid coordinate is in its axis's - /// configured transform's domain. - fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + /// Checks the axis count and that every raw grid coordinate (cast to `f64`) is + /// in its axis's configured transform's domain. + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { self.transforms .validate_len(data.ndim(), "GridTransform", "transforms")?; - let transformed_grid: Vec> = data + let transformed_grid: Vec> = data .grid .iter() .enumerate() - .map(|(dim, grid)| self.transform_axis(dim, grid.view())) + .map(|(dim, grid)| self.transform_axis(dim, grid_to_f64(grid.view()).view())) .collect::>()?; let values = self.transformed_values_view(data.values.view()); - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: transformed_grid.iter().map(|g| g.view()).collect(), values, }; self.inner.validate(&view) } - /// Transforms the grid into `grid_cache`, then initializes `inner` against a - /// transient view zipping `grid_cache` with `data.values`. - /// - /// A raw grid is always strictly increasing, but a decreasing transform (e.g. - /// `Reciprocal`) would otherwise leave that axis of `grid_cache` decreasing; - /// `Transform::is_increasing` flags that case so the transformed axis (and the - /// matching `values` axis) can be reversed back to ascending, matching every - /// downstream strategy's ascending-grid assumption. - fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + /// Transforms the (`f64`-cast) grid into `grid_cache`, then initializes `inner` + /// against a transient view zipping `grid_cache` with `data.values` (untouched; + /// `GridTransform` is grid-side only). + fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { self.grid_cache = data .grid .iter() .enumerate() - .map(|(dim, grid)| self.transform_axis(dim, grid.view())) + .map(|(dim, grid)| self.transform_axis(dim, grid_to_f64(grid.view()).view())) .collect::>()?; let values = self.transformed_values_view(data.values.view()); - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: self.grid_cache.iter().map(|g| g.view()).collect(), values, }; @@ -303,89 +352,64 @@ where fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { self.check_point_domain(point)?; - let transformed_point: Vec = point + let transformed_point: Vec = point .iter() .enumerate() .map(|(dim, &x)| self.transforms[dim].forward(x)) .collect(); let values = self.transformed_values_view(data.values.view()); - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: self.grid_cache.iter().map(|g| g.view()).collect(), values, }; self.inner.interpolate(&view, &transformed_point) } - /// Forward-transforms into `grid_cache`'s coordinate space, then delegates the - /// actual wrap to `inner.interpolate_wrapped` rather than wrapping here itself: - /// wrapping doesn't commute with a nonlinear transform, so it must happen in the - /// *final* (innermost) transformed space where the periodic strategy actually - /// lives, mirroring how `ValuesTransform::interpolate_wrapped` defers to `inner` - /// for the same reason (composing two `GridTransform`s and wrapping at the outer - /// layer's space, then forward-transforming again, uses the wrong period). For a - /// non-transform `inner`, `StrategyND::interpolate_wrapped`'s default wraps - /// directly against `grid_cache`, reproducing this layer's own wrap exactly. fn interpolate_wrapped( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result + data: &InterpDataNDBase, + point: &[f64], + ) -> Result where - D::Elem: Num + Euclid + Copy, + Dg::Elem: NumCast + Copy, { self.check_point_domain(point)?; - let transformed_point: Vec = point + let transformed_point: Vec = point .iter() .enumerate() .map(|(dim, &x)| self.transforms[dim].forward(x)) .collect(); let values = self.transformed_values_view(data.values.view()); - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: self.grid_cache.iter().map(|g| g.view()).collect(), values, }; self.inner.interpolate_wrapped(&view, &transformed_point) } - /// Skips the domain check `interpolate` does, and calls `inner.interpolate_fast` - /// rather than `inner.interpolate`, so "fast" propagates through nested - /// `GridTransform`/`ValuesTransform` layers instead of stopping at the outermost - /// one. - /// - /// An out-of-domain point is the caller's problem here, same as any other - /// unchecked `_fast` method: `forward`-ing it produces `NaN` (`Log`/`Sqrt`) or - /// `+-inf` (`Reciprocal`, only at `x == 0`). `NaN` poisons the grid search's - /// comparisons and panics; `+-inf` compares normally, so it's treated like an - /// ordinary out-of-bounds extrapolation query and silently produces `NaN` output - /// instead. Expected, not a bug: check [`Transform::in_domain`] yourself first if - /// you need a guarantee either way. - fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[D::Elem]) -> D::Elem { - let transformed_point: Vec = point + fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[f64]) -> Dv::Elem { + let transformed_point: Vec = point .iter() .enumerate() .map(|(dim, &x)| self.transforms[dim].forward(x)) .collect(); let values = self.transformed_values_view(data.values.view()); - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: self.grid_cache.iter().map(|g| g.view()).collect(), values, }; self.inner.interpolate_fast(&view, &transformed_point) } - /// Domain-checks every point in the batch before transforming, aggregating - /// every violation across the *whole batch* into one - /// [`InterpolateError::GridTransformDomain`] instead of erroring on the first - /// one, mirroring how `Extrapolate::Error` aggregates out-of-bounds points. fn batch_interpolate_into( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - out: &mut [D::Elem], + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], ) -> Result<(), InterpolateError> { if out.len() != points.len() { return Err(InterpolateError::OutputLength { @@ -394,7 +418,7 @@ where }); } self.check_batch_domain(points.iter().copied())?; - let transformed_points: Vec> = points + let transformed_points: Vec> = points .iter() .map(|point| { point @@ -404,10 +428,9 @@ where .collect() }) .collect(); - let transformed_refs: Vec<&[D::Elem]> = - transformed_points.iter().map(Vec::as_slice).collect(); + let transformed_refs: Vec<&[f64]> = transformed_points.iter().map(Vec::as_slice).collect(); let values = self.transformed_values_view(data.values.view()); - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: self.grid_cache.iter().map(|g| g.view()).collect(), values, }; @@ -415,19 +438,10 @@ where .batch_interpolate_into(&view, &transformed_refs, out) } - /// Checks this layer's own domain per point (not bailing on the first violation), - /// then recurses into `inner` with only the *outer-valid* points' forward - /// transforms, remapping `inner`'s failure indices back to their true batch - /// position: an outer-invalid point can't be forward-transformed meaningfully, but - /// its presence must not hide `inner`'s own violations for the *other* points, so - /// a `GridTransform` nested inside another one still gets every one of its - /// violations aggregated, not just the outer layer's. Exposed here so generic - /// callers (e.g. `Extrapolate::Wrap`'s batch dispatch) can pre-scan a batch - /// without knowing the concrete strategy type. - fn check_batch_domain(&self, points: &[&[D::Elem]]) -> Result<(), InterpolateError> { + fn check_batch_domain(&self, points: &[&[f64]]) -> Result<(), InterpolateError> { let mut failures = Vec::new(); let mut valid_indices = Vec::new(); - let mut transformed_points: Vec> = Vec::new(); + let mut transformed_points: Vec> = Vec::new(); for (index, &point) in points.iter().enumerate() { let point_failures = self.point_domain_failures(index, point); if point_failures.is_empty() { @@ -443,8 +457,7 @@ where failures.extend(point_failures); } } - let transformed_refs: Vec<&[D::Elem]> = - transformed_points.iter().map(Vec::as_slice).collect(); + let transformed_refs: Vec<&[f64]> = transformed_points.iter().map(Vec::as_slice).collect(); match self.inner.check_batch_domain(&transformed_refs) { Ok(()) => {} Err(InterpolateError::GridTransformDomain(inner_failures)) => { @@ -468,27 +481,36 @@ where } } -impl StrategyND for ValuesTransform +impl StrategyND for ValuesTransform where - D: Data + RawDataClone + Clone, - D::Elem: Float + Debug, - S: for<'a> StrategyND> + Clone + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: NumCast + Copy + Debug + PartialEq, + // Both sides are view-repr'd here: `data.grid[i].view()` produces a view + // regardless of what `Dg` itself is, matching how the pre-split code bounded + // this the same way on its single shared type param. + S: for<'a> StrategyND, ViewRepr<&'a f64>> + Clone + Debug, { - /// Checks that every data value is in the configured transform's domain. - fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { - let transformed_values = self.transform_values(data.values.view())?; - let view = InterpDataNDView { + /// Checks that every data value (cast to `f64`) is in the configured transform's + /// domain. + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + let values_f64: ArrayD = data.values.map(|&x| to_f64(x)); + let transformed_values = self.transform_values(values_f64.view())?; + let view = InterpDataNDBase { grid: data.grid.iter().map(|g| g.view()).collect(), values: transformed_values.view(), }; self.inner.validate(&view) } - /// Transforms `data.values` into `values_cache`, then initializes `inner` - /// against a transient view zipping `data.grid` (untouched) with `values_cache`. - fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { - self.values_cache = self.transform_values(data.values.view())?; - let view = InterpDataNDView { + /// Transforms `data.values` (cast to `f64`) into `values_cache`, then + /// initializes `inner` against a transient view zipping `data.grid` (untouched; + /// `ValuesTransform` is value-side only) with `values_cache`. + fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + let values_f64: ArrayD = data.values.map(|&x| to_f64(x)); + self.values_cache = self.transform_values(values_f64.view())?; + let view = InterpDataNDBase { grid: data.grid.iter().map(|g| g.view()).collect(), values: self.values_cache.view(), }; @@ -497,55 +519,57 @@ where fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { - let view = InterpDataNDView { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { + let view = InterpDataNDBase { grid: data.grid.iter().map(|g| g.view()).collect(), values: self.values_cache.view(), }; let result = self.inner.interpolate(&view, point)?; - Ok(self.transform.inverse(result)) + from_f64_checked(self.transform.inverse(result)) } /// Hands `point` unmodified to `inner.interpolate_wrapped`, so a nested /// `GridTransform` handles the actual raw-space wrap; must not wrap here itself. fn interpolate_wrapped( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result + data: &InterpDataNDBase, + point: &[f64], + ) -> Result where - D::Elem: Num + Euclid + Copy, + Dg::Elem: NumCast + Copy, { - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: data.grid.iter().map(|g| g.view()).collect(), values: self.values_cache.view(), }; let result = self.inner.interpolate_wrapped(&view, point)?; - Ok(self.transform.inverse(result)) + from_f64_checked(self.transform.inverse(result)) } /// Calls `inner.interpolate_fast` rather than `inner.interpolate`, so "fast" /// propagates through nested `GridTransform`/`ValuesTransform` layers instead of /// stopping at the outermost one. - fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[D::Elem]) -> D::Elem { - let view = InterpDataNDView { + /// + /// Truncates rather than rounds for an integer `Dv`, same as `from_f64_checked`. + fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[f64]) -> Dv::Elem { + let view = InterpDataNDBase { grid: data.grid.iter().map(|g| g.view()).collect(), values: self.values_cache.view(), }; let result = self.inner.interpolate_fast(&view, point); - self.transform.inverse(result) + num_traits::cast(self.transform.inverse(result)) + .expect("inverse-transformed value doesn't fit in value type") } - /// Delegates the whole batch to `inner` via a transient view over `values_cache`, - /// so a nested `GridTransform`'s batch-domain aggregation isn't lost to the - /// point-by-point default; inverse-transforms every output afterward. + /// Delegates the whole batch to `inner` via a transient `f64` view over + /// `values_cache`, then casts + inverse-transforms every output afterward. fn batch_interpolate_into( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - out: &mut [D::Elem], + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], ) -> Result<(), InterpolateError> { if out.len() != points.len() { return Err(InterpolateError::OutputLength { @@ -553,20 +577,22 @@ where found: out.len(), }); } - let view = InterpDataNDView { + let view = InterpDataNDBase { grid: data.grid.iter().map(|g| g.view()).collect(), values: self.values_cache.view(), }; - self.inner.batch_interpolate_into(&view, points, out)?; - for o in out.iter_mut() { - *o = self.transform.inverse(*o); + let mut scratch = vec![0f64; points.len()]; + self.inner + .batch_interpolate_into(&view, points, &mut scratch)?; + for (o, v) in out.iter_mut().zip(scratch) { + *o = from_f64_checked(self.transform.inverse(v))?; } Ok(()) } /// Forwards to `inner`, so a nested `GridTransform`'s domain is still checked /// through a wrapping `ValuesTransform`. - fn check_batch_domain(&self, points: &[&[D::Elem]]) -> Result<(), InterpolateError> { + fn check_batch_domain(&self, points: &[&[f64]]) -> Result<(), InterpolateError> { self.inner.check_batch_domain(points) } diff --git a/src/interpolator/n/tests.rs b/src/interpolator/n/tests.rs index e61ec4f..94b5877 100644 --- a/src/interpolator/n/tests.rs +++ b/src/interpolator/n/tests.rs @@ -954,3 +954,29 @@ fn test_dyn_interpolator_heterogeneous_storage() { .downcast_ref::>() .is_none()); } + +/// Tg/Tv split prototype (issue #57): an integer grid (`Dg::Elem = i32`, no `Float` +/// bound at all) paired with `u8` values through `Linear`. Reached through +/// `InterpNDBase`'s inherent `interpolate_f64` directly, since `Dg != Dv` can't go +/// through the `Interpolator` trait (see `InterpNDBase`'s own docs). +#[test] +fn test_tg_tv_split_integer_grid_u8_values() { + let interp: InterpNDBase, OwnedRepr, strategy::Linear> = InterpNDBase::new( + vec![array![0i32, 10, 20]], + array![10u8, 20, 30].into_dyn(), + strategy::Linear, + Extrapolate::Error, + ) + .unwrap(); + + // Blended point: halfway between grid[0]=0 (value 10) and grid[1]=10 (value 20), + // blended in f64 precision and cast down to u8 once at the end. + assert_eq!(interp.interpolate_f64(&[5.0]).unwrap(), 15u8); + + // Exact grid hit: a direct, uncast `u8` read, no precision lost to a blend that + // never happens. + assert_eq!(interp.interpolate_f64(&[10.0]).unwrap(), 20u8); + + // Out-of-bounds query still goes through the usual `Extrapolate::Error` path. + assert!(interp.interpolate_f64(&[25.0]).is_err()); +} diff --git a/src/lib.rs b/src/lib.rs index 90e0101..9c2d882 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -73,7 +73,7 @@ pub(crate) use ndarray::prelude::*; pub(crate) use ndarray::{Data, IntoDimension, Ix, OwnedRepr, RawDataClone, Slice, ViewRepr}; pub use num_traits; -pub(crate) use num_traits::{clamp, Euclid, Float, Num, One, Zero}; +pub(crate) use num_traits::{clamp, Euclid, Float, Num, NumCast, One, Zero}; pub(crate) use core::ops::Sub; diff --git a/src/strategy/enums/mod.rs b/src/strategy/enums/mod.rs index 91cabdc..f831689 100644 --- a/src/strategy/enums/mod.rs +++ b/src/strategy/enums/mod.rs @@ -411,7 +411,7 @@ mod tests { #[test] fn test_nd() { - let mut interp: InterpND<_, strategy::enums::StrategyNDEnum> = InterpND::new( + let mut interp: InterpND<_, strategy::enums::StrategyNDEnum> = InterpND::new( vec![ array![0.05, 0.10, 0.15], array![0.10, 0.20, 0.30], diff --git a/src/strategy/enums/n.rs b/src/strategy/enums/n.rs index c124a21..de0a00e 100644 --- a/src/strategy/enums/n.rs +++ b/src/strategy/enums/n.rs @@ -1,21 +1,384 @@ use super::*; -strategy_enum_impl!( - StrategyNDEnum, - StrategyND, - InterpDataNDBase, - &[D::Elem], - &[&[D::Elem]], - [ - (Nearest, strategy::Nearest), - (Step, strategy::Step), - (Linear, strategy::Linear), - (LinearUniform, strategy::LinearUniform), - (CubicC2, strategy::CubicC2), - (GridTransform, strategy::GridTransform>>), - (ValuesTransform, strategy::ValuesTransform>>), - ] -); +/// See [enums module](super) documentation. +/// +/// Hand-written (not generated by the `strategy_enum_impl!` macro used in this +/// module's parent, unlike `Strategy1D`/`2D`/`3DEnum`) and **not** generic over a +/// grid/value type at all: every built-in ND strategy either does no arithmetic at +/// all (`Nearest`, `Step`, genuinely `Dg`/`Dv`-independent) or works in the shared +/// `f64` blend precision (`Linear`, `LinearUniform`, `CubicC2`, `GridTransform`, +/// `ValuesTransform`, see [`StrategyND`]'s docs), so nothing stored here is actually +/// `Tg`- or `Tv`-typed. +/// `impl StrategyND for StrategyNDEnum` below is a blanket impl over +/// every `Dg`/`Dv` pair meeting the bound, rather than one type per instantiation. +#[allow(missing_docs)] +#[derive(Debug, Clone, PartialEq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +#[cfg_attr(feature = "serde", serde(untagged))] +#[non_exhaustive] +pub enum StrategyNDEnum { + Nearest(strategy::Nearest), + Step(strategy::Step), + Linear(strategy::Linear), + LinearUniform(strategy::LinearUniform), + CubicC2(strategy::CubicC2), + GridTransform(strategy::GridTransform>), + ValuesTransform(strategy::ValuesTransform>), +} + +impl From for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::Nearest) -> Self { + Self::Nearest(strategy) + } +} +impl From for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::Step) -> Self { + Self::Step(strategy) + } +} +impl From for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::Linear) -> Self { + Self::Linear(strategy) + } +} +impl From for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::LinearUniform) -> Self { + Self::LinearUniform(strategy) + } +} +impl From> for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::CubicC2) -> Self { + Self::CubicC2(strategy) + } +} +impl From>> for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::GridTransform>) -> Self { + Self::GridTransform(strategy) + } +} +impl From>> for StrategyNDEnum { + #[inline] + fn from(strategy: strategy::ValuesTransform>) -> Self { + Self::ValuesTransform(strategy) + } +} + +impl StrategyND for StrategyNDEnum +where + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug + NumCast + PartialOrd + Copy, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug + NumCast + Copy, +{ + #[inline] + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + match self { + Self::Nearest(s) => StrategyND::::validate(s, data), + Self::Step(s) => StrategyND::::validate(s, data), + Self::Linear(s) => StrategyND::::validate(s, data), + Self::LinearUniform(s) => StrategyND::::validate(s, data), + Self::CubicC2(s) => StrategyND::::validate(s, data), + Self::GridTransform(s) => StrategyND::::validate(s, data), + Self::ValuesTransform(s) => StrategyND::::validate(s, data), + } + } + + #[inline] + fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + match self { + Self::Nearest(s) => StrategyND::::init(s, data), + Self::Step(s) => StrategyND::::init(s, data), + Self::Linear(s) => StrategyND::::init(s, data), + Self::LinearUniform(s) => StrategyND::::init(s, data), + Self::CubicC2(s) => StrategyND::::init(s, data), + Self::GridTransform(s) => StrategyND::::init(s, data), + Self::ValuesTransform(s) => StrategyND::::init(s, data), + } + } + + #[inline] + fn interpolate( + &self, + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { + match self { + Self::Nearest(s) => StrategyND::::interpolate(s, data, point), + Self::Step(s) => StrategyND::::interpolate(s, data, point), + Self::Linear(s) => StrategyND::::interpolate(s, data, point), + Self::LinearUniform(s) => StrategyND::::interpolate(s, data, point), + Self::CubicC2(s) => StrategyND::::interpolate(s, data, point), + Self::GridTransform(s) => StrategyND::::interpolate(s, data, point), + Self::ValuesTransform(s) => StrategyND::::interpolate(s, data, point), + } + } + + #[inline] + fn interpolate_wrapped( + &self, + data: &InterpDataNDBase, + point: &[f64], + ) -> Result + where + Dg::Elem: NumCast + Copy, + { + match self { + Self::Nearest(s) => StrategyND::::interpolate_wrapped(s, data, point), + Self::Step(s) => StrategyND::::interpolate_wrapped(s, data, point), + Self::Linear(s) => StrategyND::::interpolate_wrapped(s, data, point), + Self::LinearUniform(s) => StrategyND::::interpolate_wrapped(s, data, point), + Self::CubicC2(s) => StrategyND::::interpolate_wrapped(s, data, point), + Self::GridTransform(s) => StrategyND::::interpolate_wrapped(s, data, point), + Self::ValuesTransform(s) => StrategyND::::interpolate_wrapped(s, data, point), + } + } + + #[inline] + fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[f64]) -> Dv::Elem { + match self { + Self::Nearest(s) => StrategyND::::interpolate_fast(s, data, point), + Self::Step(s) => StrategyND::::interpolate_fast(s, data, point), + Self::Linear(s) => StrategyND::::interpolate_fast(s, data, point), + Self::LinearUniform(s) => StrategyND::::interpolate_fast(s, data, point), + Self::CubicC2(s) => StrategyND::::interpolate_fast(s, data, point), + Self::GridTransform(s) => StrategyND::::interpolate_fast(s, data, point), + Self::ValuesTransform(s) => StrategyND::::interpolate_fast(s, data, point), + } + } + + #[inline] + fn batch_interpolate_into( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], + ) -> Result<(), InterpolateError> { + match self { + Self::Nearest(s) => StrategyND::::batch_interpolate_into(s, data, points, out), + Self::Step(s) => StrategyND::::batch_interpolate_into(s, data, points, out), + Self::Linear(s) => StrategyND::::batch_interpolate_into(s, data, points, out), + Self::LinearUniform(s) => { + StrategyND::::batch_interpolate_into(s, data, points, out) + } + Self::CubicC2(s) => StrategyND::::batch_interpolate_into(s, data, points, out), + Self::GridTransform(s) => { + StrategyND::::batch_interpolate_into(s, data, points, out) + } + Self::ValuesTransform(s) => { + StrategyND::::batch_interpolate_into(s, data, points, out) + } + } + } + + #[inline] + fn batch_interpolate_fast_into( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], + ) { + match self { + Self::Nearest(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + Self::Step(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + Self::Linear(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + Self::LinearUniform(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + Self::CubicC2(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + Self::GridTransform(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + Self::ValuesTransform(s) => { + StrategyND::::batch_interpolate_fast_into(s, data, points, out) + } + } + } + + #[inline] + fn batch_interpolate( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Result, InterpolateError> + where + Dv::Elem: Num, + { + match self { + Self::Nearest(s) => StrategyND::::batch_interpolate(s, data, points), + Self::Step(s) => StrategyND::::batch_interpolate(s, data, points), + Self::Linear(s) => StrategyND::::batch_interpolate(s, data, points), + Self::LinearUniform(s) => StrategyND::::batch_interpolate(s, data, points), + Self::CubicC2(s) => StrategyND::::batch_interpolate(s, data, points), + Self::GridTransform(s) => StrategyND::::batch_interpolate(s, data, points), + Self::ValuesTransform(s) => StrategyND::::batch_interpolate(s, data, points), + } + } + + #[inline] + fn batch_interpolate_fast( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Vec + where + Dv::Elem: Num + Copy, + { + match self { + Self::Nearest(s) => StrategyND::::batch_interpolate_fast(s, data, points), + Self::Step(s) => StrategyND::::batch_interpolate_fast(s, data, points), + Self::Linear(s) => StrategyND::::batch_interpolate_fast(s, data, points), + Self::LinearUniform(s) => StrategyND::::batch_interpolate_fast(s, data, points), + Self::CubicC2(s) => StrategyND::::batch_interpolate_fast(s, data, points), + Self::GridTransform(s) => StrategyND::::batch_interpolate_fast(s, data, points), + Self::ValuesTransform(s) => { + StrategyND::::batch_interpolate_fast(s, data, points) + } + } + } + + #[inline] + fn allow_extrapolate(&self) -> bool { + match self { + Self::Nearest(s) => StrategyND::::allow_extrapolate(s), + Self::Step(s) => StrategyND::::allow_extrapolate(s), + Self::Linear(s) => StrategyND::::allow_extrapolate(s), + Self::LinearUniform(s) => StrategyND::::allow_extrapolate(s), + Self::CubicC2(s) => StrategyND::::allow_extrapolate(s), + Self::GridTransform(s) => StrategyND::::allow_extrapolate(s), + Self::ValuesTransform(s) => StrategyND::::allow_extrapolate(s), + } + } + + #[inline] + fn check_batch_domain(&self, points: &[&[f64]]) -> Result<(), InterpolateError> { + match self { + Self::Nearest(s) => StrategyND::::check_batch_domain(s, points), + Self::Step(s) => StrategyND::::check_batch_domain(s, points), + Self::Linear(s) => StrategyND::::check_batch_domain(s, points), + Self::LinearUniform(s) => StrategyND::::check_batch_domain(s, points), + Self::CubicC2(s) => StrategyND::::check_batch_domain(s, points), + Self::GridTransform(s) => StrategyND::::check_batch_domain(s, points), + Self::ValuesTransform(s) => StrategyND::::check_batch_domain(s, points), + } + } +} + +/// See [enums module](super) documentation. `Box` here wraps the concrete enum (not +/// `dyn Trait`), breaking the infinite size a wrapper strategy's `inner: Self` field +/// would otherwise have if it named the enum directly (e.g. `GridTransform`'s +/// `inner`); serde's blanket `Box` impl covers it, so unlike `Box` this +/// is fully serde-compatible. +impl StrategyND for Box +where + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug + NumCast + PartialOrd + Copy, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug + NumCast + Copy, +{ + #[inline] + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + (**self).validate(data) + } + + #[inline] + fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + (**self).init(data) + } + + #[inline] + fn interpolate( + &self, + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { + (**self).interpolate(data, point) + } + + #[inline] + fn interpolate_wrapped( + &self, + data: &InterpDataNDBase, + point: &[f64], + ) -> Result + where + Dg::Elem: NumCast + Copy, + { + (**self).interpolate_wrapped(data, point) + } + + #[inline] + fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[f64]) -> Dv::Elem { + (**self).interpolate_fast(data, point) + } + + #[inline] + fn batch_interpolate_into( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], + ) -> Result<(), InterpolateError> { + (**self).batch_interpolate_into(data, points, out) + } + + #[inline] + fn batch_interpolate_fast_into( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], + ) { + (**self).batch_interpolate_fast_into(data, points, out) + } + + #[inline] + fn batch_interpolate( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Result, InterpolateError> + where + Dv::Elem: Num, + { + (**self).batch_interpolate(data, points) + } + + #[inline] + fn batch_interpolate_fast( + &self, + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Vec + where + Dv::Elem: Num + Copy, + { + (**self).batch_interpolate_fast(data, points) + } + + #[inline] + fn allow_extrapolate(&self) -> bool { + StrategyND::::allow_extrapolate(&**self) + } + + #[inline] + fn check_batch_domain(&self, points: &[&[f64]]) -> Result<(), InterpolateError> { + StrategyND::::check_batch_domain(&**self, points) + } +} #[cfg(test)] mod tests { @@ -26,15 +389,15 @@ mod tests { #[cfg(feature = "serde")] fn test_serde() { assert_eq!( - serde_json::to_string(&StrategyNDEnum::::from(Linear)).unwrap(), + serde_json::to_string(&StrategyNDEnum::from(Linear)).unwrap(), serde_json::to_string(&Linear).unwrap(), ); assert_eq!( - serde_json::to_string(&StrategyNDEnum::::from(Nearest)).unwrap(), + serde_json::to_string(&StrategyNDEnum::from(Nearest)).unwrap(), serde_json::to_string(&Nearest).unwrap(), ); assert_eq!( - serde_json::to_string(&StrategyNDEnum::::from(Step::from( + serde_json::to_string(&StrategyNDEnum::from(Step::from( strategy::step::StepDirection::Lower ))) .unwrap(), diff --git a/src/strategy/traits.rs b/src/strategy/traits.rs index 041e471..e8f3f03 100644 --- a/src/strategy/traits.rs +++ b/src/strategy/traits.rs @@ -4,7 +4,9 @@ use super::*; /// Generates a fixed-dimensionality strategy trait (`Strategy1D`/`2D`/`3D`) plus its /// `impl … for Box` forwarding impl. `StrategyND` takes points as `&[D::Elem]` -/// instead of `&[D::Elem; N]`, so it can't share this shape and stays hand-written below. +/// instead of `&[D::Elem; N]`, and is split into separate grid/value type parameters +/// (`Dg`/`Dv`, with `point: &[f64]`, see [`StrategyND`]'s docs), so it can't share +/// this shape and stays hand-written below. macro_rules! fixed_strategy_trait { ($Trait:ident, $InterpData:ident, $N:literal, $doc:literal) => { #[doc = $doc] @@ -323,16 +325,28 @@ fixed_strategy_trait!( ); /// N-D interpolation strategy. -pub trait StrategyND: Debug + DynClone +/// +/// Split into a grid type `Dg` and a value type `Dv` (ninterp's Tg/Tv split, +/// prototyped here first, see issue #57). The query `point` is always `&[f64]`: +/// a query point falls *between* grid points, so it can't be `Dg`-typed once `Dg` +/// isn't assumed `Float` (an integer or date-like grid can't represent a fractional +/// position), and `f64` serves as the shared continuous precision both `Dg` and `Dv` +/// cast into/out of for any strategy that needs real arithmetic (fractional +/// position, blending). A strategy that does no arithmetic (`Nearest`, `Step`) still +/// takes `point: &[f64]` for a uniform trait shape, casting `Dg` into `f64` only to +/// compare. +pub trait StrategyND: Debug + DynClone where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, { /// Validate strategy state against interpolation data. Pure check, no mutation. /// /// Default no-op. Override for invariant checks that don't require precomputed /// state (grid uniformity, direction-count matching, etc). - fn validate(&self, _data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn validate(&self, _data: &InterpDataNDBase) -> Result<(), ValidateError> { Ok(()) } @@ -341,7 +355,7 @@ where /// Default no-op. Override only when the strategy caches something derived from /// `data` (e.g. precomputed spline coefficients). Unlike [`StrategyND::validate`], /// this may do real, non-trivial calculation. - fn init(&mut self, _data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn init(&mut self, _data: &InterpDataNDBase) -> Result<(), ValidateError> { Ok(()) } @@ -352,34 +366,36 @@ where /// panics on non-contiguous storage (possible with `Interp*View`). See /// [`crate::strategy::utils`] for ready-made per-axis search helpers (bracket search, /// exact-match short-circuit, step-direction lookup, uniform-grid fast path) built from - /// the same primitives the built-in strategies use. + /// the same primitives the built-in strategies use, including cast-aware variants for + /// when `Dg::Elem` isn't `f64` itself. fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result; + data: &InterpDataNDBase, + point: &[f64], + ) -> Result; /// Resolves an out-of-bounds point under [`Extrapolate::Wrap`](`crate::interpolator::Extrapolate::Wrap`), then interpolates. /// - /// Default wraps in `data.grid`'s own (raw) coordinate space. Override only if - /// this strategy's working coordinate space differs from `data.grid`'s. + /// Default wraps in `data.grid`'s own (raw) coordinate space, casting each grid + /// bound to `f64` to match `point`. Override only if this strategy's working + /// coordinate space differs from `data.grid`'s. fn interpolate_wrapped( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result + data: &InterpDataNDBase, + point: &[f64], + ) -> Result where - D::Elem: Num + Euclid + Copy, + Dg::Elem: NumCast + Copy, { - let wrapped: Vec = point + let wrapped: Vec = point .iter() .enumerate() .map(|(dim, &pt)| { - wrap( - pt, - *data.grid[dim].first().unwrap(), - *data.grid[dim].last().unwrap(), - ) + let first: f64 = num_traits::cast(*data.grid[dim].first().unwrap()) + .expect("grid element must cast to f64"); + let last: f64 = num_traits::cast(*data.grid[dim].last().unwrap()) + .expect("grid element must cast to f64"); + wrap(pt, first, last) }) .collect(); self.interpolate(data, &wrapped) @@ -392,7 +408,7 @@ where /// strategy's checked path does real internal fallible work beyond producing /// the final `Ok(...)`; otherwise the default already compiles to the same thing. #[inline] - fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[D::Elem]) -> D::Elem { + fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[f64]) -> Dv::Elem { self.interpolate(data, point) .expect("interpolate_fast: invalid point or data") } @@ -406,9 +422,9 @@ where /// no strategy shipped in this crate does that today. fn batch_interpolate_into( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - out: &mut [D::Elem], + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], ) -> Result<(), InterpolateError> { if out.len() != points.len() { return Err(InterpolateError::OutputLength { @@ -429,9 +445,9 @@ where /// [`StrategyND::interpolate_fast`]. fn batch_interpolate_fast_into( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - out: &mut [D::Elem], + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], ) { assert_eq!( out.len(), @@ -451,15 +467,15 @@ where /// for a locate sweep instead of one binary search per point). fn batch_interpolate( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - ) -> Result, InterpolateError> + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Result, InterpolateError> where - D::Elem: Num, + Dv::Elem: Num, { let mut out = Vec::with_capacity(points.len()); for _ in 0..points.len() { - out.push(D::Elem::zero()); + out.push(Dv::Elem::zero()); } self.batch_interpolate_into(data, points, &mut out)?; Ok(out) @@ -474,15 +490,15 @@ where /// with no additional override needed. fn batch_interpolate_fast( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - ) -> Vec + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Vec where - D::Elem: Num + Copy, + Dv::Elem: Num + Copy, { let mut out = Vec::with_capacity(points.len()); for _ in 0..points.len() { - out.push(D::Elem::zero()); + out.push(Dv::Elem::zero()); } self.batch_interpolate_fast_into(data, points, &mut out); out @@ -492,11 +508,11 @@ where /// (distinct from `data.grid`'s bounds), aggregating every violation across the /// whole batch instead of just the first. /// - /// Default no-op: most strategies accept any `D::Elem`. Override only if `self` + /// Default no-op: most strategies accept any point. Override only if `self` /// (or a wrapped inner strategy) restricts the domain further, e.g. /// [`GridTransform`]'s configured [`Transform`]. Called as a pre-scan before /// doing any actual interpolation work, so it must stay cheap. - fn check_batch_domain(&self, _points: &[&[D::Elem]]) -> Result<(), InterpolateError> { + fn check_batch_domain(&self, _points: &[&[f64]]) -> Result<(), InterpolateError> { Ok(()) } @@ -504,29 +520,31 @@ where fn allow_extrapolate(&self) -> bool; } -clone_trait_object!( StrategyND); +clone_trait_object!( StrategyND); -impl StrategyND for Box> +impl StrategyND for Box> where - D: Data + RawDataClone + Clone, - D::Elem: PartialEq + Debug, + Dg: Data + RawDataClone + Clone, + Dg::Elem: PartialEq + Debug, + Dv: Data + RawDataClone + Clone, + Dv::Elem: PartialEq + Debug, { #[inline] - fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn validate(&self, data: &InterpDataNDBase) -> Result<(), ValidateError> { (**self).validate(data) } #[inline] - fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { + fn init(&mut self, data: &InterpDataNDBase) -> Result<(), ValidateError> { (**self).init(data) } #[inline] fn interpolate( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result { + data: &InterpDataNDBase, + point: &[f64], + ) -> Result { (**self).interpolate(data, point) } @@ -536,33 +554,33 @@ where } #[inline] - fn check_batch_domain(&self, points: &[&[D::Elem]]) -> Result<(), InterpolateError> { + fn check_batch_domain(&self, points: &[&[f64]]) -> Result<(), InterpolateError> { (**self).check_batch_domain(points) } #[inline] fn interpolate_wrapped( &self, - data: &InterpDataNDBase, - point: &[D::Elem], - ) -> Result + data: &InterpDataNDBase, + point: &[f64], + ) -> Result where - D::Elem: Num + Euclid + Copy, + Dg::Elem: NumCast + Copy, { (**self).interpolate_wrapped(data, point) } #[inline] - fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[D::Elem]) -> D::Elem { + fn interpolate_fast(&self, data: &InterpDataNDBase, point: &[f64]) -> Dv::Elem { (**self).interpolate_fast(data, point) } #[inline] fn batch_interpolate_into( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - out: &mut [D::Elem], + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], ) -> Result<(), InterpolateError> { (**self).batch_interpolate_into(data, points, out) } @@ -570,9 +588,9 @@ where #[inline] fn batch_interpolate_fast_into( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - out: &mut [D::Elem], + data: &InterpDataNDBase, + points: &[&[f64]], + out: &mut [Dv::Elem], ) { (**self).batch_interpolate_fast_into(data, points, out) } @@ -580,11 +598,11 @@ where #[inline] fn batch_interpolate( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - ) -> Result, InterpolateError> + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Result, InterpolateError> where - D::Elem: Num, + Dv::Elem: Num, { (**self).batch_interpolate(data, points) } @@ -592,11 +610,11 @@ where #[inline] fn batch_interpolate_fast( &self, - data: &InterpDataNDBase, - points: &[&[D::Elem]], - ) -> Vec + data: &InterpDataNDBase, + points: &[&[f64]], + ) -> Vec where - D::Elem: Num + Copy, + Dv::Elem: Num + Copy, { (**self).batch_interpolate_fast(data, points) } diff --git a/src/strategy/utils.rs b/src/strategy/utils.rs index eb93e0f..e586046 100644 --- a/src/strategy/utils.rs +++ b/src/strategy/utils.rs @@ -42,6 +42,66 @@ pub fn locate_lower_index(grid: ArrayView1, point: &T) -> usiz } } +/// [`locate_lower_index`], but for a `grid` whose element type differs from `point`'s: +/// casts each probed grid element to `f64` via [`NumCast`] for the comparison instead +/// of requiring both to share a type. Used where the query point has been normalized +/// to `f64` (ninterp's Tg/Tv split; see the `n` module) ahead of a grid whose own +/// element type may not even be [`Float`]. +/// +/// # Panics +/// Panics if any grid element fails to cast to `f64` (not expected for any numeric `T`). +pub fn locate_lower_index_cast( + grid: ArrayView1, + point: f64, +) -> usize { + let cast = |x: T| -> f64 { num_traits::cast(x).expect("grid element must cast to f64") }; + let first = cast(*grid.first().unwrap()); + let last = cast(*grid.last().unwrap()); + if point < first { + return 0; + } + if point >= last { + return grid.len() - 2; + } + + let mut low = 0; + let mut high = grid.len() - 1; + + while low < high { + let mid = low + (high - low) / 2; + + if cast(grid[mid]) >= point { + high = mid; + } else { + low = mid + 1; + } + } + + if low > 0 && cast(grid[low]) >= point { + low - 1 + } else { + low + } +} + +/// [`exact_index`], but for a `grid` whose element type differs from `point`'s: casts +/// each probed grid element to `f64` via [`NumCast`] for the comparison. See +/// [`locate_lower_index_cast`]. +pub fn exact_index_cast( + grid: ArrayView1, + lower: usize, + point: f64, +) -> Option { + let cast = |x: T| -> f64 { num_traits::cast(x).expect("grid element must cast to f64") }; + if cast(grid[lower]) == point { + Some(lower) + } else if cast(grid[lower + 1]) == point { + Some(lower + 1) + } else { + None + } +} + /// Per-axis locate for linear-family strategies: either an exact grid hit, /// or an interior interpolation position. pub enum AxisLocation { @@ -105,6 +165,36 @@ pub fn locate_step_index( } } +/// [`locate_step_index`], but for a `grid` whose element type differs from `point`'s: +/// casts each probed grid element to `f64` via [`NumCast`] for the comparison. See +/// [`locate_lower_index_cast`]. +pub fn locate_step_index_cast( + dir: StepDirection, + grid: ArrayView1, + point: f64, +) -> usize { + let cast = |x: T| -> f64 { num_traits::cast(x).expect("grid element must cast to f64") }; + match dir { + StepDirection::Lower => { + let x_l = locate_lower_index_cast(grid, point); + if point == cast(*grid.last().unwrap()) { + grid.len() - 1 + } else if point == cast(grid[x_l + 1]) { + x_l + 1 + } else { + x_l + } + } + StepDirection::Upper => { + if point == cast(*grid.first().unwrap()) { + 0 + } else { + locate_lower_index_cast(grid, point) + 1 + } + } + } +} + /// Returns the exact grid index if `point` lies on `grid[lower]` or `grid[lower+1]`, else `None`. /// /// Used to short-circuit interpolation when a query point coincides with a grid coordinate. diff --git a/tests/serde_strategies.rs b/tests/serde_strategies.rs index e67b6a7..20f1fee 100644 --- a/tests/serde_strategies.rs +++ b/tests/serde_strategies.rs @@ -71,4 +71,27 @@ macro_rules! enum_round_trip_test { enum_round_trip_test!(strategy_1d_enum_round_trips_every_variant, Strategy1DEnum); enum_round_trip_test!(strategy_2d_enum_round_trips_every_variant, Strategy2DEnum); enum_round_trip_test!(strategy_3d_enum_round_trips_every_variant, Strategy3DEnum); -enum_round_trip_test!(strategy_nd_enum_round_trips_every_variant, StrategyNDEnum); + +/// Same coverage as [`enum_round_trip_test!`], hand-written for `StrategyNDEnum`: +/// it's the one enum not covered by that macro, since it's non-generic (the Tg/Tv +/// split prototype, see `src/strategy/enums/n.rs`, hardcodes every ND strategy's +/// blend precision to `f64` rather than parameterizing the enum by it), unlike +/// `Strategy1D`/`2D`/`3DEnum`. +#[test] +fn strategy_nd_enum_round_trips_every_variant() { + round_trip(&StrategyNDEnum::from(Nearest)); + round_trip(&StrategyNDEnum::from(Linear)); + round_trip(&StrategyNDEnum::from(LinearUniform)); + round_trip(&StrategyNDEnum::from(Step::lower())); + round_trip(&StrategyNDEnum::from(Step::upper())); + round_trip(&StrategyNDEnum::from(Step::new(vec![ + StepDirection::Lower, + StepDirection::Upper, + ]))); + for bc in cubic_c2_variants() { + round_trip(&StrategyNDEnum::from(bc)); + } + let inner: Box = Box::new(StrategyNDEnum::from(Linear)); + round_trip(&StrategyNDEnum::from(GridTransform::log(inner.clone()))); + round_trip(&StrategyNDEnum::from(ValuesTransform::log(inner))); +}