implement ndg evaluate

Author: espdev

src/ndg.rs

@@ -1,6 +1,7 @@
 mod validate;
 mod make;
 mod evaluate;
+mod util;
 
 use ndarray::{
     NdFloat,
@@ -191,7 +192,10 @@ impl<'a, T, D> GridCubicSmoothingSpline<'a, T, D>
     ///
     pub fn evaluate(&self, xi: &[ArrayView1<'a, T>]) -> Result<Array<T, D>> {
         let xi = xi.to_vec();
-        self.evaluate_spline(&xi)
+        self.evaluate_validate(&xi)?;
+        let yi = self.evaluate_spline(&xi);
+
+        Ok(yi)
     }
 
     /// Returns the ref to smoothing parameters vector or None

src/ndg/evaluate.rs

@@ -1,8 +1,8 @@
 use ndarray::{NdFloat, Dimension, Array, ArrayView1};
 use almost::AlmostEqual;
 
-use crate::Result;
-use super::{NdGridSpline, GridCubicSmoothingSpline};
+use crate::{NdSpline, ndarrayext::to_2d_simple};
+use super::{NdGridSpline, GridCubicSmoothingSpline, util::permute_axes};
 
 
 impl<'a, T, D> NdGridSpline<'a, T, D>
@@ -12,7 +12,42 @@ impl<'a, T, D> NdGridSpline<'a, T, D>
 {
     /// Implements evaluating the spline on the given mesh of Xi-sites
     pub(super) fn evaluate_spline(&self, xi: &[ArrayView1<'a, T>]) -> Array<T, D> {
-        unimplemented!();
+        let mut coeffs = self.coeffs.to_owned();
+        let mut coeffs_shape = coeffs.shape().to_vec();
+
+        let ndim_m1 = self.ndim - 1;
+        let permuted_axes = permute_axes::<D>(self.ndim);
+
+        for ax in (0..self.ndim).rev() {
+            let xi_ax = xi[ax];
+
+            let coeffs_2d = {
+                let coeffs_2d = to_2d_simple(coeffs.view()).unwrap();
+
+                NdSpline::evaluate_spline(
+                    self.order[ax],
+                    self.pieces[ax],
+                    self.breaks[ax],
+                    coeffs_2d,
+                    xi_ax,
+                )
+            };
+
+            let mut shape = D::zeros(self.ndim);
+            shape[ndim_m1] = xi_ax.len();
+            for i in 0..ndim_m1 {
+                shape[i] = coeffs_shape[i];
+            }
+
+            coeffs = coeffs_2d
+                .into_shape(shape).unwrap()
+                .permuted_axes(permuted_axes.clone())
+                .to_owned();
+
+            coeffs_shape = coeffs.shape().to_vec();
+        }
+
+        coeffs
     }
 }
 
@@ -22,7 +57,7 @@ impl<'a, T, D> GridCubicSmoothingSpline<'a, T, D>
         T: NdFloat + AlmostEqual + Default,
         D: Dimension
 {
-    pub(super) fn evaluate_spline(&self, xi: &[ArrayView1<'a, T>]) -> Result<Array<T, D>> {
-        unimplemented!();
+    pub(super) fn evaluate_spline(&self, xi: &[ArrayView1<'a, T>]) -> Array<T, D> {
+        self.spline.as_ref().unwrap().evaluate_spline(&xi)
     }
 }

src/ndg/make.rs

@@ -4,7 +4,7 @@ use almost::AlmostEqual;
 use crate::{Result, CubicSmoothingSpline};
 use crate::ndarrayext::to_2d_simple;
 
-use super::{GridCubicSmoothingSpline, NdGridSpline};
+use super::{GridCubicSmoothingSpline, NdGridSpline, util::permute_axes};
 
 
 impl<'a, T, D> GridCubicSmoothingSpline<'a, T, D>
@@ -22,11 +22,7 @@ impl<'a, T, D> GridCubicSmoothingSpline<'a, T, D>
 
         let mut smooth: Vec<Option<T>> = vec![None; ndim];
 
-        let mut permute_axes = D::zeros(ndim);
-        permute_axes[0] = ndim_m1;
-        for ax in 0..ndim_m1 {
-            permute_axes[ax+1] = ax;
-        }
+        let permuted_axes = permute_axes::<D>(ndim);
 
         for ax in (0..ndim).rev() {
             let x = breaks[ax].view();
@@ -53,7 +49,7 @@ impl<'a, T, D> GridCubicSmoothingSpline<'a, T, D>
 
                 spline.coeffs()
                     .into_shape(new_shape).unwrap()
-                    .permuted_axes(permute_axes.clone())
+                    .permuted_axes(permuted_axes.clone())
                     .to_owned()
             };
 

src/ndg/util.rs

@@ -0,0 +1,16 @@
+use ndarray::Dimension;
+
+
+pub(super) fn permute_axes<D>(ndim: usize) -> D
+    where
+        D: Dimension
+{
+    let mut permute_axes = D::zeros(ndim);
+
+    permute_axes[0] = ndim - 1;
+    for ax in 0..(ndim - 1) {
+        permute_axes[ax + 1] = ax;
+    }
+
+    permute_axes
+}

src/ndg/validate.rs

@@ -25,6 +25,22 @@ impl<'a, T, D> GridCubicSmoothingSpline<'a, T, D>
 
         Ok(())
     }
+
+    pub(super) fn evaluate_validate(&self, xi: &[ArrayView1<'a, T>]) -> Result<()> {
+        let x_len = self.x.len();
+        let xi_len = xi.len();
+
+        if xi_len != x_len {
+            return Err(
+                InvalidInputData(
+                    format!("The number of `xi` vectors ({}) is not equal to the number of dimensions ({})",
+                            xi_len, x_len)
+                )
+            )
+        }
+
+        Ok(())
+    }
 }
 
 

src/umv.rs

@@ -94,7 +94,13 @@ impl<'a, T> NdSpline<'a, T>
 
     /// Evaluates the spline on the given data sites
     pub fn evaluate(&self, xi: ArrayView1<'a, T>) -> Array2<T> {
-        self.evaluate_spline(xi)
+        Self::evaluate_spline(
+            self.order,
+            self.pieces,
+            self.breaks.view(),
+            self.coeffs.view(),
+            xi,
+        )
     }
 }
 
@@ -278,7 +284,7 @@ impl<'a, T, D> CubicSmoothingSpline<'a, T, D>
     /// - If reshaping Y data to 2-d view has failed
     ///
     pub fn make(mut self) -> Result<Self> {
-        self.make_validate_data()?;
+        self.make_validate()?;
         self.make_spline()?;
         Ok(self)
     }
@@ -294,7 +300,7 @@ impl<'a, T, D> CubicSmoothingSpline<'a, T, D>
         where X: AsArray<'a, T>
     {
         let xi = xi.into();
-        self.evaluate_validate_data(xi)?;
+        self.evaluate_validate(xi)?;
 
         let yi = self.evaluate_spline(xi)?;
         Ok(yi)

src/umv/evaluate.rs

@@ -1,4 +1,4 @@
-use ndarray::{NdFloat, Dimension, Array, Array1, Array2, ArrayView1, Axis, s, stack};
+use ndarray::{NdFloat, Dimension, Array, Array1, Array2, ArrayView1, Axis, s, stack, ArrayView2};
 use almost::AlmostEqual;
 
 use crate::{Result, ndarrayext};
@@ -9,10 +9,16 @@ impl<'a, T> NdSpline<'a, T>
     where T: NdFloat + AlmostEqual
 {
     /// Implements evaluating the spline on the given mesh of Xi-sites
-    pub(super) fn evaluate_spline(&self, xi: ArrayView1<'a, T>) -> Array2<T> {
+    pub(crate) fn evaluate_spline(
+        order: usize,
+        pieces: usize,
+        breaks: ArrayView1<'_, T>,
+        coeffs: ArrayView2<'_, T>,
+        xi: ArrayView1<'a, T>) -> Array2<T>
+    {
 
         let edges = {
-            let mesh = self.breaks.slice(s![1..-1]);
+            let mesh = breaks.slice(s![1..-1]);
             let one = Array1::<T>::ones((1, ));
             let left_bound = &one * T::neg_infinity();
             let right_bound = &one * T::infinity();
@@ -24,7 +30,7 @@ impl<'a, T> NdSpline<'a, T>
 
         // Go to local coordinates
         let xi = {
-            let indexed_breaks = indices.mapv(|i| self.breaks[i]);
+            let indexed_breaks = indices.mapv(|i| breaks[i]);
             &xi - &indexed_breaks
         };
 
@@ -35,7 +41,7 @@ impl<'a, T> NdSpline<'a, T>
         let get_indexed_coeffs = |inds: &Array1<usize>| {
             // Returns Nx1 2-d array of coeffs by given index
             let coeffs_by_index = |&index| {
-                self.coeffs.slice(s![.., index]).insert_axis(Axis(1))
+                coeffs.slice(s![.., index]).insert_axis(Axis(1))
             };
 
             // Get the M-sized vector of coeffs values Nx1 arrays
@@ -51,8 +57,8 @@ impl<'a, T> NdSpline<'a, T>
         // Vectorized computing the spline pieces (polynoms) on the given data sites
         let mut values = get_indexed_coeffs(&indices);
 
-        for _ in 1..self.order {
-            indices += self.pieces;
+        for _ in 1..order {
+            indices += pieces;
             values = values * &xi + get_indexed_coeffs(&indices);
         }
 
@@ -75,7 +81,7 @@ impl<'a, T, D> CubicSmoothingSpline<'a, T, D>
             shape[i] = s
         }
 
-        let yi_2d = self.spline.as_ref().unwrap().evaluate_spline(xi);
+        let yi_2d = self.spline.as_ref().unwrap().evaluate(xi);
         let yi = ndarrayext::from_2d(&yi_2d, shape, axis)?.to_owned();
 
         Ok(yi)

src/umv/validate.rs

@@ -19,7 +19,7 @@ use crate::{
 impl<'a, T, D> CubicSmoothingSpline<'a, T, D>
     where T: NdFloat + Default + AlmostEqual, D: Dimension
 {
-    pub(super) fn make_validate_data(&self) -> Result<()> {
+    pub(super) fn make_validate(&self) -> Result<()> {
         let x_size = self.x.len();
 
         if x_size < 2 {
@@ -80,7 +80,15 @@ impl<'a, T, D> CubicSmoothingSpline<'a, T, D>
         Ok(())
     }
 
-    pub(super) fn evaluate_validate_data(&self, xi: ArrayView1<'a, T>) -> Result<()> {
+    pub(super) fn evaluate_validate(&self, xi: ArrayView1<'a, T>) -> Result<()> {
+        if xi.len() < 1 {
+            return Err(
+                InvalidInputData(
+                    "The size of data vectors must be greater or equal to 1".to_string()
+                )
+            )
+        }
+
         if self.spline.is_none() {
             return Err(
                 InvalidInputData(

tests/ndg_evaluate.rs

@@ -0,0 +1,25 @@
+use ndarray::array;
+use approx::assert_abs_diff_eq;
+
+use csaps::GridCubicSmoothingSpline;
+
+
+#[test]
+fn test_make_surface() {
+    let x0 = array![1., 2., 3.];
+    let x1 = array![1., 2., 3., 4.];
+
+    let x = vec![x0.view(), x1.view()];
+
+    let y = array![
+        [1., 2., 3., 4.],
+        [5., 6., 7., 8.],
+        [9., 10., 11., 12.],
+    ];
+
+    let yi = GridCubicSmoothingSpline::new(&x, &y)
+        .make().unwrap()
+        .evaluate(&x).unwrap();
+
+    assert_abs_diff_eq!(yi, y);
+}