diff --git a/Cargo.toml b/Cargo.toml index 6d8aea4..74694df 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,3 +3,11 @@ members = ["ndinterp", "ndinterp_capi"] # "ndinterp_py"] [profile.release] lto = true + +[workspace.package] +authors = ["Alessandro Candido ", "Juan Cruz-Martinez ", "Christopher Schwan "] +edition = "2021" +license = "GPL-3.0-or-later" +repository = "https://github.com/AleCandido/ndinterp" +categories = ["science"] +description = "N-dimensional interpolation library" diff --git a/ndinterp/Cargo.toml b/ndinterp/Cargo.toml index 0157157..43328ae 100644 --- a/ndinterp/Cargo.toml +++ b/ndinterp/Cargo.toml @@ -1,16 +1,18 @@ [package] name = "ndinterp" version = "0.0.1" -authors = ["Alessandro Candido "] -edition = "2021" -license = "GPL-3.0-or-later" -repository = "https://github.com/AleCandido/ndinterp" +edition.workspace = true +authors.workspace = true +license.workspace = true +repository.workspace = true +categories.workspace = true readme = "README.md" -categories = ["science"] description = "N-dimensional interpolation library" +keywords = ["math", "science"] [dependencies] ndarray = "0.15.4" ndarray-linalg = "0.14.1" petgraph = "0.6.2" thiserror = "1.0.40" +itertools = "0.11.0" diff --git a/ndinterp/examples/alphas.rs b/ndinterp/examples/alphas.rs index 8017ca2..86f75ed 100644 --- a/ndinterp/examples/alphas.rs +++ b/ndinterp/examples/alphas.rs @@ -40,7 +40,7 @@ fn main() { ]; let grid = Grid { - input: logq2, + xgrid: vec![logq2.to_vec()], values: alpha_s_vals, }; diff --git a/ndinterp/src/grid.rs b/ndinterp/src/grid.rs index a4c67f8..cad9507 100644 --- a/ndinterp/src/grid.rs +++ b/ndinterp/src/grid.rs @@ -10,19 +10,34 @@ //! y = f(x1, x2, x3...) //! use crate::interpolate::InterpolationError; -use ndarray::Array1; +use itertools::izip; +use ndarray::{Array, ArrayView1, Dimension}; // 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, - pub values: Array1, +pub struct Grid { + /// Arrays with the input vectors (x_i) + pub xgrid: Vec>, + + /// Output points + pub values: Array, +} + +/// A grid slice is always 1-Dimensional +#[derive(Debug)] +pub struct GridSlice<'a> { + /// A reference to one of the input vectors of the grid + pub x: &'a Vec, + /// A view of the slice of values corresponding to x + pub y: ArrayView1<'a, f64>, } -impl Grid { +impl<'a> GridSlice<'a> { // 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 @@ -33,9 +48,9 @@ impl Grid { /// input at position index as the ratio between the differences dy/dx computed as: /// 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 dy = self.values[index] - self.values[index - 1]; + pub fn derivative_at(&'a self, index: usize) -> f64 { + let dx = self.x[index] - self.x[index - 1]; + let dy = self.y[index] - self.y[index - 1]; dy / dx } @@ -44,53 +59,66 @@ impl Grid { /// /// Dx_{i} = \Delta x_{i} = x_{i} - x_{i-} /// y'_{i} = 1/2 * ( (y_{i+1}-y_{i})/Dx_{i+1} + (y_{i}-y_{i-1})/Dx_{i} ) - pub fn central_derivative_at(&self, index: usize) -> f64 { + pub fn central_derivative_at(&'a self, index: usize) -> f64 { let dy_f = self.derivative_at(index + 1); let dy_b = self.derivative_at(index); 0.5 * (dy_f + dy_b) } +} - /// Find the index of the last value in the input such that input(idx) < query +impl Grid { + /// Find the index of the last value in the input xgrid such that xgrid(idx) < query /// If the query is outside the grid returns an extrapolation error - pub fn closest_below(&self, query: f64) -> Result { - if query > self.input[self.input.len() - 1] { - Err(InterpolationError::ExtrapolationAbove(query)) - } else if query < self.input[0] { - Err(InterpolationError::ExtrapolationBelow(query)) - } else { - let u_idx = self.input.iter().position(|x| x > &query).unwrap(); - let idx = u_idx - 1; - Ok(idx) + pub fn closest_below( + &self, + input_query: &[f64], + ) -> Result<[usize; N], InterpolationError> { + let mut ret = [0; N]; + + for (r, &query, igrid) in izip!(&mut ret, input_query, &self.xgrid) { + if query > *igrid.last().unwrap() { + return Err(InterpolationError::ExtrapolationAbove(query)); + } else if query < igrid[0] { + return Err(InterpolationError::ExtrapolationBelow(query)); + } + + let u_idx = igrid.partition_point(|&x| x < query); + *r = u_idx - 1; } + Ok(ret) } } #[cfg(test)] mod tests { use super::*; - use ndarray::array; + use ndarray::{array, Ix1}; - fn gen_grid() -> Grid { - let x = array![0., 1., 2., 3., 4.]; + fn gen_grid() -> Grid { + let x = vec![vec![0., 1., 2., 3., 4.]]; let y = array![4., 3., 2., 1., 1.]; - let grid = Grid { - input: x, + + Grid { + xgrid: x, values: y, - }; - grid + } } #[test] fn check_derivative() { let grid = gen_grid(); - assert_eq!(grid.central_derivative_at(1), -1.); - assert_eq!(grid.central_derivative_at(3), -0.5); + let grid_slice = GridSlice { + x: &grid.xgrid[0], + y: grid.values.view(), + }; + assert_eq!(grid_slice.central_derivative_at(1), -1.); + assert_eq!(grid_slice.central_derivative_at(3), -0.5); } #[test] fn check_index_search() { let grid = gen_grid(); - assert_eq!(grid.closest_below(0.5).unwrap(), 0); - assert_eq!(grid.closest_below(3.2).unwrap(), 3); + assert_eq!(grid.closest_below::<1>(&[0.5]).unwrap()[0], 0); + assert_eq!(grid.closest_below::<1>(&[3.2]).unwrap()[0], 3); } } diff --git a/ndinterp/src/grid/cubic.rs b/ndinterp/src/grid/cubic.rs index 488f922..1d59852 100644 --- a/ndinterp/src/grid/cubic.rs +++ b/ndinterp/src/grid/cubic.rs @@ -1,8 +1,9 @@ //! Implements cubic interpolation algorithms -use crate::grid::Grid; +use crate::grid::{Grid, GridSlice}; use crate::interpolate::InterpolationError; pub use crate::interpolate::Interpolator; +use ndarray::Ix1; /// Cubic interpolation in 1D /// @@ -18,40 +19,35 @@ 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, } -impl Interpolator for Cubic1d { - /// Use Cubic interpolation 1d to compute y(query) - /// The interpolation uses the two nearest neighbours and their derivatives computed as an - /// average of the differences above and below. - /// - /// Special cases are considered when the interpolation occurs between the first (last) two - /// bins, where the derivative would involve points outside the grids. - /// - /// Two special are considered, when the interpolation occurs between the first (last) two - /// bins, the derivative at the boundary is approximated by the forward (backward) difference - fn interpolate(&self, query: f64) -> Result { - let idx = self.grid.closest_below(query)?; - let dx = self.grid.input[idx + 1] - self.grid.input[idx]; +impl<'a> GridSlice<'a> { + /// Implements utilities for a GridSlice that can be used by cubic interpolation Nd + /// Takes as input the value being queried and its index within the given slice + fn cubic_interpolate_1d(&'a self, query: f64, idx: usize) -> Result { + // grid slice utilities are expected to be called multipled times for the same + // query and so it is convient to pass idx from the outside to avoid expensive searches + let dx = self.x[idx + 1] - self.x[idx]; // Upper and lower bounds and derivatives - let yu = self.grid.values[idx + 1]; - let yl = self.grid.values[idx]; + let yu = self.y[idx + 1]; + let yl = self.y[idx]; - let dydxu = if idx == self.grid.input.len() - 2 { - dx * self.grid.derivative_at(idx + 1) + let dydxu = if idx == self.x.len() - 2 { + dx * self.derivative_at(idx + 1) } else { - dx * self.grid.central_derivative_at(idx + 1) + dx * self.central_derivative_at(idx + 1) }; let dydxl = if idx == 0 { - dx * self.grid.derivative_at(idx + 1) + dx * self.derivative_at(idx + 1) } else { - dx * self.grid.central_derivative_at(idx) + dx * self.central_derivative_at(idx) }; - let t = (query - self.grid.input[idx]) / dx; + let t = (query - self.x[idx]) / dx; let t2 = t * t; let t3 = t2 * t; @@ -63,3 +59,26 @@ impl Interpolator for Cubic1d { Ok(p0 + p1 + m0 + m1) } } + +impl Interpolator for Cubic1d { + /// Use Cubic interpolation 1d to compute y(query) + /// The interpolation uses the two nearest neighbours and their derivatives computed as an + /// average of the differences above and below. + /// + /// Special cases are considered when the interpolation occurs between the first (last) two + /// bins, where the derivative would involve points outside the grids. + /// + /// Two special are considered, when the interpolation occurs between the first (last) two + /// bins, the derivative at the boundary is approximated by the forward (backward) difference + fn interpolate(&self, query: f64) -> Result { + let raw_idx = self.grid.closest_below::<1>(&[query])?; + let idx = raw_idx[0]; + + let grid_sl = GridSlice { + x: &self.grid.xgrid[0], + y: self.grid.values.view(), + }; + + grid_sl.cubic_interpolate_1d(query, idx) + } +} diff --git a/ndinterp/src/interpolate.rs b/ndinterp/src/interpolate.rs index dfe548b..f858fb3 100644 --- a/ndinterp/src/interpolate.rs +++ b/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 { + /// Produce the result of the inteprolation given a (nd) point 'query' fn interpolate(&self, query: T) -> Result; } -/// ---- deal with the stuff below later ---- -pub trait Interpolate { - type Point: Metric; - - fn interpolate(&self, query: &Self::Point) -> f64; -} - -pub struct Input { - pub point: Point, - pub value: f64, -} - -impl From<(Point, f64)> for Input { - fn from(item: (Point, f64)) -> Self { - Self { - point: item.0, - value: item.1, - } - } -} - -impl Input { - pub fn stack(points: Vec, values: Vec) -> Vec { - 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 { +// pub point: Point, +// pub value: f64, +//} +// +//impl From<(Point, f64)> for Input { +// fn from(item: (Point, f64)) -> Self { +// Self { +// point: item.0, +// value: item.1, +// } +// } +//} +// +//impl Input { +// pub fn stack(points: Vec, values: Vec) -> Vec { +// zip(points.into_iter(), values.into_iter()) +// .map(|t| t.into()) +// .collect() +// } +//} diff --git a/ndinterp/src/lib.rs b/ndinterp/src/lib.rs index 0915396..e3c2911 100644 --- a/ndinterp/src/lib.rs +++ b/ndinterp/src/lib.rs @@ -1,10 +1,8 @@ -#![warn(clippy::all, clippy::cargo, clippy::nursery, clippy::pedantic)] +#![warn(clippy::all, clippy::cargo)] #![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; diff --git a/ndinterp_capi/Cargo.toml b/ndinterp_capi/Cargo.toml index 8da9009..9a33db6 100644 --- a/ndinterp_capi/Cargo.toml +++ b/ndinterp_capi/Cargo.toml @@ -1,14 +1,21 @@ [package] name = "ndinterp_capi" version = "0.0.1" -edition = "2021" +edition.workspace = true +authors.workspace = true +license.workspace = true +repository.workspace = true +categories.workspace = true +readme = "README.md" +description = "C bindings for ndinterp" +keywords = ["math", "science"] [dependencies] ndarray = "0.15.4" ndinterp = { path = "../ndinterp/" } [lib] -name = "ndinterp" +name = "ndinterp_capi" crate-type = ["cdylib"] [features] diff --git a/ndinterp_capi/benches/Makefile b/ndinterp_capi/benches/Makefile index 43e8360..78ba1cc 100644 --- a/ndinterp_capi/benches/Makefile +++ b/ndinterp_capi/benches/Makefile @@ -1,14 +1,14 @@ CXX = c++ CXXFLAGS = -std=c++17 # Assuming the lib has been installed with cinstall -NDFLAGS = $(shell pkg-config ndinterp --libs) $(shell pkg-config ndinterp --cflags) +NDFLAGS = $(shell pkg-config ndinterp_capi --libs) $(shell pkg-config ndinterp_capi --cflags) LHAFLAGS = $(shell lhapdf-config --libs --cflags) lharun: lhacheck ./lhacheck lhacheck: lhacheck.cpp - $(CXX) $(CXXFLAGS) $< $(NDFLAGS) $(LHAFLAGS) -o $@ + $(CXX) $(CXXFLAGS) -g $< $(NDFLAGS) $(LHAFLAGS) -o $@ clean: rm -f lhacheck diff --git a/ndinterp_capi/benches/lhacheck.cpp b/ndinterp_capi/benches/lhacheck.cpp index 47b5150..d0083d7 100644 --- a/ndinterp_capi/benches/lhacheck.cpp +++ b/ndinterp_capi/benches/lhacheck.cpp @@ -1,5 +1,5 @@ #include "LHAPDF/LHAPDF.h" -#include "ndinterp.h" +#include "ndinterp_capi.h" #include #include #include diff --git a/ndinterp_capi/src/lib.rs b/ndinterp_capi/src/lib.rs index 654aa3e..6b0e251 100644 --- a/ndinterp_capi/src/lib.rs +++ b/ndinterp_capi/src/lib.rs @@ -1,11 +1,14 @@ -#![warn(clippy::all, clippy::cargo, clippy::nursery, clippy::pedantic)] +//! C interface for ndinterp +#![warn(clippy::all, clippy::cargo)] #![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 @@ -13,19 +16,22 @@ pub struct Cubic1d; /// /// # Safety /// -/// This function is only safe to call as long as `input_c` and `values_c` are arrays with sizes +/// This function is only safe to call as long as `xgrid_c` and `values_c` are arrays with sizes /// larger or equal to `size`. #[no_mangle] pub unsafe extern "C" fn create_cubic_interpolator1d( - input_c: *const f64, + xgrid_c: *const f64, values_c: *const f64, size: usize, ) -> Box { - let input = ArrayView1::from_shape_ptr(size, input_c); + // Use slice instead of vec, so that rust doesn't release the memory coming from C++ + let slice_input = unsafe { slice::from_raw_parts(xgrid_c, size) }; + // Make a copy of the data into a vector (of vectors) for rust to own + let xgrid = vec![slice_input.to_vec()]; let values = ArrayView1::from_shape_ptr(size, values_c); let grid = grid::Grid { - input: input.into_owned(), + xgrid, values: values.into_owned(), }; let cubic_interpolator = grid::cubic::Cubic1d { grid };