Update the code for n-d

ndinterp/examples/alphas.rs

@@ -40,7 +40,7 @@ fn main() {
     ];
 
     let grid = Grid {
-        input: logq2,
+        input: vec![logq2.to_vec()],
         values: alpha_s_vals,
     };
 

ndinterp/src/grid.rs

@@ -10,19 +10,24 @@
 //!     y = f(x1, x2, x3...)
 //!
 use crate::interpolate::InterpolationError;
-use ndarray::Array1;
+use ndarray::{Array, Dimension, Ix1};
 
 // Make public the families of interpolation algorithms implemented for grids
 pub mod cubic;
 
+/// A grid is made of two components:
+///     A d-dimensional vector of 1-dimensional sorted vectors for the input points
+///     A d-dimensional array for the grid values of
 #[derive(Debug)]
-pub struct Grid {
-    /// A grid is made of two (1-dimensional) sorted arrays.
-    pub input: Array1<f64>,
-    pub values: Array1<f64>,
+pub struct Grid<D: Dimension> {
+    /// Arrays with the input vectors
+    pub input: Vec<Vec<f64>>,
+
+    /// Output points
+    pub values: Array<f64, D>,
 }
 
-impl Grid {
+impl Grid<Ix1> {
     // TODO: at the moment we are using here the derivatives that LHAPDF is using for the
     // interpolation in alpha_s, these are probably enough for this use case but not in general
     // - [ ] Implement a more robust form of the derivative
@@ -34,7 +39,7 @@ impl Grid {
     ///     dy = y_{i} - y_{i-1}
     ///     dx = x_{i} - x_{x-1}
     pub fn derivative_at(&self, index: usize) -> f64 {
-        let dx = self.input[index] - self.input[index - 1];
+        let dx = self.input[0][index] - self.input[0][index - 1];
         let dy = self.values[index] - self.values[index - 1];
         dy / dx
     }
@@ -53,12 +58,13 @@ impl Grid {
     /// Find the index of the last value in the input such that input(idx) < query
     /// If the query is outside the grid returns an extrapolation error
     pub fn closest_below(&self, query: f64) -> Result<usize, InterpolationError> {
-        if query > self.input[self.input.len() - 1] {
+        if query > self.input[0][self.input[0].len() - 1] {
             Err(InterpolationError::ExtrapolationAbove(query))
-        } else if query < self.input[0] {
+        } else if query < self.input[0][0] {
             Err(InterpolationError::ExtrapolationBelow(query))
         } else {
-            let u_idx = self.input.iter().position(|x| x > &query).unwrap();
+            let u_idx = self.input[0].partition_point(|&x| x < query);
+            //            let u_idx = self.input.iter().position(|x| x > &query).unwrap();
             let idx = u_idx - 1;
             Ok(idx)
         }
@@ -70,8 +76,8 @@ mod tests {
     use super::*;
     use ndarray::array;
 
-    fn gen_grid() -> Grid {
-        let x = array![0., 1., 2., 3., 4.];
+    fn gen_grid() -> Grid<Ix1> {
+        let x = vec![vec![0., 1., 2., 3., 4.]];
         let y = array![4., 3., 2., 1., 1.];
         let grid = Grid {
             input: x,

ndinterp/src/grid/cubic.rs

@@ -3,6 +3,7 @@
 use crate::grid::Grid;
 use crate::interpolate::InterpolationError;
 pub use crate::interpolate::Interpolator;
+use ndarray::Ix1;
 
 /// Cubic interpolation in 1D
 ///
@@ -18,7 +19,8 @@ pub use crate::interpolate::Interpolator;
 /// with hij the Hermite basis functions
 #[derive(Debug)]
 pub struct Cubic1d {
-    pub grid: Grid,
+    /// The grid object contains all necessary information to perform the interpolation
+    pub grid: Grid<Ix1>,
 }
 
 impl Interpolator<f64> for Cubic1d {
@@ -33,13 +35,13 @@ impl Interpolator<f64> for Cubic1d {
     /// bins, the derivative at the boundary is approximated by the forward (backward) difference
     fn interpolate(&self, query: f64) -> Result<f64, InterpolationError> {
         let idx = self.grid.closest_below(query)?;
-        let dx = self.grid.input[idx + 1] - self.grid.input[idx];
+        let dx = self.grid.input[0][idx + 1] - self.grid.input[0][idx];
 
         // Upper and lower bounds and derivatives
         let yu = self.grid.values[idx + 1];
         let yl = self.grid.values[idx];
 
-        let dydxu = if idx == self.grid.input.len() - 2 {
+        let dydxu = if idx == self.grid.input[0].len() - 2 {
             dx * self.grid.derivative_at(idx + 1)
         } else {
             dx * self.grid.central_derivative_at(idx + 1)
@@ -51,7 +53,7 @@ impl Interpolator<f64> for Cubic1d {
             dx * self.grid.central_derivative_at(idx)
         };
 
-        let t = (query - self.grid.input[idx]) / dx;
+        let t = (query - self.grid.input[0][idx]) / dx;
         let t2 = t * t;
         let t3 = t2 * t;
 

ndinterp/src/interpolate.rs

@@ -1,46 +1,49 @@
-use std::iter::zip;
+//! This module implements interpolation rutines
 use thiserror::Error;
 
-use crate::metric::Metric;
-
+/// Errors encountered during interpolation
 #[derive(Debug, Error)]
 pub enum InterpolationError {
+    /// Raised when the queried value is above the maximum
     #[error("The value queried ({0}) is above the maximum")]
     ExtrapolationAbove(f64),
 
+    /// Raised when the queried value is below the minimum
     #[error("The value queried ({0}) is below the minimum")]
     ExtrapolationBelow(f64),
 }
 
+/// Methods which all interpolator must implement
 pub trait Interpolator<T> {
+    /// Produce the result of the inteprolation given a (nd) point 'query'
     fn interpolate(&self, query: T) -> Result<f64, InterpolationError>;
 }
 
-/// ---- deal with the stuff below later ----
-pub trait Interpolate {
-    type Point: Metric;
-
-    fn interpolate(&self, query: &Self::Point) -> f64;
-}
-
-pub struct Input<Point: Metric> {
-    pub point: Point,
-    pub value: f64,
-}
-
-impl<Point: Metric> From<(Point, f64)> for Input<Point> {
-    fn from(item: (Point, f64)) -> Self {
-        Self {
-            point: item.0,
-            value: item.1,
-        }
-    }
-}
-
-impl<Point: Metric> Input<Point> {
-    pub fn stack(points: Vec<Point>, values: Vec<f64>) -> Vec<Self> {
-        zip(points.into_iter(), values.into_iter())
-            .map(|t| t.into())
-            .collect()
-    }
-}
+///// ---- deal with the stuff below later ----
+//pub trait Interpolate {
+//    type Point: Metric;
+//
+//    fn interpolate(&self, query: &Self::Point) -> f64;
+//}
+//
+//pub struct Input<Point: Metric> {
+//    pub point: Point,
+//    pub value: f64,
+//}
+//
+//impl<Point: Metric> From<(Point, f64)> for Input<Point> {
+//    fn from(item: (Point, f64)) -> Self {
+//        Self {
+//            point: item.0,
+//            value: item.1,
+//        }
+//    }
+//}
+//
+//impl<Point: Metric> Input<Point> {
+//    pub fn stack(points: Vec<Point>, values: Vec<f64>) -> Vec<Self> {
+//        zip(points.into_iter(), values.into_iter())
+//            .map(|t| t.into())
+//            .collect()
+//    }
+//}

ndinterp/src/lib.rs

@@ -1,10 +1,8 @@
 #![warn(clippy::all, clippy::cargo, clippy::nursery, clippy::pedantic)]
 #![warn(missing_docs)]
 #![doc = include_str!("../README.md")]
-#![warn(clippy::all, clippy::pedantic, clippy::restriction, clippy::cargo)]
-#![warn(missing_docs)]
 
 pub mod grid;
 pub mod interpolate;
-pub mod metric;
-pub mod scatter;
+//pub mod metric;
+//pub mod scatter;

ndinterp_capi/benches/Makefile

@@ -8,7 +8,7 @@ lharun: lhacheck
 	./lhacheck
 
 lhacheck: lhacheck.cpp
-	$(CXX) $(CXXFLAGS) $< $(NDFLAGS) $(LHAFLAGS) -o $@
+	$(CXX) $(CXXFLAGS) -g $< $(NDFLAGS) $(LHAFLAGS) -o $@
 
 clean:
 	rm -f lhacheck

ndinterp_capi/src/lib.rs

@@ -1,11 +1,14 @@
+//! C interface for ndinterp
 #![warn(clippy::all, clippy::cargo, clippy::nursery, clippy::pedantic)]
 #![warn(missing_docs)]
 
+use core::slice;
+
 use ndarray::ArrayView1;
-/// C interface for ndinterp
 use ndinterp::grid;
 use ndinterp::interpolate::Interpolator;
 
+/// cubic1d inteprolator
 pub struct Cubic1d;
 
 /// Creates a cubic1d interpolator given the nodes
@@ -21,11 +24,14 @@ pub unsafe extern "C" fn create_cubic_interpolator1d(
     values_c: *const f64,
     size: usize,
 ) -> Box<grid::cubic::Cubic1d> {
-    let input = ArrayView1::from_shape_ptr(size, input_c);
+    // Use slice instead of vec so that rust doesn't take ownership of the data and releases
+    // the burden is on the function calling them
+    let slice_input = unsafe { slice::from_raw_parts(input_c, size) } ;
+    let input = vec![slice_input.to_vec()];
     let values = ArrayView1::from_shape_ptr(size, values_c);
 
     let grid = grid::Grid {
-        input: input.into_owned(),
+        input,
         values: values.into_owned(),
     };
     let cubic_interpolator = grid::cubic::Cubic1d { grid };