From 8c78f6fef081545bc81b5b6f9c4aa6badf7e138a Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Tue, 2 Sep 2025 15:56:58 +0200 Subject: [PATCH 01/21] Crude implementation of Chebyshev batch interpolations --- neopdf/src/gridpdf.rs | 161 +++++++++++++- neopdf/src/pdf.rs | 17 ++ neopdf/src/strategy.rs | 331 ++++++++++++++++++++++++++++- neopdf/tests/pdf.rs | 27 +++ neopdf_capi/src/include/NeoPDF.hpp | 17 ++ neopdf_capi/src/lib.rs | 40 ++++ 6 files changed, 591 insertions(+), 2 deletions(-) diff --git a/neopdf/src/gridpdf.rs b/neopdf/src/gridpdf.rs index 443f4e8c..019b710c 100644 --- a/neopdf/src/gridpdf.rs +++ b/neopdf/src/gridpdf.rs @@ -7,7 +7,8 @@ use core::panic; -use ndarray::{Array1, Array2}; +use ndarray::{s, Array1, Array2}; +use ninterp::data::{InterpData2D, InterpData3D}; use rayon::prelude::*; use serde::{Deserialize, Serialize}; use thiserror::Error; @@ -16,6 +17,7 @@ use super::alphas::AlphaS; use super::interpolator::{DynInterpolator, InterpolatorFactory}; use super::metadata::{InterpolatorType, MetaData}; use super::parser::SubgridData; +use super::strategy::LogChebyshevBatchInterpolation; use super::subgrid::{ParamRange, RangeParameters, SubGrid}; /// Errors that can occur during PDF grid operations. @@ -338,6 +340,163 @@ impl GridPDF { Array2::from_shape_vec(grid_shape, data).unwrap() } + /// Interpolates PDF values for multiple points in parallel using Chebyshev batch interpolation. + /// + /// # Arguments + /// + /// * `flavor_id` - The flavor ID. + /// * `points` - A slice containing the collection of knots to interpolate on. + /// A knot is a collection of points containing `(nucleon, alphas, x, Q2)`. + /// + /// # Returns + /// + /// A `Vec` of interpolated PDF values. + pub fn xfxq2_cheby_batch(&self, flavor_id: i32, points: &[&[f64]]) -> Result, Error> { + if points.is_empty() { + return Ok(Vec::new()); + } + + let subgrid_idx = self.knot_array.find_subgrid(points[0]).ok_or_else(|| { + let (x, q2) = self.get_x_q2(points[0]); + Error::SubgridNotFound { x, q2 } + })?; + + let pid_idx = self.knot_array.pid_index(flavor_id).ok_or_else(|| { + Error::InterpolationError(format!("Invalid flavor ID: {}", flavor_id)) + })?; + + if !matches!(self.info.interpolator_type, InterpolatorType::LogChebyshev) { + return Err(Error::InterpolationError( + "xfxq2_cheby_batch only supports LogChebyshev interpolator".to_string(), + )); + } + + let subgrid = &self.knot_array.subgrids[subgrid_idx]; + let use_log = true; // LogChebyshev is always log + + let log_points: Vec> = points + .iter() + .map(|p| { + p.iter() + .map(|&v| if use_log { v.ln() } else { v }) + .collect() + }) + .collect(); + + let results = match subgrid.interpolation_config() { + super::interpolator::InterpolationConfig::TwoD => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<2>::default(); + let grid_slice = subgrid.grid_slice(pid_idx).to_owned(); + let data = InterpData2D::new( + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + grid_slice, + ) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + batch_interpolator + .init(&data) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + + let points_2d: Vec<[f64; 2]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 2D")) + .collect(); + + batch_interpolator + .interpolate(&data, &points_2d) + .map_err(|e| Error::InterpolationError(e.to_string())) + } + super::interpolator::InterpolationConfig::ThreeDNucleons => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![.., 0, pid_idx, 0, .., ..]).to_owned(); + let reshaped_data = grid_data + .into_shape_with_order(( + subgrid.nucleons.len(), + subgrid.xs.len(), + subgrid.q2s.len(), + )) + .expect("Failed to reshape 3D data"); + let data = InterpData3D::new( + subgrid.nucleons.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + batch_interpolator + .init(&data) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + let points_3d: Vec<[f64; 3]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + batch_interpolator + .interpolate(&data, &points_3d) + .map_err(|e| Error::InterpolationError(e.to_string())) + } + super::interpolator::InterpolationConfig::ThreeDAlphas => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![0, .., pid_idx, 0, .., ..]).to_owned(); + let reshaped_data = grid_data + .into_shape_with_order(( + subgrid.alphas.len(), + subgrid.xs.len(), + subgrid.q2s.len(), + )) + .expect("Failed to reshape 3D data"); + let data = InterpData3D::new( + subgrid.alphas.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + batch_interpolator + .init(&data) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + let points_3d: Vec<[f64; 3]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + batch_interpolator + .interpolate(&data, &points_3d) + .map_err(|e| Error::InterpolationError(e.to_string())) + } + super::interpolator::InterpolationConfig::ThreeDKt => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![0, 0, pid_idx, .., .., ..]).to_owned(); + let reshaped_data = grid_data + .into_shape_with_order((subgrid.kts.len(), subgrid.xs.len(), subgrid.q2s.len())) + .expect("Failed to reshape 3D data"); + let data = InterpData3D::new( + subgrid.kts.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + batch_interpolator + .init(&data) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + let points_3d: Vec<[f64; 3]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + batch_interpolator + .interpolate(&data, &points_3d) + .map_err(|e| Error::InterpolationError(e.to_string())) + } + _ => Err(Error::InterpolationError( + "Unsupported dimension for xfxq2_cheby_batch".to_string(), + )), + }?; + + Ok(results + .into_iter() + .map(|r| self.apply_force_positive(r)) + .collect()) + } + /// Get the values of the momentum fraction `x` and momentum scale `Q2`. /// /// # Arguments diff --git a/neopdf/src/pdf.rs b/neopdf/src/pdf.rs index bd6893ef..c7f9a33f 100644 --- a/neopdf/src/pdf.rs +++ b/neopdf/src/pdf.rs @@ -287,6 +287,23 @@ impl PDF { self.grid_pdf.xfxq2s(pids, slice_points) } + /// Interpolates the PDF value (xf) for multiple points using Chebyshev batch interpolation. + /// + /// Abstraction to the `GridPDF::xfxq2_cheby_batch` method. + /// + /// # Arguments + /// + /// * `pid` - The flavor ID. + /// * `points` - A slice containing the collection of knots to interpolate on. + /// A knot is a collection of points containing `(nucleon, alphas, x, Q2)`. + /// + /// # Returns + /// + /// A `Vec` of interpolated PDF values. + pub fn xfxq2_cheby_batch(&self, pid: i32, points: &[&[f64]]) -> Vec { + self.grid_pdf.xfxq2_cheby_batch(pid, points).unwrap() + } + /// Interpolates the strong coupling constant `alpha_s` for a given Q2. /// /// Abstraction to the `GridPDF::alphas_q2` method. diff --git a/neopdf/src/strategy.rs b/neopdf/src/strategy.rs index f7d1a42c..27fb12f3 100644 --- a/neopdf/src/strategy.rs +++ b/neopdf/src/strategy.rs @@ -14,7 +14,7 @@ //! All interpolation strategies are designed to work with `ninterp`'s data structures and traits, //! ensuring compatibility and extensibility. -use ndarray::{Data, RawDataClone}; +use ndarray::{Array2, Axis, Data, RawDataClone}; use ninterp::data::{InterpData1D, InterpData2D, InterpData3D}; use ninterp::error::{InterpolateError, ValidateError}; use ninterp::strategy::traits::{Strategy1D, Strategy2D, Strategy3D}; @@ -1210,6 +1210,335 @@ where } } +/// Implements a global N-dimensional batch interpolation using Chebyshev polynomials +/// with logarithmic coordinate scaling. +/// +/// This strategy is optimized for interpolating multiple points at once by leveraging +/// matrix operations with `ndarray`. +/// +/// TODO: Potentially merge this with `LogChebyshevInterpolation`. +#[derive(Debug, Clone)] +pub struct LogChebyshevBatchInterpolation { + // Pre-computed weights for the barycentric formula for each dimension. + weights: [Vec; DIM], + // Grid points in the t-domain [-1, 1] for each dimension. + t_coords: [Vec; DIM], +} + +impl Default for LogChebyshevBatchInterpolation { + fn default() -> Self { + Self { + weights: std::array::from_fn(|_| Vec::new()), + t_coords: std::array::from_fn(|_| Vec::new()), + } + } +} + +impl Serialize for LogChebyshevBatchInterpolation { + fn serialize(&self, serializer: S) -> Result + where + S: serde::Serializer, + { + use serde::ser::SerializeStruct; + let mut state = serializer.serialize_struct("LogChebyshevBatchInterpolation", 2)?; + state.serialize_field("weights", &self.weights.as_slice())?; + state.serialize_field("t_coords", &self.t_coords.as_slice())?; + state.end() + } +} + +impl<'de, const DIM: usize> Deserialize<'de> for LogChebyshevBatchInterpolation { + fn deserialize(deserializer: D) -> Result + where + D: serde::Deserializer<'de>, + { + #[derive(Deserialize)] + struct Helper { + weights: Vec>, + t_coords: Vec>, + } + + let helper = Helper::deserialize(deserializer)?; + let weights = helper.weights.try_into().map_err(|v: Vec>| { + serde::de::Error::invalid_length(v.len(), &"an array of the correct length") + })?; + let t_coords = helper.t_coords.try_into().map_err(|v: Vec>| { + serde::de::Error::invalid_length(v.len(), &"an array of the correct length") + })?; + + Ok(Self { weights, t_coords }) + } +} + +impl LogChebyshevBatchInterpolation { + /// Computes the barycentric weights for a given set of Chebyshev points. + /// The formula for the weights is `w_j = (-1)^j * delta_j`, where `delta_j` + /// is 1/2 for the first and last points, and 1 otherwise. + fn compute_barycentric_weights(n: usize) -> Vec { + let mut weights = vec![1.0; n]; + (0..n).for_each(|j| { + if j % 2 == 1 { + weights[j] = -1.0; + } + }); + weights[0] *= 0.5; + + if n > 1 { + weights[n - 1] *= 0.5; + } + + weights + } + + /// Compute barycentric coefficients for multiple points in batch + fn barycentric_coefficients( + t_values: &[f64], + t_coords: &[f64], + weights: &[f64], + ) -> Array2 { + let num_points = t_values.len(); + let num_coords = t_coords.len(); + let mut coeffs = Array2::::zeros((num_points, num_coords)); + + for (p, &t) in t_values.iter().enumerate() { + let mut found_exact = false; + for (j, &t_j) in t_coords.iter().enumerate() { + if (t - t_j).abs() < 1e-15 { + coeffs[[p, j]] = 1.0; + found_exact = true; + break; + } + } + + if !found_exact { + let mut terms = Vec::with_capacity(num_coords); + for (j, &t_j) in t_coords.iter().enumerate() { + terms.push(weights[j] / (t - t_j)); + } + + let sum: f64 = terms.iter().sum(); + + for (j, &term) in terms.iter().enumerate() { + coeffs[[p, j]] = term / sum; + } + } + } + + coeffs + } +} + +impl LogChebyshevBatchInterpolation<1> { + pub fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> + where + D: Data + RawDataClone + Clone, + { + let x_coords = data.grid[0].as_slice().unwrap(); + let n = x_coords.len(); + if n < 2 { + return Err(ValidateError::Other( + "LogChebyshevBatchInterpolation requires at least 2 grid points.".to_string(), + )); + } + + self.t_coords[0] = (0..n) + .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) + .collect(); + + self.weights[0] = Self::compute_barycentric_weights(n); + + Ok(()) + } + + pub fn interpolate( + &self, + data: &InterpData1D, + points: &[[f64; 1]], + ) -> Result, InterpolateError> + where + D: Data + RawDataClone + Clone, + { + let x_coords = data.grid[0].as_slice().unwrap(); + let f_values = data.values.to_owned(); + + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + + let mut t_x_vals = Vec::with_capacity(points.len()); + for &[x] in points { + if x < x_min || x > x_max { + return Err(InterpolateError::Other(format!( + "Input point x={} is outside the grid range [{}, {}]", + x, x_min, x_max + ))); + } + let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; + t_x_vals.push(t_x); + } + + let c_x = Self::barycentric_coefficients(&t_x_vals, &self.t_coords[0], &self.weights[0]); + let results = c_x.dot(&f_values); + + Ok(results.to_vec()) + } +} + +impl LogChebyshevBatchInterpolation<2> { + pub fn init(&mut self, data: &InterpData2D) -> Result<(), ValidateError> + where + D: Data + RawDataClone + Clone, + { + for dim in 0..2 { + let x_coords = data.grid[dim].as_slice().unwrap(); + let n = x_coords.len(); + if n < 2 { + return Err(ValidateError::Other( + "LogChebyshevBatchInterpolation requires at least 2 grid points per dimension." + .to_string(), + )); + } + self.t_coords[dim] = (0..n) + .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) + .collect(); + self.weights[dim] = Self::compute_barycentric_weights(n); + } + + Ok(()) + } + + pub fn interpolate( + &self, + data: &InterpData2D, + points: &[[f64; 2]], + ) -> Result, InterpolateError> + where + D: Data + RawDataClone + Clone, + { + let x_coords = data.grid[0].as_slice().unwrap(); + let y_coords = data.grid[1].as_slice().unwrap(); + + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + let y_min = *y_coords.first().unwrap(); + let y_max = *y_coords.last().unwrap(); + + let mut t_x_vals = Vec::with_capacity(points.len()); + let mut t_y_vals = Vec::with_capacity(points.len()); + + for &[x, y] in points { + if x < x_min || x > x_max { + return Err(InterpolateError::Other(format!( + "Input point x={} is outside the grid range [{}, {}]", + x, x_min, x_max + ))); + } + if y < y_min || y > y_max { + return Err(InterpolateError::Other(format!( + "Input point y={} is outside the grid range [{}, {}]", + y, y_min, y_max + ))); + } + + let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; + let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; + + t_x_vals.push(t_x); + t_y_vals.push(t_y); + } + + let c_x = Self::barycentric_coefficients(&t_x_vals, &self.t_coords[0], &self.weights[0]); + let c_y = Self::barycentric_coefficients(&t_y_vals, &self.t_coords[1], &self.weights[1]); + let v = data.values.to_owned(); + let results = (&c_x.dot(&v) * &c_y).sum_axis(Axis(1)); + + Ok(results.to_vec()) + } +} + +impl LogChebyshevBatchInterpolation<3> { + pub fn init(&mut self, data: &InterpData3D) -> Result<(), ValidateError> + where + D: Data + RawDataClone + Clone, + { + for dim in 0..3 { + let x_coords = data.grid[dim].as_slice().unwrap(); + let n = x_coords.len(); + if n < 2 { + return Err(ValidateError::Other( + "LogChebyshevBatchInterpolation requires at least 2 grid points per dimension." + .to_string(), + )); + } + self.t_coords[dim] = (0..n) + .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) + .collect(); + self.weights[dim] = Self::compute_barycentric_weights(n); + } + Ok(()) + } + + pub fn interpolate( + &self, + data: &InterpData3D, + points: &[[f64; 3]], + ) -> Result, InterpolateError> + where + D: Data + RawDataClone + Clone, + { + let x_coords = data.grid[0].as_slice().unwrap(); + let y_coords = data.grid[1].as_slice().unwrap(); + let z_coords = data.grid[2].as_slice().unwrap(); + + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + let y_min = *y_coords.first().unwrap(); + let y_max = *y_coords.last().unwrap(); + let z_min = *z_coords.first().unwrap(); + let z_max = *z_coords.last().unwrap(); + + let mut t_x_vals = Vec::with_capacity(points.len()); + let mut t_y_vals = Vec::with_capacity(points.len()); + let mut t_z_vals = Vec::with_capacity(points.len()); + + for &[x, y, z] in points { + if x < x_min || x > x_max || y < y_min || y > y_max || z < z_min || z > z_max { + return Err(InterpolateError::Other( + "Point outside grid range".to_string(), + )); + } + t_x_vals.push(2.0 * (x - x_min) / (x_max - x_min) - 1.0); + t_y_vals.push(2.0 * (y - y_min) / (y_max - y_min) - 1.0); + t_z_vals.push(2.0 * (z - z_min) / (z_max - z_min) - 1.0); + } + + let c_x = Self::barycentric_coefficients(&t_x_vals, &self.t_coords[0], &self.weights[0]); + let c_y = Self::barycentric_coefficients(&t_y_vals, &self.t_coords[1], &self.weights[1]); + let c_z = Self::barycentric_coefficients(&t_z_vals, &self.t_coords[2], &self.weights[2]); + + let v = &data.values; + + let num_points = points.len(); + let (nx, ny, nz) = (x_coords.len(), y_coords.len(), z_coords.len()); + + let v_flat = v.to_owned().into_shape_with_order((nx, ny * nz)).unwrap(); + let temp1 = c_x.dot(&v_flat); + let temp1_3d = temp1.into_shape_with_order((num_points, ny, nz)).unwrap(); + + let mut results = Vec::with_capacity(num_points); + for p in 0..num_points { + let temp_slice = temp1_3d.index_axis(Axis(0), p); + let cy_slice = c_y.row(p); + let cz_slice = c_z.row(p); + + let temp2 = cy_slice.dot(&temp_slice); + let result = cz_slice.dot(&temp2); + results.push(result); + } + + Ok(results) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/neopdf/tests/pdf.rs b/neopdf/tests/pdf.rs index dfc8ef9e..12c7051e 100644 --- a/neopdf/tests/pdf.rs +++ b/neopdf/tests/pdf.rs @@ -301,3 +301,30 @@ fn test_pdf_download() { assert!(gluon_as.is_finite()); assert!(gluon_xf.is_finite()); } + +#[test] +#[ignore] +pub fn test_xfxq2_cheby_batch() { + // TODO: Use proper grid with Chebyshev interpolations + let pdf = PDF::load("NNPDF40_nnlo_as_01180", 0); + + let ids: Vec = (-4..=4).filter(|&x| x != 0).collect(); + let xs = [1e-5, 1e-3, 1e-3, 1.0]; + let q2s = [5.0, 10.0, 100.0]; + + let flatten_points: Vec> = xs + .iter() + .flat_map(|&x| q2s.iter().map(move |&q2| vec![x, q2])) + .collect(); + let points_interp: Vec<&[f64]> = flatten_points.iter().map(Vec::as_slice).collect(); + let slice_points: &[&[f64]] = &points_interp; + + for &pid in &ids { + let results_batch = pdf.xfxq2_cheby_batch(pid, slice_points); + let results_single: Vec = slice_points.iter().map(|p| pdf.xfxq2(pid, p)).collect(); + + for (res_b, res_s) in results_batch.iter().zip(results_single.iter()) { + assert!((res_b - res_s).abs() < LOW_PRECISION); + } + } +} diff --git a/neopdf_capi/src/include/NeoPDF.hpp b/neopdf_capi/src/include/NeoPDF.hpp index 66d0c79b..4cc26f7e 100644 --- a/neopdf_capi/src/include/NeoPDF.hpp +++ b/neopdf_capi/src/include/NeoPDF.hpp @@ -157,6 +157,23 @@ class NeoPDF { return neopdf_pdf_xfxq2_nd(this->raw, pid, params.data(), params.size()); } + /** @brief Compute the `xf` value for a generic set of parameters using batch Chebyshev interpolation. */ + std::vector + xfxQ2_cheby_batch(int pid, const std::vector> &points) const { + std::vector c_points(points.size()); + std::vector lengths(points.size()); + for (size_t i = 0; i < points.size(); ++i) { + c_points[i] = points[i].data(); + lengths[i] = points[i].size(); + } + + std::vector results(points.size()); + neopdf_pdf_xfxq2_cheby_batch(this->raw, pid, c_points.data(), lengths.data(), + points.size(), results.data()); + + return results; + } + /** @brief Compute the value of `alphas` at the Q2 value. */ double alphasQ2(double q2) const { return neopdf_pdf_alphas_q2(this->raw, q2); diff --git a/neopdf_capi/src/lib.rs b/neopdf_capi/src/lib.rs index d4b590a9..d154cd44 100644 --- a/neopdf_capi/src/lib.rs +++ b/neopdf_capi/src/lib.rs @@ -315,6 +315,46 @@ pub unsafe extern "C" fn neopdf_pdf_xfxq2_nd( pdf_obj.xfxq2(id, params) } +/// Interpolates PDF values for multiple points in parallel using Chebyshev batch interpolation. +/// +/// # Safety +/// +/// The `pdf` pointer must be a valid pointer to a `NeoPDF` object. +/// The `points` pointer must be a valid pointer to an array of pointers to `c_double`. +/// The `lengths` pointer must be a valid pointer to an array of `usize`. +/// `num_points` must be the correct number of points. +/// The `results` pointer must be valid for writing `num_points` elements. +/// +/// # Panics +/// +/// TODO +#[no_mangle] +pub unsafe extern "C" fn neopdf_pdf_xfxq2_cheby_batch( + pdf: *mut NeoPDFWrapper, + pid: i32, + points: *const *const c_double, + lengths: *const usize, + num_points: usize, + results: *mut c_double, +) { + assert!(!pdf.is_null()); + let pdf_obj = unsafe { &(*pdf).0 }; + + let points_slices: &[*const c_double] = unsafe { slice::from_raw_parts(points, num_points) }; + let lengths_slice: &[usize] = unsafe { slice::from_raw_parts(lengths, num_points) }; + + let rust_points: Vec<&[f64]> = points_slices + .iter() + .zip(lengths_slice) + .map(|(&p, &l)| unsafe { slice::from_raw_parts(p, l) }) + .collect(); + + let res_vec = pdf_obj.xfxq2_cheby_batch(pid, &rust_points); + + let results_slice = unsafe { slice::from_raw_parts_mut(results, res_vec.len()) }; + results_slice.copy_from_slice(&res_vec); +} + /// Clip the interpolated values if they turned out negatives. /// /// # Panics From 55444a1bfe3a73e34657a62c5d42f2c7d47ebf05 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Tue, 2 Sep 2025 16:14:48 +0200 Subject: [PATCH 02/21] Remove boundary check in `LogChebyshevInterpolation` for extrapolation --- neopdf/src/strategy.rs | 50 +++++------------------------------------- 1 file changed, 6 insertions(+), 44 deletions(-) diff --git a/neopdf/src/strategy.rs b/neopdf/src/strategy.rs index 27fb12f3..9474f362 100644 --- a/neopdf/src/strategy.rs +++ b/neopdf/src/strategy.rs @@ -1031,13 +1031,6 @@ where let x_min = *x_coords.first().unwrap(); let x_max = *x_coords.last().unwrap(); - if x < x_min || x > x_max { - return Err(InterpolateError::Other(format!( - "Input point {} is outside the grid range [{}, {}]", - x, x_min, x_max - ))); - } - if (x_max - x_min).abs() < 1e-15 { return Ok(f_values[0]); } @@ -1089,21 +1082,8 @@ where let x_min = *x_coords.first().unwrap(); let x_max = *x_coords.last().unwrap(); - if x < x_min || x > x_max { - return Err(InterpolateError::Other(format!( - "Input point x={} is outside the grid range [{}, {}]", - x, x_min, x_max - ))); - } - let y_min = *y_coords.first().unwrap(); let y_max = *y_coords.last().unwrap(); - if y < y_min || y > y_max { - return Err(InterpolateError::Other(format!( - "Input point y={} is outside the grid range [{}, {}]", - y, y_min, y_max - ))); - } let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; @@ -1160,30 +1140,10 @@ where let x_min = *x_coords.first().unwrap(); let x_max = *x_coords.last().unwrap(); - if x < x_min || x > x_max { - return Err(InterpolateError::Other(format!( - "Input point x={} is outside the grid range [{}, {}]", - x, x_min, x_max - ))); - } - let y_min = *y_coords.first().unwrap(); let y_max = *y_coords.last().unwrap(); - if y < y_min || y > y_max { - return Err(InterpolateError::Other(format!( - "Input point y={} is outside the grid range [{}, {}]", - y, y_min, y_max - ))); - } - let z_min = *z_coords.first().unwrap(); let z_max = *z_coords.last().unwrap(); - if z < z_min || z > z_max { - return Err(InterpolateError::Other(format!( - "Input point z={} is outside the grid range [{}, {}]", - z, z_min, z_max - ))); - } let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; @@ -1425,6 +1385,7 @@ impl LogChebyshevBatchInterpolation<2> { let mut t_x_vals = Vec::with_capacity(points.len()); let mut t_y_vals = Vec::with_capacity(points.len()); + // TODO: Remove this for extrapolation -> speed gain for &[x, y] in points { if x < x_min || x > x_max { return Err(InterpolateError::Other(format!( @@ -1500,6 +1461,7 @@ impl LogChebyshevBatchInterpolation<3> { let mut t_y_vals = Vec::with_capacity(points.len()); let mut t_z_vals = Vec::with_capacity(points.len()); + // TODO: Remove this for extrapolation -> speed gain for &[x, y, z] in points { if x < x_min || x > x_max || y < y_min || y > y_max || z < z_min || z > z_max { return Err(InterpolateError::Other( @@ -1746,19 +1708,19 @@ mod tests { let alphas_cubic = AlphaSCubicInterpolation; // Test within interpolation range - let result = alphas_cubic.interpolate(&data, &[2.25f64.ln()]).unwrap(); // Q=1.5 + let result = alphas_cubic.interpolate(&data, &[2.25f64.ln()]).unwrap(); assert!(result > 0.1 && result < 0.14); // Test at grid point - let result = alphas_cubic.interpolate(&data, &[4.0f64.ln()]).unwrap(); // Q=2.0 + let result = alphas_cubic.interpolate(&data, &[4.0f64.ln()]).unwrap(); assert_close(result, 0.11, EPSILON); // Test extrapolation below range - let result = alphas_cubic.interpolate(&data, &[0.5f64.ln()]).unwrap(); // Q=sqrt(0.5) + let result = alphas_cubic.interpolate(&data, &[0.5f64.ln()]).unwrap(); assert!(result < 0.1); // Test extrapolation above range - let result = alphas_cubic.interpolate(&data, &[30.0f64.ln()]).unwrap(); // Q=sqrt(30) + let result = alphas_cubic.interpolate(&data, &[30.0f64.ln()]).unwrap(); assert_close(result, 0.14, EPSILON); } From e902acee7eb8ac796965cdaf182286260d09af8f Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Tue, 2 Sep 2025 16:39:29 +0200 Subject: [PATCH 03/21] Clean up a bit --- neopdf/src/gridpdf.rs | 124 ++----------------------------------- neopdf/src/interpolator.rs | 117 +++++++++++++++++++++++++++++++++- 2 files changed, 122 insertions(+), 119 deletions(-) diff --git a/neopdf/src/gridpdf.rs b/neopdf/src/gridpdf.rs index 019b710c..4e8e6193 100644 --- a/neopdf/src/gridpdf.rs +++ b/neopdf/src/gridpdf.rs @@ -7,8 +7,7 @@ use core::panic; -use ndarray::{s, Array1, Array2}; -use ninterp::data::{InterpData2D, InterpData3D}; +use ndarray::{Array1, Array2}; use rayon::prelude::*; use serde::{Deserialize, Serialize}; use thiserror::Error; @@ -17,7 +16,6 @@ use super::alphas::AlphaS; use super::interpolator::{DynInterpolator, InterpolatorFactory}; use super::metadata::{InterpolatorType, MetaData}; use super::parser::SubgridData; -use super::strategy::LogChebyshevBatchInterpolation; use super::subgrid::{ParamRange, RangeParameters, SubGrid}; /// Errors that can occur during PDF grid operations. @@ -367,129 +365,19 @@ impl GridPDF { if !matches!(self.info.interpolator_type, InterpolatorType::LogChebyshev) { return Err(Error::InterpolationError( - "xfxq2_cheby_batch only supports LogChebyshev interpolator".to_string(), + "xfxq2_cheby_batch only supports `LogChebyshev` interpolator".to_string(), )); } let subgrid = &self.knot_array.subgrids[subgrid_idx]; - let use_log = true; // LogChebyshev is always log - let log_points: Vec> = points .iter() - .map(|p| { - p.iter() - .map(|&v| if use_log { v.ln() } else { v }) - .collect() - }) + .map(|p| p.iter().map(|&v| v.ln()).collect()) .collect(); - let results = match subgrid.interpolation_config() { - super::interpolator::InterpolationConfig::TwoD => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<2>::default(); - let grid_slice = subgrid.grid_slice(pid_idx).to_owned(); - let data = InterpData2D::new( - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - grid_slice, - ) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - batch_interpolator - .init(&data) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - - let points_2d: Vec<[f64; 2]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 2D")) - .collect(); - - batch_interpolator - .interpolate(&data, &points_2d) - .map_err(|e| Error::InterpolationError(e.to_string())) - } - super::interpolator::InterpolationConfig::ThreeDNucleons => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); - let grid_data = subgrid.grid.slice(s![.., 0, pid_idx, 0, .., ..]).to_owned(); - let reshaped_data = grid_data - .into_shape_with_order(( - subgrid.nucleons.len(), - subgrid.xs.len(), - subgrid.q2s.len(), - )) - .expect("Failed to reshape 3D data"); - let data = InterpData3D::new( - subgrid.nucleons.mapv(f64::ln), - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - reshaped_data, - ) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - batch_interpolator - .init(&data) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - let points_3d: Vec<[f64; 3]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 3D")) - .collect(); - batch_interpolator - .interpolate(&data, &points_3d) - .map_err(|e| Error::InterpolationError(e.to_string())) - } - super::interpolator::InterpolationConfig::ThreeDAlphas => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); - let grid_data = subgrid.grid.slice(s![0, .., pid_idx, 0, .., ..]).to_owned(); - let reshaped_data = grid_data - .into_shape_with_order(( - subgrid.alphas.len(), - subgrid.xs.len(), - subgrid.q2s.len(), - )) - .expect("Failed to reshape 3D data"); - let data = InterpData3D::new( - subgrid.alphas.mapv(f64::ln), - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - reshaped_data, - ) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - batch_interpolator - .init(&data) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - let points_3d: Vec<[f64; 3]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 3D")) - .collect(); - batch_interpolator - .interpolate(&data, &points_3d) - .map_err(|e| Error::InterpolationError(e.to_string())) - } - super::interpolator::InterpolationConfig::ThreeDKt => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); - let grid_data = subgrid.grid.slice(s![0, 0, pid_idx, .., .., ..]).to_owned(); - let reshaped_data = grid_data - .into_shape_with_order((subgrid.kts.len(), subgrid.xs.len(), subgrid.q2s.len())) - .expect("Failed to reshape 3D data"); - let data = InterpData3D::new( - subgrid.kts.mapv(f64::ln), - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - reshaped_data, - ) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - batch_interpolator - .init(&data) - .map_err(|e| Error::InterpolationError(e.to_string()))?; - let points_3d: Vec<[f64; 3]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 3D")) - .collect(); - batch_interpolator - .interpolate(&data, &points_3d) - .map_err(|e| Error::InterpolationError(e.to_string())) - } - _ => Err(Error::InterpolationError( - "Unsupported dimension for xfxq2_cheby_batch".to_string(), - )), - }?; + let results = + super::interpolator::chebyshev_batch_interpolate(subgrid, pid_idx, log_points) + .map_err(Error::InterpolationError)?; Ok(results .into_iter() diff --git a/neopdf/src/interpolator.rs b/neopdf/src/interpolator.rs index 411ac764..a080e963 100644 --- a/neopdf/src/interpolator.rs +++ b/neopdf/src/interpolator.rs @@ -12,6 +12,7 @@ //! The [`SubGrid`] struct is defined in `subgrid.rs`. use ndarray::{s, OwnedRepr}; +use ninterp::data::{InterpData2D, InterpData3D}; use ninterp::error::InterpolateError; use ninterp::interpolator::{ Extrapolate, Interp2D, Interp2DOwned, Interp3D, Interp3DOwned, InterpND, InterpNDOwned, @@ -23,7 +24,7 @@ use ninterp::strategy::Linear; use super::metadata::InterpolatorType; use super::strategy::{ BilinearInterpolation, LogBicubicInterpolation, LogBilinearInterpolation, - LogChebyshevInterpolation, LogTricubicInterpolation, + LogChebyshevBatchInterpolation, LogChebyshevInterpolation, LogTricubicInterpolation, }; use super::subgrid::SubGrid; @@ -457,6 +458,120 @@ impl InterpolatorFactory { } } +pub fn chebyshev_batch_interpolate( + subgrid: &SubGrid, + pid_idx: usize, + log_points: Vec>, +) -> Result, String> { + match subgrid.interpolation_config() { + InterpolationConfig::TwoD => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<2>::default(); + let grid_slice = subgrid.grid_slice(pid_idx).to_owned(); + let data = InterpData2D::new( + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + grid_slice, + ) + .map_err(|e| e.to_string())?; + + batch_interpolator.init(&data).map_err(|e| e.to_string())?; + + let points_2d: Vec<[f64; 2]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 2D")) + .collect(); + + batch_interpolator + .interpolate(&data, &points_2d) + .map_err(|e| e.to_string()) + } + + InterpolationConfig::ThreeDNucleons => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![.., 0, pid_idx, 0, .., ..]).to_owned(); + let reshaped_data = grid_data + .into_shape_with_order(( + subgrid.nucleons.len(), + subgrid.xs.len(), + subgrid.q2s.len(), + )) + .expect("Failed to reshape 3D data"); + + let data = InterpData3D::new( + subgrid.nucleons.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| e.to_string())?; + + batch_interpolator.init(&data).map_err(|e| e.to_string())?; + + let points_3d: Vec<[f64; 3]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + + batch_interpolator + .interpolate(&data, &points_3d) + .map_err(|e| e.to_string()) + } + InterpolationConfig::ThreeDAlphas => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![0, .., pid_idx, 0, .., ..]).to_owned(); + let reshaped_data = grid_data + .into_shape_with_order((subgrid.alphas.len(), subgrid.xs.len(), subgrid.q2s.len())) + .expect("Failed to reshape 3D data"); + + let data = InterpData3D::new( + subgrid.alphas.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| e.to_string())?; + + batch_interpolator.init(&data).map_err(|e| e.to_string())?; + + let points_3d: Vec<[f64; 3]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + + batch_interpolator + .interpolate(&data, &points_3d) + .map_err(|e| e.to_string()) + } + InterpolationConfig::ThreeDKt => { + let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![0, 0, pid_idx, .., .., ..]).to_owned(); + let reshaped_data = grid_data + .into_shape_with_order((subgrid.kts.len(), subgrid.xs.len(), subgrid.q2s.len())) + .expect("Failed to reshape 3D data"); + + let data = InterpData3D::new( + subgrid.kts.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| e.to_string())?; + + batch_interpolator.init(&data).map_err(|e| e.to_string())?; + + let points_3d: Vec<[f64; 3]> = log_points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + + batch_interpolator + .interpolate(&data, &points_3d) + .map_err(|e| e.to_string()) + } + _ => Err("Unsupported dimension for xfxq2_cheby_batch".to_string()), + } +} + #[cfg(test)] mod tests { use super::*; From 75029e0d4662b3287a4e0e58361400b41fc3ad5c Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Tue, 2 Sep 2025 17:27:56 +0200 Subject: [PATCH 04/21] Remove boundary check to allow extrapolations --- neopdf/src/strategy.rs | 37 +++---------------------------------- 1 file changed, 3 insertions(+), 34 deletions(-) diff --git a/neopdf/src/strategy.rs b/neopdf/src/strategy.rs index 9474f362..5127a3aa 100644 --- a/neopdf/src/strategy.rs +++ b/neopdf/src/strategy.rs @@ -1326,14 +1326,7 @@ impl LogChebyshevBatchInterpolation<1> { let mut t_x_vals = Vec::with_capacity(points.len()); for &[x] in points { - if x < x_min || x > x_max { - return Err(InterpolateError::Other(format!( - "Input point x={} is outside the grid range [{}, {}]", - x, x_min, x_max - ))); - } - let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; - t_x_vals.push(t_x); + t_x_vals.push(2.0 * (x - x_min) / (x_max - x_min) - 1.0); } let c_x = Self::barycentric_coefficients(&t_x_vals, &self.t_coords[0], &self.weights[0]); @@ -1384,27 +1377,9 @@ impl LogChebyshevBatchInterpolation<2> { let mut t_x_vals = Vec::with_capacity(points.len()); let mut t_y_vals = Vec::with_capacity(points.len()); - - // TODO: Remove this for extrapolation -> speed gain for &[x, y] in points { - if x < x_min || x > x_max { - return Err(InterpolateError::Other(format!( - "Input point x={} is outside the grid range [{}, {}]", - x, x_min, x_max - ))); - } - if y < y_min || y > y_max { - return Err(InterpolateError::Other(format!( - "Input point y={} is outside the grid range [{}, {}]", - y, y_min, y_max - ))); - } - - let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; - let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; - - t_x_vals.push(t_x); - t_y_vals.push(t_y); + t_x_vals.push(2.0 * (x - x_min) / (x_max - x_min) - 1.0); + t_y_vals.push(2.0 * (y - y_min) / (y_max - y_min) - 1.0); } let c_x = Self::barycentric_coefficients(&t_x_vals, &self.t_coords[0], &self.weights[0]); @@ -1461,13 +1436,7 @@ impl LogChebyshevBatchInterpolation<3> { let mut t_y_vals = Vec::with_capacity(points.len()); let mut t_z_vals = Vec::with_capacity(points.len()); - // TODO: Remove this for extrapolation -> speed gain for &[x, y, z] in points { - if x < x_min || x > x_max || y < y_min || y > y_max || z < z_min || z > z_max { - return Err(InterpolateError::Other( - "Point outside grid range".to_string(), - )); - } t_x_vals.push(2.0 * (x - x_min) / (x_max - x_min) - 1.0); t_y_vals.push(2.0 * (y - y_min) / (y_max - y_min) - 1.0); t_z_vals.push(2.0 * (z - z_min) / (z_max - z_min) - 1.0); From e1c9ff9e85c0f93ec4a0514511ee5275d02839f8 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Wed, 3 Sep 2025 00:05:42 +0200 Subject: [PATCH 05/21] Add Unit Tests to `LogChebyshevBatchInterpolation` --- neopdf/src/strategy.rs | 89 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/neopdf/src/strategy.rs b/neopdf/src/strategy.rs index 5127a3aa..67beb733 100644 --- a/neopdf/src/strategy.rs +++ b/neopdf/src/strategy.rs @@ -1939,4 +1939,93 @@ mod tests { assert_close(result, expected, EPSILON); } + + #[test] + fn test_log_chebyshev_batch_interpolation_1d() { + let n = 21; + let x_min: f64 = 0.1; + let x_max: f64 = 10.0; + let x_coords = create_cheby_grid(n, x_min, x_max); + + let f_values: Vec = x_coords.iter().map(|&x| x.ln()).collect(); + let data = create_test_data_1d(x_coords.iter().map(|v| v.ln()).collect(), f_values); + let mut cheby = LogChebyshevBatchInterpolation::<1>::default(); + cheby.init(&data).unwrap(); + + let test_points = [[2.5f64.ln()], [5.0f64.ln()], [7.5f64.ln()]]; + let expected: Vec = test_points.iter().map(|p| p[0]).collect(); + let results = cheby.interpolate(&data, &test_points).unwrap(); + + for (res, exp) in results.iter().zip(expected.iter()) { + assert_close(*res, *exp, EPSILON); + } + } + + #[test] + fn test_log_chebyshev_batch_interpolation_2d() { + let n = 11; + let x_coords = create_cheby_grid(n, 0.1, 10.0); + let y_coords = create_cheby_grid(n, 0.1, 10.0); + + let f_values: Vec = x_coords + .iter() + .flat_map(|&x| y_coords.iter().map(move |&y| x.ln() + y.ln())) + .collect(); + + let data = create_test_data_2d( + x_coords.iter().map(|v| v.ln()).collect(), + y_coords.iter().map(|v| v.ln()).collect(), + f_values, + ); + let mut cheby = LogChebyshevBatchInterpolation::<2>::default(); + cheby.init(&data).unwrap(); + + let test_points = [ + [2.5f64.ln(), 3.5f64.ln()], + [5.0f64.ln(), 6.0f64.ln()], + [7.5f64.ln(), 8.5f64.ln()], + ]; + let expected: Vec = test_points.iter().map(|p| p[0] + p[1]).collect(); + let results = cheby.interpolate(&data, &test_points).unwrap(); + + for (res, exp) in results.iter().zip(expected.iter()) { + assert_close(*res, *exp, EPSILON); + } + } + + #[test] + fn test_log_chebyshev_batch_interpolation_3d() { + let n = 7; + let x_coords = create_cheby_grid(n, 0.1, 10.0); + let y_coords = create_cheby_grid(n, 0.1, 10.0); + let z_coords = create_cheby_grid(n, 0.1, 10.0); + + let f_values: Vec = x_coords + .iter() + .cartesian_product(y_coords.iter()) + .cartesian_product(z_coords.iter()) + .map(|((&x, &y), &z)| x.ln() + y.ln() + z.ln()) + .collect(); + + let data = create_test_data_3d( + x_coords.iter().map(|v| v.ln()).collect(), + y_coords.iter().map(|v| v.ln()).collect(), + z_coords.iter().map(|v| v.ln()).collect(), + f_values, + ); + let mut cheby = LogChebyshevBatchInterpolation::<3>::default(); + cheby.init(&data).unwrap(); + + let test_points = [ + [2.5f64.ln(), 3.5f64.ln(), 4.5f64.ln()], + [5.0f64.ln(), 6.0f64.ln(), 7.0f64.ln()], + [7.5f64.ln(), 8.5f64.ln(), 9.5f64.ln()], + ]; + let expected: Vec = test_points.iter().map(|p| p[0] + p[1] + p[2]).collect(); + let results = cheby.interpolate(&data, &test_points).unwrap(); + + for (res, exp) in results.iter().zip(expected.iter()) { + assert_close(*res, *exp, EPSILON); + } + } } From 0c9dfb3b7234b57441bfae6be6cef762fc9413f9 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Wed, 3 Sep 2025 08:37:31 +0200 Subject: [PATCH 06/21] Move logics for Batch Chebyshev into `interpolator.rs` --- neopdf/src/gridpdf.rs | 11 +- neopdf/src/interpolator.rs | 228 +++++++++++++++++++------------------ 2 files changed, 125 insertions(+), 114 deletions(-) diff --git a/neopdf/src/gridpdf.rs b/neopdf/src/gridpdf.rs index 4e8e6193..b9f9de95 100644 --- a/neopdf/src/gridpdf.rs +++ b/neopdf/src/gridpdf.rs @@ -365,7 +365,7 @@ impl GridPDF { if !matches!(self.info.interpolator_type, InterpolatorType::LogChebyshev) { return Err(Error::InterpolationError( - "xfxq2_cheby_batch only supports `LogChebyshev` interpolator".to_string(), + "xfxq2_cheby_batch only supports LogChebyshev interpolator".to_string(), )); } @@ -375,9 +375,12 @@ impl GridPDF { .map(|p| p.iter().map(|&v| v.ln()).collect()) .collect(); - let results = - super::interpolator::chebyshev_batch_interpolate(subgrid, pid_idx, log_points) - .map_err(Error::InterpolationError)?; + let batch_interpolator = InterpolatorFactory::create_batch_interpolator(subgrid, pid_idx) + .map_err(Error::InterpolationError)?; + + let results = batch_interpolator + .interpolate(log_points) + .map_err(|e| Error::InterpolationError(e.to_string()))?; Ok(results .into_iter() diff --git a/neopdf/src/interpolator.rs b/neopdf/src/interpolator.rs index a080e963..43e886bc 100644 --- a/neopdf/src/interpolator.rs +++ b/neopdf/src/interpolator.rs @@ -116,6 +116,40 @@ where } } +/// An enum to dispatch batch interpolation to the correct Chebyshev interpolator. +pub enum BatchInterpolator { + Chebyshev2D( + LogChebyshevBatchInterpolation<2>, + InterpData2D>, + ), + Chebyshev3D( + LogChebyshevBatchInterpolation<3>, + InterpData3D>, + ), +} + +impl BatchInterpolator { + /// Interpolates a batch of points. + pub fn interpolate(&self, points: Vec>) -> Result, InterpolateError> { + match self { + BatchInterpolator::Chebyshev2D(strategy, data) => { + let points_2d: Vec<[f64; 2]> = points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 2D")) + .collect(); + strategy.interpolate(data, &points_2d) + } + BatchInterpolator::Chebyshev3D(strategy, data) => { + let points_3d: Vec<[f64; 3]> = points + .into_iter() + .map(|p| p.try_into().expect("Invalid point dimension for 3D")) + .collect(); + strategy.interpolate(data, &points_3d) + } + } + } +} + /// Factory for creating dynamic interpolators based on interpolation type and grid dimensions. pub struct InterpolatorFactory; @@ -456,119 +490,93 @@ impl InterpolatorFactory { _ => panic!("Unsupported 5D interpolator: {:?}", interp_type), } } -} -pub fn chebyshev_batch_interpolate( - subgrid: &SubGrid, - pid_idx: usize, - log_points: Vec>, -) -> Result, String> { - match subgrid.interpolation_config() { - InterpolationConfig::TwoD => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<2>::default(); - let grid_slice = subgrid.grid_slice(pid_idx).to_owned(); - let data = InterpData2D::new( - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - grid_slice, - ) - .map_err(|e| e.to_string())?; - - batch_interpolator.init(&data).map_err(|e| e.to_string())?; - - let points_2d: Vec<[f64; 2]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 2D")) - .collect(); - - batch_interpolator - .interpolate(&data, &points_2d) - .map_err(|e| e.to_string()) - } + pub fn create_batch_interpolator( + subgrid: &SubGrid, + pid_idx: usize, + ) -> Result { + match subgrid.interpolation_config() { + InterpolationConfig::TwoD => { + let mut strategy = LogChebyshevBatchInterpolation::<2>::default(); + let grid_slice = subgrid.grid_slice(pid_idx).to_owned(); - InterpolationConfig::ThreeDNucleons => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); - let grid_data = subgrid.grid.slice(s![.., 0, pid_idx, 0, .., ..]).to_owned(); - let reshaped_data = grid_data - .into_shape_with_order(( - subgrid.nucleons.len(), - subgrid.xs.len(), - subgrid.q2s.len(), - )) - .expect("Failed to reshape 3D data"); - - let data = InterpData3D::new( - subgrid.nucleons.mapv(f64::ln), - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - reshaped_data, - ) - .map_err(|e| e.to_string())?; - - batch_interpolator.init(&data).map_err(|e| e.to_string())?; - - let points_3d: Vec<[f64; 3]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 3D")) - .collect(); - - batch_interpolator - .interpolate(&data, &points_3d) - .map_err(|e| e.to_string()) - } - InterpolationConfig::ThreeDAlphas => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); - let grid_data = subgrid.grid.slice(s![0, .., pid_idx, 0, .., ..]).to_owned(); - let reshaped_data = grid_data - .into_shape_with_order((subgrid.alphas.len(), subgrid.xs.len(), subgrid.q2s.len())) - .expect("Failed to reshape 3D data"); - - let data = InterpData3D::new( - subgrid.alphas.mapv(f64::ln), - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - reshaped_data, - ) - .map_err(|e| e.to_string())?; - - batch_interpolator.init(&data).map_err(|e| e.to_string())?; - - let points_3d: Vec<[f64; 3]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 3D")) - .collect(); - - batch_interpolator - .interpolate(&data, &points_3d) - .map_err(|e| e.to_string()) - } - InterpolationConfig::ThreeDKt => { - let mut batch_interpolator = LogChebyshevBatchInterpolation::<3>::default(); - let grid_data = subgrid.grid.slice(s![0, 0, pid_idx, .., .., ..]).to_owned(); - let reshaped_data = grid_data - .into_shape_with_order((subgrid.kts.len(), subgrid.xs.len(), subgrid.q2s.len())) - .expect("Failed to reshape 3D data"); - - let data = InterpData3D::new( - subgrid.kts.mapv(f64::ln), - subgrid.xs.mapv(f64::ln), - subgrid.q2s.mapv(f64::ln), - reshaped_data, - ) - .map_err(|e| e.to_string())?; - - batch_interpolator.init(&data).map_err(|e| e.to_string())?; - - let points_3d: Vec<[f64; 3]> = log_points - .into_iter() - .map(|p| p.try_into().expect("Invalid point dimension for 3D")) - .collect(); - - batch_interpolator - .interpolate(&data, &points_3d) - .map_err(|e| e.to_string()) + let data = InterpData2D::new( + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + grid_slice, + ) + .map_err(|e| e.to_string())?; + strategy.init(&data).map_err(|e| e.to_string())?; + + Ok(BatchInterpolator::Chebyshev2D(strategy, data)) + } + InterpolationConfig::ThreeDNucleons => { + let mut strategy = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![.., 0, pid_idx, 0, .., ..]).to_owned(); + + let reshaped_data = grid_data + .into_shape_with_order(( + subgrid.nucleons.len(), + subgrid.xs.len(), + subgrid.q2s.len(), + )) + .expect("Failed to reshape 3D data"); + + let data = InterpData3D::new( + subgrid.nucleons.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| e.to_string())?; + strategy.init(&data).map_err(|e| e.to_string())?; + + Ok(BatchInterpolator::Chebyshev3D(strategy, data)) + } + InterpolationConfig::ThreeDAlphas => { + let mut strategy = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![0, .., pid_idx, 0, .., ..]).to_owned(); + + let reshaped_data = grid_data + .into_shape_with_order(( + subgrid.alphas.len(), + subgrid.xs.len(), + subgrid.q2s.len(), + )) + .expect("Failed to reshape 3D data"); + + let data = InterpData3D::new( + subgrid.alphas.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| e.to_string())?; + strategy.init(&data).map_err(|e| e.to_string())?; + + Ok(BatchInterpolator::Chebyshev3D(strategy, data)) + } + InterpolationConfig::ThreeDKt => { + let mut strategy = LogChebyshevBatchInterpolation::<3>::default(); + let grid_data = subgrid.grid.slice(s![0, 0, pid_idx, .., .., ..]).to_owned(); + + let reshaped_data = grid_data + .into_shape_with_order((subgrid.kts.len(), subgrid.xs.len(), subgrid.q2s.len())) + .expect("Failed to reshape 3D data"); + + let data = InterpData3D::new( + subgrid.kts.mapv(f64::ln), + subgrid.xs.mapv(f64::ln), + subgrid.q2s.mapv(f64::ln), + reshaped_data, + ) + .map_err(|e| e.to_string())?; + strategy.init(&data).map_err(|e| e.to_string())?; + + Ok(BatchInterpolator::Chebyshev3D(strategy, data)) + } + _ => Err("Unsupported dimension for batch interpolation".to_string()), } - _ => Err("Unsupported dimension for xfxq2_cheby_batch".to_string()), } } From 357cbf3ce7463ee8cec8600770bb78cfde4b3f65 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Wed, 3 Sep 2025 08:53:43 +0200 Subject: [PATCH 07/21] Propagate changes into the Fortran API --- neopdf_fapi/neopdf.f90 | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/neopdf_fapi/neopdf.f90 b/neopdf_fapi/neopdf.f90 index 02d2e786..c53a1292 100644 --- a/neopdf_fapi/neopdf.f90 +++ b/neopdf_fapi/neopdf.f90 @@ -222,6 +222,16 @@ function c_neopdf_pdf_xfxq2_nd(pdf, id, params, num_params) bind(c, name="neopdf real (c_double) :: c_neopdf_pdf_xfxq2_nd end function + subroutine c_neopdf_pdf_xfxq2_cheby_batch(pdf, flavor_id, points, num_points, num_dims, results) bind(c, name="neopdf_pdf_xfxq2_cheby_batch") + use iso_c_binding + type(c_ptr), value :: pdf + integer(c_int), value :: flavor_id + real(c_double), intent(in) :: points(*) + integer(c_size_t), value :: num_points + integer(c_size_t), value :: num_dims + real(c_double), intent(out) :: results(*) + end subroutine + subroutine c_neopdf_pdf_set_force_positive(pdf, option) bind(c, name="neopdf_pdf_set_force_positive") use iso_c_binding type (c_ptr), value :: pdf @@ -510,6 +520,20 @@ function neopdf_pdf_xfxq2_nd(pdf, id, params) result(res) res = c_neopdf_pdf_xfxq2_nd(pdf%ptr, id, params, int(size(params), c_size_t)) end function + subroutine neopdf_pdf_xfxq2_cheby_batch(pdf, flavor_id, points, results) + implicit none + type(neopdf_pdf), intent(in) :: pdf + integer, intent(in) :: flavor_id + real(dp), intent(in) :: points(:,:) + real(dp), intent(out) :: results(:) + integer(c_size_t) :: num_dims, num_points + + num_dims = size(points, 1) + num_points = size(points, 2) + + call c_neopdf_pdf_xfxq2_cheby_batch(pdf%ptr, flavor_id, points, num_points, num_dims, results) + end subroutine + subroutine neopdf_pdf_set_force_positive(pdf, option) implicit none type (neopdf_pdf), intent(in) :: pdf From d4ad1c69ac408054f5ce68468a30623cfc384087 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Wed, 3 Sep 2025 09:11:55 +0200 Subject: [PATCH 08/21] Propagate the changes into the Python API --- neopdf_pyapi/src/pdf.rs | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/neopdf_pyapi/src/pdf.rs b/neopdf_pyapi/src/pdf.rs index 59802332..17ad0f30 100644 --- a/neopdf_pyapi/src/pdf.rs +++ b/neopdf_pyapi/src/pdf.rs @@ -403,6 +403,33 @@ impl PyPDF { self.pdf.xfxq2(id, ¶ms) } + /// Interpolates the PDF value (xf) for a list containg a set of parameters. + /// + /// Parameters + /// ---------- + /// id : int + /// The flavor ID (e.g., 21 for gluon, 1 for d-quark). + /// params: list[list[float]] + /// A list containing the list of points. Each element in the list + /// is in turn a list containing the parameters that the grids depends + /// on. If the PDF grid only contains `x` and `Q2` dependence then its + /// value is `[x, q2]`; if it contains either the `A` and `alpha_s` + /// dependence, then its value is `[A, x, q2]` or `[alpha_s, x, q2]` + /// respectively; if it contains both, then `[A, alpha_s, x, q2]`. + /// + /// Returns + /// ------- + /// float + /// The interpolated PDF value. Returns 0.0 if extrapolation is + /// attempted and not allowed. + #[must_use] + #[pyo3(name = "xfxQ2_Chebyshev_batch")] + #[allow(clippy::needless_pass_by_value)] + pub fn xfxq2_cheby_batch(&self, id: i32, params: Vec>) -> Vec { + let slices: Vec<&[f64]> = params.iter().map(Vec::as_slice).collect(); + self.pdf.xfxq2_cheby_batch(id, &slices) + } + /// Interpolates the PDF value (xf) for lists of flavors, x-values, /// and Q2-values. /// From e34116fa79eea676c4ab6b66d7d785e83e2d79b2 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Tue, 9 Sep 2025 12:14:37 +0200 Subject: [PATCH 09/21] Make sequential Chebyshev ~20 times faster --- neopdf/src/strategy.rs | 381 +++++++++++++++++++++++++++------------- neopdf_cli/tests/pdf.rs | 2 +- 2 files changed, 258 insertions(+), 125 deletions(-) diff --git a/neopdf/src/strategy.rs b/neopdf/src/strategy.rs index 67beb733..e2c3ed08 100644 --- a/neopdf/src/strategy.rs +++ b/neopdf/src/strategy.rs @@ -21,6 +21,9 @@ use ninterp::strategy::traits::{Strategy1D, Strategy2D, Strategy3D}; use serde::{Deserialize, Serialize}; use std::f64::consts::PI; +// Tolerance for coincidence detection +const COINCIDENCE_TOL: f64 = 1e-15; + use super::utils; /// Implements bilinear interpolation for 2D data. @@ -881,12 +884,17 @@ where /// number of grid points in each dimension. /// - **Grid Requirement**: For optimal stability and to avoid Runge's phenomenon, the grid /// points should correspond to the roots or extrema of Chebyshev polynomials. + #[derive(Debug, Clone)] pub struct LogChebyshevInterpolation { // Pre-computed weights for the barycentric formula for each dimension. weights: [Vec; DIM], // Grid points in the t-domain [-1, 1] for each dimension. t_coords: [Vec; DIM], + // Pre-computed bounds for fast bounds checking and coordinate transformation. + bounds: [(f64, f64); DIM], + // Pre-computed scaling factors for coordinate transformation. + scale_factors: [f64; DIM], } impl Default for LogChebyshevInterpolation { @@ -894,6 +902,8 @@ impl Default for LogChebyshevInterpolation { Self { weights: std::array::from_fn(|_| Vec::new()), t_coords: std::array::from_fn(|_| Vec::new()), + bounds: [(0.0, 0.0); DIM], + scale_factors: [1.0; DIM], } } } @@ -904,9 +914,11 @@ impl Serialize for LogChebyshevInterpolation { S: serde::Serializer, { use serde::ser::SerializeStruct; - let mut state = serializer.serialize_struct("LogChebyshevInterpolation", 2)?; + let mut state = serializer.serialize_struct("LogChebyshevInterpolation", 4)?; state.serialize_field("weights", &self.weights.as_slice())?; state.serialize_field("t_coords", &self.t_coords.as_slice())?; + state.serialize_field("bounds", &self.bounds.as_slice())?; + state.serialize_field("scale_factors", &self.scale_factors.as_slice())?; state.end() } } @@ -920,17 +932,34 @@ impl<'de, const DIM: usize> Deserialize<'de> for LogChebyshevInterpolation struct Helper { weights: Vec>, t_coords: Vec>, + bounds: Vec<(f64, f64)>, + scale_factors: Vec, } let helper = Helper::deserialize(deserializer)?; + let weights = helper.weights.try_into().map_err(|v: Vec>| { serde::de::Error::invalid_length(v.len(), &"an array of the correct length") })?; + let t_coords = helper.t_coords.try_into().map_err(|v: Vec>| { serde::de::Error::invalid_length(v.len(), &"an array of the correct length") })?; - Ok(Self { weights, t_coords }) + let bounds = helper.bounds.try_into().map_err(|v: Vec<(f64, f64)>| { + serde::de::Error::invalid_length(v.len(), &"an array of the correct length") + })?; + + let scale_factors = helper.scale_factors.try_into().map_err(|v: Vec| { + serde::de::Error::invalid_length(v.len(), &"an array of the correct length") + })?; + + Ok(Self { + weights, + t_coords, + bounds, + scale_factors, + }) } } @@ -952,96 +981,242 @@ impl LogChebyshevInterpolation { weights } - /// Computes normalized barycentric coefficients for interpolation - /// Returns a vector of coefficients that sum to 1 - fn barycentric_coefficients(t: f64, t_coords: &[f64], weights: &[f64]) -> Vec { - let mut coeffs = vec![0.0; t_coords.len()]; - + fn interpolate_1d(&self, t: f64, t_coords: &[f64], weights: &[f64], f_value: F) -> f64 + where + F: Fn(usize) -> f64, + { + // Check for coincidence first for (j, &t_j) in t_coords.iter().enumerate() { - if (t - t_j).abs() < 1e-15 { - // t is exactly at grid point j - return unit vector - coeffs[j] = 1.0; - return coeffs; + if (t - t_j).abs() < COINCIDENCE_TOL { + return f_value(j); } } - let mut terms = Vec::with_capacity(t_coords.len()); + let mut numer = 0.0; + let mut denom = 0.0; + for (j, &t_j) in t_coords.iter().enumerate() { - terms.push(weights[j] / (t - t_j)); + let term = weights[j] / (t - t_j); + numer += term * f_value(j); + denom += term; + } + + numer / denom + } + + #[allow(clippy::too_many_arguments)] + fn interpolate_2d( + &self, + t1: f64, + t2: f64, + t1_coords: &[f64], + t2_coords: &[f64], + weights1: &[f64], + weights2: &[f64], + f_value: F, + ) -> f64 + where + F: Fn(usize, usize) -> f64, + { + let mut coincident_1 = None; + let mut coincident_2 = None; + + for (j, &t1_j) in t1_coords.iter().enumerate() { + if (t1 - t1_j).abs() < COINCIDENCE_TOL { + coincident_1 = Some(j); + break; + } } - let sum: f64 = terms.iter().sum(); - for (j, &term) in terms.iter().enumerate() { - coeffs[j] = term / sum; + for (k, &t2_k) in t2_coords.iter().enumerate() { + if (t2 - t2_k).abs() < COINCIDENCE_TOL { + coincident_2 = Some(k); + break; + } } - coeffs - } + match (coincident_1, coincident_2) { + (Some(j), Some(k)) => f_value(j, k), + (Some(j), None) => self.interpolate_1d(t2, t2_coords, weights2, |k| f_value(j, k)), + (None, Some(k)) => self.interpolate_1d(t1, t1_coords, weights1, |j| f_value(j, k)), + (None, None) => { + let mut numer = 0.0; + let mut denom = 0.0; - /// Legacy barycentric interpolation method (kept for compatibility) - fn barycentric_interpolate(t: f64, t_coords: &[f64], f_values: &[f64], weights: &[f64]) -> f64 { - let mut numer = 0.0; - let mut denom = 0.0; + for (j, &t1_j) in t1_coords.iter().enumerate() { + let term1 = weights1[j] / (t1 - t1_j); - for (j, &t_j) in t_coords.iter().enumerate() { - if (t - t_j).abs() < 1e-15 { - return f_values[j]; + for (k, &t2_k) in t2_coords.iter().enumerate() { + let term2 = weights2[k] / (t2 - t2_k); + let weight = term1 * term2; + + numer += weight * f_value(j, k); + denom += weight; + } + } + + numer / denom } + } + } - let term = weights[j] / (t - t_j); - numer += term * f_values[j]; - denom += term; + #[allow(clippy::too_many_arguments)] + fn interpolate_3d( + &self, + t1: f64, + t2: f64, + t3: f64, + t1_coords: &[f64], + t2_coords: &[f64], + t3_coords: &[f64], + weights1: &[f64], + weights2: &[f64], + weights3: &[f64], + f_value: F, + ) -> f64 + where + F: Fn(usize, usize, usize) -> f64, + { + let mut coincident_dims = [None; 3]; + let t_vals = [t1, t2, t3]; + let t_coords_all = [t1_coords, t2_coords, t3_coords]; + + for (dim, (&t_val, t_coords)) in t_vals.iter().zip(t_coords_all.iter()).enumerate() { + for (idx, &t_coord) in t_coords.iter().enumerate() { + if (t_val - t_coord).abs() < COINCIDENCE_TOL { + coincident_dims[dim] = Some(idx); + break; + } + } } - numer / denom + match coincident_dims { + [Some(i), Some(j), Some(k)] => f_value(i, j, k), + [Some(i), Some(j), None] => { + self.interpolate_1d(t3, t3_coords, weights3, |k| f_value(i, j, k)) + } + [Some(i), None, Some(k)] => { + self.interpolate_1d(t2, t2_coords, weights2, |j| f_value(i, j, k)) + } + [None, Some(j), Some(k)] => { + self.interpolate_1d(t1, t1_coords, weights1, |i| f_value(i, j, k)) + } + [Some(i), None, None] => { + self.interpolate_2d(t2, t3, t2_coords, t3_coords, weights2, weights3, |j, k| { + f_value(i, j, k) + }) + } + [None, Some(j), None] => { + self.interpolate_2d(t1, t3, t1_coords, t3_coords, weights1, weights3, |i, k| { + f_value(i, j, k) + }) + } + [None, None, Some(k)] => { + self.interpolate_2d(t1, t2, t1_coords, t2_coords, weights1, weights2, |i, j| { + f_value(i, j, k) + }) + } + [None, None, None] => { + let mut result = 0.0; + let mut total_weight = 0.0; + + for (i, &t1_i) in t1_coords.iter().enumerate() { + let w1 = weights1[i] / (t1 - t1_i); + + for (j, &t2_j) in t2_coords.iter().enumerate() { + let w2 = weights2[j] / (t2 - t2_j); + let w12 = w1 * w2; + + for (k, &t3_k) in t3_coords.iter().enumerate() { + let w3 = weights3[k] / (t3 - t3_k); + let weight = w12 * w3; + + result += weight * f_value(i, j, k); + total_weight += weight; + } + } + } + + result / total_weight + } + } } -} -impl Strategy1D for LogChebyshevInterpolation<1> -where - D: Data + RawDataClone + Clone, -{ - fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> { - let x_coords = data.grid[0].as_slice().unwrap(); + fn init_dimension(&mut self, dim: usize, x_coords: &[f64]) -> Result<(), ValidateError> { let n = x_coords.len(); if n < 2 { return Err(ValidateError::Other( - "LogChebyshevInterpolation requires at least 2 grid points.".to_string(), + "LogChebyshevInterpolation requires at least 2 grid points per dimension." + .to_string(), )); } - self.t_coords[0] = (0..n) + self.t_coords[dim] = (0..n) .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) .collect(); + self.weights[dim] = Self::compute_barycentric_weights(n); - self.weights[0] = Self::compute_barycentric_weights(n); + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + + if (x_max - x_min).abs() < 1e-15 { + return Err(ValidateError::Other( + "Grid coordinates have zero range".to_string(), + )); + } + + self.bounds[dim] = (x_min, x_max); + self.scale_factors[dim] = 1.0 / (x_max - x_min); Ok(()) } + fn check_bounds_and_transform_generic( + &self, + point: &[f64], + dim_count: usize, + ) -> Result, InterpolateError> { + let mut transformed = vec![0.0; dim_count]; + + for i in 0..dim_count { + let val = point[i]; + let (min, max) = self.bounds[i]; + + if val < min || val > max { + return Err(InterpolateError::Other(format!( + "Input point coordinate {} = {} is outside the grid range [{}, {}]", + i, val, min, max + ))); + } + + transformed[i] = 2.0 * (val - min) * self.scale_factors[i] - 1.0; + } + + Ok(transformed) + } +} + +impl Strategy1D for LogChebyshevInterpolation<1> +where + D: Data + RawDataClone + Clone, +{ + fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> { + let x_coords = data.grid[0].as_slice().unwrap(); + self.init_dimension(0, x_coords) + } + fn interpolate( &self, data: &InterpData1D, point: &[f64; 1], ) -> Result { - let x = point[0]; - let x_coords = data.grid[0].as_slice().unwrap(); - let f_values = data.values.as_slice().unwrap(); - - let x_min = *x_coords.first().unwrap(); - let x_max = *x_coords.last().unwrap(); + let t_coords = self.check_bounds_and_transform_generic(&point[..], 1)?; + let t = t_coords[0]; - if (x_max - x_min).abs() < 1e-15 { - return Ok(f_values[0]); - } - let t = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; + let f_values = data.values.as_slice().unwrap(); + let result = self.interpolate_1d(t, &self.t_coords[0], &self.weights[0], |j| f_values[j]); - Ok(Self::barycentric_interpolate( - t, - &self.t_coords[0], - f_values, - &self.weights[0], - )) + Ok(result) } fn allow_extrapolate(&self) -> bool { @@ -1056,17 +1231,7 @@ where fn init(&mut self, data: &InterpData2D) -> Result<(), ValidateError> { for dim in 0..2 { let x_coords = data.grid[dim].as_slice().unwrap(); - let n = x_coords.len(); - if n < 2 { - return Err(ValidateError::Other( - "LogChebyshevInterpolation requires at least 2 grid points per dimension." - .to_string(), - )); - } - self.t_coords[dim] = (0..n) - .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) - .collect(); - self.weights[dim] = Self::compute_barycentric_weights(n); + self.init_dimension(dim, x_coords)?; } Ok(()) } @@ -1076,27 +1241,18 @@ where data: &InterpData2D, point: &[f64; 2], ) -> Result { - let [x, y] = *point; - let x_coords = data.grid[0].as_slice().unwrap(); - let y_coords = data.grid[1].as_slice().unwrap(); - - let x_min = *x_coords.first().unwrap(); - let x_max = *x_coords.last().unwrap(); - let y_min = *y_coords.first().unwrap(); - let y_max = *y_coords.last().unwrap(); + let t_coords = self.check_bounds_and_transform_generic(&point[..], 2)?; + let [t_x, t_y] = [t_coords[0], t_coords[1]]; - let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; - let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; - - let x_coeffs = Self::barycentric_coefficients(t_x, &self.t_coords[0], &self.weights[0]); - let y_coeffs = Self::barycentric_coefficients(t_y, &self.t_coords[1], &self.weights[1]); - - let mut result = 0.0; - for (i, &x_coeff) in x_coeffs.iter().enumerate() { - for (j, &y_coeff) in y_coeffs.iter().enumerate() { - result += x_coeff * y_coeff * data.values[[i, j]]; - } - } + let result = self.interpolate_2d( + t_x, + t_y, + &self.t_coords[0], + &self.t_coords[1], + &self.weights[0], + &self.weights[1], + |i, j| data.values[[i, j]], + ); Ok(result) } @@ -1113,17 +1269,7 @@ where fn init(&mut self, data: &InterpData3D) -> Result<(), ValidateError> { for dim in 0..3 { let x_coords = data.grid[dim].as_slice().unwrap(); - let n = x_coords.len(); - if n < 2 { - return Err(ValidateError::Other( - "LogChebyshevInterpolation requires at least 2 grid points per dimension." - .to_string(), - )); - } - self.t_coords[dim] = (0..n) - .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) - .collect(); - self.weights[dim] = Self::compute_barycentric_weights(n); + self.init_dimension(dim, x_coords)?; } Ok(()) } @@ -1133,34 +1279,21 @@ where data: &InterpData3D, point: &[f64; 3], ) -> Result { - let [x, y, z] = *point; - let x_coords = data.grid[0].as_slice().unwrap(); - let y_coords = data.grid[1].as_slice().unwrap(); - let z_coords = data.grid[2].as_slice().unwrap(); + let t_coords = self.check_bounds_and_transform_generic(&point[..], 3)?; + let [t_x, t_y, t_z] = [t_coords[0], t_coords[1], t_coords[2]]; - let x_min = *x_coords.first().unwrap(); - let x_max = *x_coords.last().unwrap(); - let y_min = *y_coords.first().unwrap(); - let y_max = *y_coords.last().unwrap(); - let z_min = *z_coords.first().unwrap(); - let z_max = *z_coords.last().unwrap(); - - let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; - let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; - let t_z = 2.0 * (z - z_min) / (z_max - z_min) - 1.0; - - let x_coeffs = Self::barycentric_coefficients(t_x, &self.t_coords[0], &self.weights[0]); - let y_coeffs = Self::barycentric_coefficients(t_y, &self.t_coords[1], &self.weights[1]); - let z_coeffs = Self::barycentric_coefficients(t_z, &self.t_coords[2], &self.weights[2]); - - let mut result = 0.0; - for (i, &x_coeff) in x_coeffs.iter().enumerate() { - for (j, &y_coeff) in y_coeffs.iter().enumerate() { - for (k, &z_coeff) in z_coeffs.iter().enumerate() { - result += x_coeff * y_coeff * z_coeff * data.values[[i, j, k]]; - } - } - } + let result = self.interpolate_3d( + t_x, + t_y, + t_z, + &self.t_coords[0], + &self.t_coords[1], + &self.t_coords[2], + &self.weights[0], + &self.weights[1], + &self.weights[2], + |i, j, k| data.values[[i, j, k]], + ); Ok(result) } diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index ebdde642..f078e689 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -171,5 +171,5 @@ fn xfxq2_kt_tmdlib() { ]) .assert() .success() - .stdout("0.0655544\n"); + .stdout("0.07892252798564643\n"); } From f2f8c50f81bfdb8d84f2a677e672f4765e97f765 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Wed, 10 Sep 2025 00:16:28 +0200 Subject: [PATCH 10/21] Revert "Make sequential Chebyshev ~20 times faster" This reverts commit e34116fa79eea676c4ab6b66d7d785e83e2d79b2. --- neopdf/src/strategy.rs | 381 +++++++++++++--------------------------- neopdf_cli/tests/pdf.rs | 2 +- 2 files changed, 125 insertions(+), 258 deletions(-) diff --git a/neopdf/src/strategy.rs b/neopdf/src/strategy.rs index e2c3ed08..67beb733 100644 --- a/neopdf/src/strategy.rs +++ b/neopdf/src/strategy.rs @@ -21,9 +21,6 @@ use ninterp::strategy::traits::{Strategy1D, Strategy2D, Strategy3D}; use serde::{Deserialize, Serialize}; use std::f64::consts::PI; -// Tolerance for coincidence detection -const COINCIDENCE_TOL: f64 = 1e-15; - use super::utils; /// Implements bilinear interpolation for 2D data. @@ -884,17 +881,12 @@ where /// number of grid points in each dimension. /// - **Grid Requirement**: For optimal stability and to avoid Runge's phenomenon, the grid /// points should correspond to the roots or extrema of Chebyshev polynomials. - #[derive(Debug, Clone)] pub struct LogChebyshevInterpolation { // Pre-computed weights for the barycentric formula for each dimension. weights: [Vec; DIM], // Grid points in the t-domain [-1, 1] for each dimension. t_coords: [Vec; DIM], - // Pre-computed bounds for fast bounds checking and coordinate transformation. - bounds: [(f64, f64); DIM], - // Pre-computed scaling factors for coordinate transformation. - scale_factors: [f64; DIM], } impl Default for LogChebyshevInterpolation { @@ -902,8 +894,6 @@ impl Default for LogChebyshevInterpolation { Self { weights: std::array::from_fn(|_| Vec::new()), t_coords: std::array::from_fn(|_| Vec::new()), - bounds: [(0.0, 0.0); DIM], - scale_factors: [1.0; DIM], } } } @@ -914,11 +904,9 @@ impl Serialize for LogChebyshevInterpolation { S: serde::Serializer, { use serde::ser::SerializeStruct; - let mut state = serializer.serialize_struct("LogChebyshevInterpolation", 4)?; + let mut state = serializer.serialize_struct("LogChebyshevInterpolation", 2)?; state.serialize_field("weights", &self.weights.as_slice())?; state.serialize_field("t_coords", &self.t_coords.as_slice())?; - state.serialize_field("bounds", &self.bounds.as_slice())?; - state.serialize_field("scale_factors", &self.scale_factors.as_slice())?; state.end() } } @@ -932,34 +920,17 @@ impl<'de, const DIM: usize> Deserialize<'de> for LogChebyshevInterpolation struct Helper { weights: Vec>, t_coords: Vec>, - bounds: Vec<(f64, f64)>, - scale_factors: Vec, } let helper = Helper::deserialize(deserializer)?; - let weights = helper.weights.try_into().map_err(|v: Vec>| { serde::de::Error::invalid_length(v.len(), &"an array of the correct length") })?; - let t_coords = helper.t_coords.try_into().map_err(|v: Vec>| { serde::de::Error::invalid_length(v.len(), &"an array of the correct length") })?; - let bounds = helper.bounds.try_into().map_err(|v: Vec<(f64, f64)>| { - serde::de::Error::invalid_length(v.len(), &"an array of the correct length") - })?; - - let scale_factors = helper.scale_factors.try_into().map_err(|v: Vec| { - serde::de::Error::invalid_length(v.len(), &"an array of the correct length") - })?; - - Ok(Self { - weights, - t_coords, - bounds, - scale_factors, - }) + Ok(Self { weights, t_coords }) } } @@ -981,242 +952,96 @@ impl LogChebyshevInterpolation { weights } - fn interpolate_1d(&self, t: f64, t_coords: &[f64], weights: &[f64], f_value: F) -> f64 - where - F: Fn(usize) -> f64, - { - // Check for coincidence first + /// Computes normalized barycentric coefficients for interpolation + /// Returns a vector of coefficients that sum to 1 + fn barycentric_coefficients(t: f64, t_coords: &[f64], weights: &[f64]) -> Vec { + let mut coeffs = vec![0.0; t_coords.len()]; + for (j, &t_j) in t_coords.iter().enumerate() { - if (t - t_j).abs() < COINCIDENCE_TOL { - return f_value(j); + if (t - t_j).abs() < 1e-15 { + // t is exactly at grid point j - return unit vector + coeffs[j] = 1.0; + return coeffs; } } - let mut numer = 0.0; - let mut denom = 0.0; - + let mut terms = Vec::with_capacity(t_coords.len()); for (j, &t_j) in t_coords.iter().enumerate() { - let term = weights[j] / (t - t_j); - numer += term * f_value(j); - denom += term; - } - - numer / denom - } - - #[allow(clippy::too_many_arguments)] - fn interpolate_2d( - &self, - t1: f64, - t2: f64, - t1_coords: &[f64], - t2_coords: &[f64], - weights1: &[f64], - weights2: &[f64], - f_value: F, - ) -> f64 - where - F: Fn(usize, usize) -> f64, - { - let mut coincident_1 = None; - let mut coincident_2 = None; - - for (j, &t1_j) in t1_coords.iter().enumerate() { - if (t1 - t1_j).abs() < COINCIDENCE_TOL { - coincident_1 = Some(j); - break; - } + terms.push(weights[j] / (t - t_j)); } - for (k, &t2_k) in t2_coords.iter().enumerate() { - if (t2 - t2_k).abs() < COINCIDENCE_TOL { - coincident_2 = Some(k); - break; - } + let sum: f64 = terms.iter().sum(); + for (j, &term) in terms.iter().enumerate() { + coeffs[j] = term / sum; } - match (coincident_1, coincident_2) { - (Some(j), Some(k)) => f_value(j, k), - (Some(j), None) => self.interpolate_1d(t2, t2_coords, weights2, |k| f_value(j, k)), - (None, Some(k)) => self.interpolate_1d(t1, t1_coords, weights1, |j| f_value(j, k)), - (None, None) => { - let mut numer = 0.0; - let mut denom = 0.0; - - for (j, &t1_j) in t1_coords.iter().enumerate() { - let term1 = weights1[j] / (t1 - t1_j); - - for (k, &t2_k) in t2_coords.iter().enumerate() { - let term2 = weights2[k] / (t2 - t2_k); - let weight = term1 * term2; - - numer += weight * f_value(j, k); - denom += weight; - } - } - - numer / denom - } - } + coeffs } - #[allow(clippy::too_many_arguments)] - fn interpolate_3d( - &self, - t1: f64, - t2: f64, - t3: f64, - t1_coords: &[f64], - t2_coords: &[f64], - t3_coords: &[f64], - weights1: &[f64], - weights2: &[f64], - weights3: &[f64], - f_value: F, - ) -> f64 - where - F: Fn(usize, usize, usize) -> f64, - { - let mut coincident_dims = [None; 3]; - let t_vals = [t1, t2, t3]; - let t_coords_all = [t1_coords, t2_coords, t3_coords]; - - for (dim, (&t_val, t_coords)) in t_vals.iter().zip(t_coords_all.iter()).enumerate() { - for (idx, &t_coord) in t_coords.iter().enumerate() { - if (t_val - t_coord).abs() < COINCIDENCE_TOL { - coincident_dims[dim] = Some(idx); - break; - } - } - } + /// Legacy barycentric interpolation method (kept for compatibility) + fn barycentric_interpolate(t: f64, t_coords: &[f64], f_values: &[f64], weights: &[f64]) -> f64 { + let mut numer = 0.0; + let mut denom = 0.0; - match coincident_dims { - [Some(i), Some(j), Some(k)] => f_value(i, j, k), - [Some(i), Some(j), None] => { - self.interpolate_1d(t3, t3_coords, weights3, |k| f_value(i, j, k)) - } - [Some(i), None, Some(k)] => { - self.interpolate_1d(t2, t2_coords, weights2, |j| f_value(i, j, k)) - } - [None, Some(j), Some(k)] => { - self.interpolate_1d(t1, t1_coords, weights1, |i| f_value(i, j, k)) - } - [Some(i), None, None] => { - self.interpolate_2d(t2, t3, t2_coords, t3_coords, weights2, weights3, |j, k| { - f_value(i, j, k) - }) - } - [None, Some(j), None] => { - self.interpolate_2d(t1, t3, t1_coords, t3_coords, weights1, weights3, |i, k| { - f_value(i, j, k) - }) - } - [None, None, Some(k)] => { - self.interpolate_2d(t1, t2, t1_coords, t2_coords, weights1, weights2, |i, j| { - f_value(i, j, k) - }) + for (j, &t_j) in t_coords.iter().enumerate() { + if (t - t_j).abs() < 1e-15 { + return f_values[j]; } - [None, None, None] => { - let mut result = 0.0; - let mut total_weight = 0.0; - - for (i, &t1_i) in t1_coords.iter().enumerate() { - let w1 = weights1[i] / (t1 - t1_i); - - for (j, &t2_j) in t2_coords.iter().enumerate() { - let w2 = weights2[j] / (t2 - t2_j); - let w12 = w1 * w2; - - for (k, &t3_k) in t3_coords.iter().enumerate() { - let w3 = weights3[k] / (t3 - t3_k); - let weight = w12 * w3; - - result += weight * f_value(i, j, k); - total_weight += weight; - } - } - } - result / total_weight - } + let term = weights[j] / (t - t_j); + numer += term * f_values[j]; + denom += term; } + + numer / denom } +} - fn init_dimension(&mut self, dim: usize, x_coords: &[f64]) -> Result<(), ValidateError> { +impl Strategy1D for LogChebyshevInterpolation<1> +where + D: Data + RawDataClone + Clone, +{ + fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> { + let x_coords = data.grid[0].as_slice().unwrap(); let n = x_coords.len(); if n < 2 { return Err(ValidateError::Other( - "LogChebyshevInterpolation requires at least 2 grid points per dimension." - .to_string(), + "LogChebyshevInterpolation requires at least 2 grid points.".to_string(), )); } - self.t_coords[dim] = (0..n) + self.t_coords[0] = (0..n) .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) .collect(); - self.weights[dim] = Self::compute_barycentric_weights(n); - - let x_min = *x_coords.first().unwrap(); - let x_max = *x_coords.last().unwrap(); - - if (x_max - x_min).abs() < 1e-15 { - return Err(ValidateError::Other( - "Grid coordinates have zero range".to_string(), - )); - } - self.bounds[dim] = (x_min, x_max); - self.scale_factors[dim] = 1.0 / (x_max - x_min); + self.weights[0] = Self::compute_barycentric_weights(n); Ok(()) } - fn check_bounds_and_transform_generic( - &self, - point: &[f64], - dim_count: usize, - ) -> Result, InterpolateError> { - let mut transformed = vec![0.0; dim_count]; - - for i in 0..dim_count { - let val = point[i]; - let (min, max) = self.bounds[i]; - - if val < min || val > max { - return Err(InterpolateError::Other(format!( - "Input point coordinate {} = {} is outside the grid range [{}, {}]", - i, val, min, max - ))); - } - - transformed[i] = 2.0 * (val - min) * self.scale_factors[i] - 1.0; - } - - Ok(transformed) - } -} - -impl Strategy1D for LogChebyshevInterpolation<1> -where - D: Data + RawDataClone + Clone, -{ - fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> { - let x_coords = data.grid[0].as_slice().unwrap(); - self.init_dimension(0, x_coords) - } - fn interpolate( &self, data: &InterpData1D, point: &[f64; 1], ) -> Result { - let t_coords = self.check_bounds_and_transform_generic(&point[..], 1)?; - let t = t_coords[0]; - + let x = point[0]; + let x_coords = data.grid[0].as_slice().unwrap(); let f_values = data.values.as_slice().unwrap(); - let result = self.interpolate_1d(t, &self.t_coords[0], &self.weights[0], |j| f_values[j]); - Ok(result) + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + + if (x_max - x_min).abs() < 1e-15 { + return Ok(f_values[0]); + } + let t = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; + + Ok(Self::barycentric_interpolate( + t, + &self.t_coords[0], + f_values, + &self.weights[0], + )) } fn allow_extrapolate(&self) -> bool { @@ -1231,7 +1056,17 @@ where fn init(&mut self, data: &InterpData2D) -> Result<(), ValidateError> { for dim in 0..2 { let x_coords = data.grid[dim].as_slice().unwrap(); - self.init_dimension(dim, x_coords)?; + let n = x_coords.len(); + if n < 2 { + return Err(ValidateError::Other( + "LogChebyshevInterpolation requires at least 2 grid points per dimension." + .to_string(), + )); + } + self.t_coords[dim] = (0..n) + .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) + .collect(); + self.weights[dim] = Self::compute_barycentric_weights(n); } Ok(()) } @@ -1241,18 +1076,27 @@ where data: &InterpData2D, point: &[f64; 2], ) -> Result { - let t_coords = self.check_bounds_and_transform_generic(&point[..], 2)?; - let [t_x, t_y] = [t_coords[0], t_coords[1]]; + let [x, y] = *point; + let x_coords = data.grid[0].as_slice().unwrap(); + let y_coords = data.grid[1].as_slice().unwrap(); - let result = self.interpolate_2d( - t_x, - t_y, - &self.t_coords[0], - &self.t_coords[1], - &self.weights[0], - &self.weights[1], - |i, j| data.values[[i, j]], - ); + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + let y_min = *y_coords.first().unwrap(); + let y_max = *y_coords.last().unwrap(); + + let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; + let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; + + let x_coeffs = Self::barycentric_coefficients(t_x, &self.t_coords[0], &self.weights[0]); + let y_coeffs = Self::barycentric_coefficients(t_y, &self.t_coords[1], &self.weights[1]); + + let mut result = 0.0; + for (i, &x_coeff) in x_coeffs.iter().enumerate() { + for (j, &y_coeff) in y_coeffs.iter().enumerate() { + result += x_coeff * y_coeff * data.values[[i, j]]; + } + } Ok(result) } @@ -1269,7 +1113,17 @@ where fn init(&mut self, data: &InterpData3D) -> Result<(), ValidateError> { for dim in 0..3 { let x_coords = data.grid[dim].as_slice().unwrap(); - self.init_dimension(dim, x_coords)?; + let n = x_coords.len(); + if n < 2 { + return Err(ValidateError::Other( + "LogChebyshevInterpolation requires at least 2 grid points per dimension." + .to_string(), + )); + } + self.t_coords[dim] = (0..n) + .map(|j| (PI * (n - 1 - j) as f64 / (n - 1) as f64).cos()) + .collect(); + self.weights[dim] = Self::compute_barycentric_weights(n); } Ok(()) } @@ -1279,21 +1133,34 @@ where data: &InterpData3D, point: &[f64; 3], ) -> Result { - let t_coords = self.check_bounds_and_transform_generic(&point[..], 3)?; - let [t_x, t_y, t_z] = [t_coords[0], t_coords[1], t_coords[2]]; + let [x, y, z] = *point; + let x_coords = data.grid[0].as_slice().unwrap(); + let y_coords = data.grid[1].as_slice().unwrap(); + let z_coords = data.grid[2].as_slice().unwrap(); - let result = self.interpolate_3d( - t_x, - t_y, - t_z, - &self.t_coords[0], - &self.t_coords[1], - &self.t_coords[2], - &self.weights[0], - &self.weights[1], - &self.weights[2], - |i, j, k| data.values[[i, j, k]], - ); + let x_min = *x_coords.first().unwrap(); + let x_max = *x_coords.last().unwrap(); + let y_min = *y_coords.first().unwrap(); + let y_max = *y_coords.last().unwrap(); + let z_min = *z_coords.first().unwrap(); + let z_max = *z_coords.last().unwrap(); + + let t_x = 2.0 * (x - x_min) / (x_max - x_min) - 1.0; + let t_y = 2.0 * (y - y_min) / (y_max - y_min) - 1.0; + let t_z = 2.0 * (z - z_min) / (z_max - z_min) - 1.0; + + let x_coeffs = Self::barycentric_coefficients(t_x, &self.t_coords[0], &self.weights[0]); + let y_coeffs = Self::barycentric_coefficients(t_y, &self.t_coords[1], &self.weights[1]); + let z_coeffs = Self::barycentric_coefficients(t_z, &self.t_coords[2], &self.weights[2]); + + let mut result = 0.0; + for (i, &x_coeff) in x_coeffs.iter().enumerate() { + for (j, &y_coeff) in y_coeffs.iter().enumerate() { + for (k, &z_coeff) in z_coeffs.iter().enumerate() { + result += x_coeff * y_coeff * z_coeff * data.values[[i, j, k]]; + } + } + } Ok(result) } diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index f078e689..ebdde642 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -171,5 +171,5 @@ fn xfxq2_kt_tmdlib() { ]) .assert() .success() - .stdout("0.07892252798564643\n"); + .stdout("0.0655544\n"); } From 7dbe88ff6f2074bdf303201c8c421a3e96a599d5 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 19 Sep 2025 12:46:40 +0200 Subject: [PATCH 11/21] Activate TMDlib workflows in the tests --- neopdf/tests/pdf.rs | 11 ++++++----- neopdf_cli/src/converter.rs | 32 +++++++++++++++++++++++++++++++- neopdf_cli/tests/converter.rs | 1 + neopdf_cli/tests/pdf.rs | 3 ++- 4 files changed, 40 insertions(+), 7 deletions(-) diff --git a/neopdf/tests/pdf.rs b/neopdf/tests/pdf.rs index 12c7051e..03097f2e 100644 --- a/neopdf/tests/pdf.rs +++ b/neopdf/tests/pdf.rs @@ -306,15 +306,16 @@ fn test_pdf_download() { #[ignore] pub fn test_xfxq2_cheby_batch() { // TODO: Use proper grid with Chebyshev interpolations - let pdf = PDF::load("NNPDF40_nnlo_as_01180", 0); + let pdf = PDF::load("MAP22_grids_FF_Km_N3LL.neopdf.lz4", 0); - let ids: Vec = (-4..=4).filter(|&x| x != 0).collect(); - let xs = [1e-5, 1e-3, 1e-3, 1.0]; + let ids: Vec = (-3..=3).filter(|&x| x != 0).collect(); + let kt = 0.1; + let xs = [1e-2, 5e-2, 1e-1, 1.0]; let q2s = [5.0, 10.0, 100.0]; let flatten_points: Vec> = xs .iter() - .flat_map(|&x| q2s.iter().map(move |&q2| vec![x, q2])) + .flat_map(|&x| q2s.iter().map(move |&q2| vec![kt, x, q2])) .collect(); let points_interp: Vec<&[f64]> = flatten_points.iter().map(Vec::as_slice).collect(); let slice_points: &[&[f64]] = &points_interp; @@ -324,7 +325,7 @@ pub fn test_xfxq2_cheby_batch() { let results_single: Vec = slice_points.iter().map(|p| pdf.xfxq2(pid, p)).collect(); for (res_b, res_s) in results_batch.iter().zip(results_single.iter()) { - assert!((res_b - res_s).abs() < LOW_PRECISION); + assert!((res_b - res_s).abs() < LOW_PRECISION, "{res_b} vs. {res_s}"); } } } diff --git a/neopdf_cli/src/converter.rs b/neopdf_cli/src/converter.rs index 760bd3ea..8eb8c5b2 100644 --- a/neopdf_cli/src/converter.rs +++ b/neopdf_cli/src/converter.rs @@ -73,7 +73,37 @@ pub enum Commands { #[arg(short, long)] output: String, }, - /// Update the metadata of the `NeoPDF` grid. + /// Update the metadata of the `NeoPDF` grid + /// + /// Here are the metadata keys that can be updated: + /// + /// - `SetDesc`: description of the set + /// - `SetIndex`: set index + /// - `NumMembers`: number of members + /// - `XMin`: minimum value of x + /// - `XMax`: maximum value of x + /// - `QMin`: minimum value of Q2 + /// - `QMax`: maximum value of Q2 + /// - `Format`: format of the set (can be any string) + /// - `Polarized`: polarization, takes a boolean + /// - `InterpolatorType`: interpolation method (bilinear, logbilinear, logbicubic, logtricubic, chebyshev, interpndlinear) + /// - `ErrorType`: set error (replica, hessian, etc.) + /// - `Particle`: Hadron PID + /// - `FlavorScheme`: flavour scheme (fixed, variable, etc.) + /// - `OrderQCD`: perturbative order + /// - `AlphaS_OrderQCD`: perturbative order of the strong coupling + /// - `MW`: W boson mass + /// - `MZ`: Z boson mass + /// - `MUp`: mass of the up quark + /// - `MDown`: mass of the down quark + /// - `MStrange`: mass of the strange quark + /// - `MCharm`: mass of the charm quark + /// - `MBottom`: mass of the bottom quark + /// - `MTop`: mass of the top quark + /// - `AlphaS_Type`: type of the `AlphaS` interpolation (ipol, etc.) + /// - `NumFlavors`: number of active quark flavours + /// - `SetType`: type of the set (spacelike, timelike) + #[command(verbatim_doc_comment)] Metadata { /// Path to the `NeoPDF` file. #[arg(short, long)] diff --git a/neopdf_cli/tests/converter.rs b/neopdf_cli/tests/converter.rs index 18c605e2..72d94529 100644 --- a/neopdf_cli/tests/converter.rs +++ b/neopdf_cli/tests/converter.rs @@ -12,6 +12,7 @@ Commands: convert Convert a single LHAPDF set to `NeoPDF` format combine-npdfs Combine multiple nuclear PDFs into a single `NeoPDF` with A dependence combine-alphas Combine multiple PDFs with different `alpha_s` values into a single `NeoPDF` + convert-tmd Convert a TMD set to `NeoPDF` format metadata Update the metadata of the `NeoPDF` grid help Print this message or the help of the given subcommand(s) diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index ebdde642..60cf97f1 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -10,6 +10,7 @@ Usage: neopdf compute Commands: xfx_q2 Evaluate xf(x, Q2, pid, ...) for a given set, member, and input values alphas_q2 Evaluate `alphasQ2` for a given set, member, and Q2 value + xfx_q2_kt Evaluate TMD PDF for a given set, member, and input values help Print this message or the help of the given subcommand(s) Options: @@ -171,5 +172,5 @@ fn xfxq2_kt_tmdlib() { ]) .assert() .success() - .stdout("0.0655544\n"); + .stdout("0.07892252798564643\n"); } From 362d92ef0e55c6991a759f54926b3741b71ed045 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 19 Sep 2025 12:51:11 +0200 Subject: [PATCH 12/21] Add more CLI tests for TMDlib --- neopdf_cli/tests/pdf.rs | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index 60cf97f1..240e9a0f 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -173,4 +173,23 @@ fn xfxq2_kt_tmdlib() { .assert() .success() .stdout("0.07892252798564643\n"); + + Command::cargo_bin("neopdf") + .unwrap() + .args([ + "compute", + "xfx_q2", + "--pdf-name", + "MAP22_grids_FF_Km_N3LL.neopdf.lz4", + "--member", + "0", + "--pid", + "2", + "1.0", + "0.1", + "10.0", + ]) + .assert() + .success() + .stdout("0.07899136744063368\n"); } From 1d2ab0a6c18b859a2340add0d5cd3add42067598 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 19 Sep 2025 12:53:35 +0200 Subject: [PATCH 13/21] Bump cache data to activate TMDlib tests --- .github/actions/cache-data/action.yml | 2 +- .github/workflows/run-rust.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/actions/cache-data/action.yml b/.github/actions/cache-data/action.yml index 94fced37..0579579f 100644 --- a/.github/actions/cache-data/action.yml +++ b/.github/actions/cache-data/action.yml @@ -10,7 +10,7 @@ runs: uses: actions/cache@v4 with: path: neopdf-data - key: data-v8 + key: data-v9 - name: Download data if cache miss if: steps.cache-data.outputs.cache-hit != 'true' run: | diff --git a/.github/workflows/run-rust.yml b/.github/workflows/run-rust.yml index 6d2c7374..cdb1396d 100644 --- a/.github/workflows/run-rust.yml +++ b/.github/workflows/run-rust.yml @@ -36,7 +36,7 @@ jobs: RUSTFLAGS: '-Cinstrument-coverage -Clink-dead-code' run: | export NEOPDF_DATA_PATH=${PWD}/neopdf-data - cargo test --workspace --exclude neopdf_pyapi --exclude neopdf_tmdlib --no-fail-fast 2> >(tee stderr 1>&2) + cargo test --workspace --exclude neopdf_pyapi --exclude neopdf_tmdlib --no-fail-fast --features=tmdlib 2> >(tee stderr 1>&2) sed -i 's/\x1B\[[0-9;]\{1,\}[A-Za-z]//g' stderr - name: Generate code coverage From 804f7440816df2534624e1c19ea2a0c9ff1615d4 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 19 Sep 2025 13:29:54 +0200 Subject: [PATCH 14/21] Ignore pure TMDlib tests for the time being --- .github/actions/cache-data/action.yml | 2 +- maintainer/download-data.sh | 9 +++++++++ neopdf_cli/tests/pdf.rs | 6 +++++- 3 files changed, 15 insertions(+), 2 deletions(-) diff --git a/.github/actions/cache-data/action.yml b/.github/actions/cache-data/action.yml index 0579579f..bac12935 100644 --- a/.github/actions/cache-data/action.yml +++ b/.github/actions/cache-data/action.yml @@ -10,7 +10,7 @@ runs: uses: actions/cache@v4 with: path: neopdf-data - key: data-v9 + key: data-v10 - name: Download data if cache miss if: steps.cache-data.outputs.cache-hit != 'true' run: | diff --git a/maintainer/download-data.sh b/maintainer/download-data.sh index 7a4761a8..7d7d99f6 100755 --- a/maintainer/download-data.sh +++ b/maintainer/download-data.sh @@ -22,6 +22,10 @@ NEOPDF_SETS=( nNNPDF30_nlo_as_0118 ) +TMDLIB_SETS=( + MAP22_grids_FF_Km_N3LL +) + # Store the data in the root of the repository cd .. test -d neopdf-data || mkdir neopdf-data @@ -37,3 +41,8 @@ done for neo in "${NEOPDF_SETS[@]}"; do wget --no-verbose --no-clobber -P neopdf-data "https://data.nnpdf.science/neopdf/data/${neo}.neopdf.lz4" done + +# Dowload TMDlib sets +for tmd in "${TMDLIB_SETS[@]}"; do + wget --no-verbose --no-clobber -P neopdf-data "https://data.nnpdf.science/neopdf/data/${tmd}.neopdf.lz4" +done diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index 240e9a0f..9f2dd305 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -153,7 +153,7 @@ fn alphasq2_neopdf() { } #[test] -#[cfg(feature = "tmdlib")] +#[ignore = "Need to find a way to cache TMDlib set."] fn xfxq2_kt_tmdlib() { Command::cargo_bin("neopdf") .unwrap() @@ -173,7 +173,11 @@ fn xfxq2_kt_tmdlib() { .assert() .success() .stdout("0.07892252798564643\n"); +} +#[test] +#[cfg(feature = "tmdlib")] +fn xfxq2_kt_neopdf_tmdlib() { Command::cargo_bin("neopdf") .unwrap() .args([ From 2be64ab542d8aba2a37974a8a256c55a6fdbb2ae Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 19 Sep 2025 13:38:03 +0200 Subject: [PATCH 15/21] Deactivate TMDlib features in Rust tests --- .github/workflows/run-rust.yml | 2 +- neopdf_cli/tests/pdf.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/run-rust.yml b/.github/workflows/run-rust.yml index cdb1396d..6d2c7374 100644 --- a/.github/workflows/run-rust.yml +++ b/.github/workflows/run-rust.yml @@ -36,7 +36,7 @@ jobs: RUSTFLAGS: '-Cinstrument-coverage -Clink-dead-code' run: | export NEOPDF_DATA_PATH=${PWD}/neopdf-data - cargo test --workspace --exclude neopdf_pyapi --exclude neopdf_tmdlib --no-fail-fast --features=tmdlib 2> >(tee stderr 1>&2) + cargo test --workspace --exclude neopdf_pyapi --exclude neopdf_tmdlib --no-fail-fast 2> >(tee stderr 1>&2) sed -i 's/\x1B\[[0-9;]\{1,\}[A-Za-z]//g' stderr - name: Generate code coverage diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index 9f2dd305..d73c8dc1 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -153,7 +153,7 @@ fn alphasq2_neopdf() { } #[test] -#[ignore = "Need to find a way to cache TMDlib set."] +#[cfg(feature = "tmdlib")] fn xfxq2_kt_tmdlib() { Command::cargo_bin("neopdf") .unwrap() From 55bc88ea0ea8f8ec5258fa10e496cba8d0eee17b Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 19 Sep 2025 13:43:52 +0200 Subject: [PATCH 16/21] Fix helpers --- neopdf_cli/tests/converter.rs | 1 - neopdf_cli/tests/pdf.rs | 1 - 2 files changed, 2 deletions(-) diff --git a/neopdf_cli/tests/converter.rs b/neopdf_cli/tests/converter.rs index 72d94529..18c605e2 100644 --- a/neopdf_cli/tests/converter.rs +++ b/neopdf_cli/tests/converter.rs @@ -12,7 +12,6 @@ Commands: convert Convert a single LHAPDF set to `NeoPDF` format combine-npdfs Combine multiple nuclear PDFs into a single `NeoPDF` with A dependence combine-alphas Combine multiple PDFs with different `alpha_s` values into a single `NeoPDF` - convert-tmd Convert a TMD set to `NeoPDF` format metadata Update the metadata of the `NeoPDF` grid help Print this message or the help of the given subcommand(s) diff --git a/neopdf_cli/tests/pdf.rs b/neopdf_cli/tests/pdf.rs index d73c8dc1..0cf9e807 100644 --- a/neopdf_cli/tests/pdf.rs +++ b/neopdf_cli/tests/pdf.rs @@ -10,7 +10,6 @@ Usage: neopdf compute Commands: xfx_q2 Evaluate xf(x, Q2, pid, ...) for a given set, member, and input values alphas_q2 Evaluate `alphasQ2` for a given set, member, and Q2 value - xfx_q2_kt Evaluate TMD PDF for a given set, member, and input values help Print this message or the help of the given subcommand(s) Options: From 893b665f1a76a4631eef6c8b926220e765f8868c Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Sat, 20 Sep 2025 22:31:42 +0200 Subject: [PATCH 17/21] Fix discrepancy between Sequential vs. Batch Chebyshev --- neopdf/tests/pdf.rs | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/neopdf/tests/pdf.rs b/neopdf/tests/pdf.rs index 03097f2e..e5e5c105 100644 --- a/neopdf/tests/pdf.rs +++ b/neopdf/tests/pdf.rs @@ -303,23 +303,31 @@ fn test_pdf_download() { } #[test] -#[ignore] pub fn test_xfxq2_cheby_batch() { - // TODO: Use proper grid with Chebyshev interpolations - let pdf = PDF::load("MAP22_grids_FF_Km_N3LL.neopdf.lz4", 0); - - let ids: Vec = (-3..=3).filter(|&x| x != 0).collect(); - let kt = 0.1; - let xs = [1e-2, 5e-2, 1e-1, 1.0]; - let q2s = [5.0, 10.0, 100.0]; + let pdf = PDF::load("MAP22_nosubgrids.neopdf.lz4", 0); + + let generate_points = |min: f64, max: f64| -> Vec { + (0..20) + .map(|i| min + i as f64 * (max - min) / 19.0) + .collect() + }; + let kts: Vec = generate_points(1e-4, 5.0); + let xs: Vec = generate_points(1e-1, 1.0); + let q2s: Vec = generate_points(1.0, 1e4); + + let mut flatten_points: Vec> = Vec::new(); + for &kt in &kts { + for &x in &xs { + for &q2 in &q2s { + flatten_points.push(vec![kt, x, q2]); + } + } + } - let flatten_points: Vec> = xs - .iter() - .flat_map(|&x| q2s.iter().map(move |&q2| vec![kt, x, q2])) - .collect(); let points_interp: Vec<&[f64]> = flatten_points.iter().map(Vec::as_slice).collect(); let slice_points: &[&[f64]] = &points_interp; + let ids: Vec = (-3..=3).filter(|&x| x != 0).collect(); for &pid in &ids { let results_batch = pdf.xfxq2_cheby_batch(pid, slice_points); let results_single: Vec = slice_points.iter().map(|p| pdf.xfxq2(pid, p)).collect(); From 4c499a32968dc695b5f9b07034db7c7832b5dedc Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Sat, 20 Sep 2025 23:41:07 +0200 Subject: [PATCH 18/21] Fix bug in batch Chebyshev --- neopdf/src/gridpdf.rs | 61 ++++++++++++++++++++++++++++++------------- neopdf/tests/pdf.rs | 20 +++++++++----- 2 files changed, 56 insertions(+), 25 deletions(-) diff --git a/neopdf/src/gridpdf.rs b/neopdf/src/gridpdf.rs index b9f9de95..7729d503 100644 --- a/neopdf/src/gridpdf.rs +++ b/neopdf/src/gridpdf.rs @@ -6,6 +6,7 @@ //! - [`GridArray`]: Stores the full set of subgrids and flavor IDs. use core::panic; +use std::collections::HashMap; use ndarray::{Array1, Array2}; use rayon::prelude::*; @@ -354,11 +355,6 @@ impl GridPDF { return Ok(Vec::new()); } - let subgrid_idx = self.knot_array.find_subgrid(points[0]).ok_or_else(|| { - let (x, q2) = self.get_x_q2(points[0]); - Error::SubgridNotFound { x, q2 } - })?; - let pid_idx = self.knot_array.pid_index(flavor_id).ok_or_else(|| { Error::InterpolationError(format!("Invalid flavor ID: {}", flavor_id)) })?; @@ -369,23 +365,52 @@ impl GridPDF { )); } - let subgrid = &self.knot_array.subgrids[subgrid_idx]; - let log_points: Vec> = points - .iter() - .map(|p| p.iter().map(|&v| v.ln()).collect()) - .collect(); + let mut subgrid_groups: HashMap> = HashMap::new(); + for (i, point) in points.iter().enumerate() { + let subgrid_idx = self.knot_array.find_subgrid(point).ok_or_else(|| { + let (x, q2) = self.get_x_q2(point); + Error::SubgridNotFound { x, q2 } + })?; + + subgrid_groups + .entry(subgrid_idx) + .or_default() + .push((i, *point)); + } - let batch_interpolator = InterpolatorFactory::create_batch_interpolator(subgrid, pid_idx) - .map_err(Error::InterpolationError)?; + let mut all_results: Vec<(usize, f64)> = Vec::new(); - let results = batch_interpolator - .interpolate(log_points) - .map_err(|e| Error::InterpolationError(e.to_string()))?; + for (subgrid_idx, group) in subgrid_groups { + let subgrid = &self.knot_array.subgrids[subgrid_idx]; - Ok(results + let (indices, group_points): (Vec<_>, Vec<_>) = group.into_iter().unzip(); + + let log_points: Vec> = group_points + .iter() + .map(|p| p.iter().map(|&v| v.ln()).collect::>()) + .collect(); + + let batch_interpolator = + InterpolatorFactory::create_batch_interpolator(subgrid, pid_idx) + .map_err(Error::InterpolationError)?; + + let results = batch_interpolator + .interpolate(log_points) + .map_err(|e| Error::InterpolationError(e.to_string()))?; + + for (original_index, result) in indices.into_iter().zip(results) { + all_results.push((original_index, result)); + } + } + + // sort the results according to the original index + all_results.sort_by_key(|&(i, _)| i); + let final_results = all_results .into_iter() - .map(|r| self.apply_force_positive(r)) - .collect()) + .map(|(_, r)| self.apply_force_positive(r)) + .collect(); + + Ok(final_results) } /// Get the values of the momentum fraction `x` and momentum scale `Q2`. diff --git a/neopdf/tests/pdf.rs b/neopdf/tests/pdf.rs index e5e5c105..be4a7a25 100644 --- a/neopdf/tests/pdf.rs +++ b/neopdf/tests/pdf.rs @@ -304,11 +304,12 @@ fn test_pdf_download() { #[test] pub fn test_xfxq2_cheby_batch() { - let pdf = PDF::load("MAP22_nosubgrids.neopdf.lz4", 0); + let pdf = PDF::load("MAP22_grids_FF_Km_N3LL.neopdf.lz4", 0); + let nb_points: usize = 20; let generate_points = |min: f64, max: f64| -> Vec { - (0..20) - .map(|i| min + i as f64 * (max - min) / 19.0) + (0..nb_points) + .map(|i| min + i as f64 * (max - min) / (nb_points as f64 - 1.0)) .collect() }; let kts: Vec = generate_points(1e-4, 5.0); @@ -330,10 +331,15 @@ pub fn test_xfxq2_cheby_batch() { let ids: Vec = (-3..=3).filter(|&x| x != 0).collect(); for &pid in &ids { let results_batch = pdf.xfxq2_cheby_batch(pid, slice_points); - let results_single: Vec = slice_points.iter().map(|p| pdf.xfxq2(pid, p)).collect(); - - for (res_b, res_s) in results_batch.iter().zip(results_single.iter()) { - assert!((res_b - res_s).abs() < LOW_PRECISION, "{res_b} vs. {res_s}"); + let results_seq: Vec = slice_points.iter().map(|p| pdf.xfxq2(pid, p)).collect(); + + for (res_b, res_s) in results_batch.iter().zip(results_seq.iter()) { + assert!( + (res_b - res_s).abs() < LOW_PRECISION, + "PID={pid}: sequential={res_s}, batch={res_b} \n Sequential={:?} \n Batch={:?}", + results_seq, + results_batch + ); } } } From 58c2e8188425d46f9638bc04e57e2cd9ae47c355 Mon Sep 17 00:00:00 2001 From: Tanjona Rabemananjara Date: Sun, 21 Sep 2025 17:04:56 +0200 Subject: [PATCH 19/21] Adress clippy warnings --- neopdf/src/gridpdf.rs | 16 +++++++++------- neopdf/tests/pdf.rs | 7 +------ 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/neopdf/src/gridpdf.rs b/neopdf/src/gridpdf.rs index 7729d503..8e75aa09 100644 --- a/neopdf/src/gridpdf.rs +++ b/neopdf/src/gridpdf.rs @@ -247,7 +247,7 @@ impl GridPDF { Some(ForcePositive::ClipNegative) => value.max(0.0), Some(ForcePositive::ClipSmall) => value.max(1e-10), Some(ForcePositive::NoClipping) => value, - None => value, + _ => value, } } @@ -289,9 +289,10 @@ impl GridPDF { Error::SubgridNotFound { x, q2 } })?; - let pid_idx = self.knot_array.pid_index(flavor_id).ok_or_else(|| { - Error::InterpolationError(format!("Invalid flavor ID: {}", flavor_id)) - })?; + let pid_idx = self + .knot_array + .pid_index(flavor_id) + .ok_or_else(|| Error::InterpolationError(format!("Invalid flavor ID: {flavor_id}")))?; let use_log = matches!( self.info.interpolator_type, @@ -355,9 +356,10 @@ impl GridPDF { return Ok(Vec::new()); } - let pid_idx = self.knot_array.pid_index(flavor_id).ok_or_else(|| { - Error::InterpolationError(format!("Invalid flavor ID: {}", flavor_id)) - })?; + let pid_idx = self + .knot_array + .pid_index(flavor_id) + .ok_or_else(|| Error::InterpolationError(format!("Invalid flavor ID: {flavor_id}")))?; if !matches!(self.info.interpolator_type, InterpolatorType::LogChebyshev) { return Err(Error::InterpolationError( diff --git a/neopdf/tests/pdf.rs b/neopdf/tests/pdf.rs index be4a7a25..6f71ea2b 100644 --- a/neopdf/tests/pdf.rs +++ b/neopdf/tests/pdf.rs @@ -334,12 +334,7 @@ pub fn test_xfxq2_cheby_batch() { let results_seq: Vec = slice_points.iter().map(|p| pdf.xfxq2(pid, p)).collect(); for (res_b, res_s) in results_batch.iter().zip(results_seq.iter()) { - assert!( - (res_b - res_s).abs() < LOW_PRECISION, - "PID={pid}: sequential={res_s}, batch={res_b} \n Sequential={:?} \n Batch={:?}", - results_seq, - results_batch - ); + assert!((res_b - res_s).abs() < LOW_PRECISION); } } } From 6c2b6db3bb213232674a8f7cbdfd2d1a99329d86 Mon Sep 17 00:00:00 2001 From: Tanjona Rabemananjara Date: Sun, 21 Sep 2025 17:17:11 +0200 Subject: [PATCH 20/21] Update change logs --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2afdcdba..0d895eb1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added a logic to compute Chebyshev interpolations in batches (https://github.com/Radonirinaunimi/neopdf/pull/64) - Added the logic to determine the Euclidean distance of a point to the closest subgrid in order to allow extrapolation. - Added version-aware serialization of the `MetaData` struct to ensure backward From 76e28c487e33ab7ba5fbe869ff3fefe93e9e2514 Mon Sep 17 00:00:00 2001 From: Tanjona Rabemananjara Date: Sun, 21 Sep 2025 17:29:38 +0200 Subject: [PATCH 21/21] Update bebchmark --- neopdf/benches/bench_pdf.rs | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/neopdf/benches/bench_pdf.rs b/neopdf/benches/bench_pdf.rs index 9041b803..a73243fc 100644 --- a/neopdf/benches/bench_pdf.rs +++ b/neopdf/benches/bench_pdf.rs @@ -10,6 +10,32 @@ fn xfxq2(c: &mut Criterion) { }); } +fn xfxq2_cheby(c: &mut Criterion) { + let pdf = PDF::load("MAP22_grids_FF_Km_N3LL.neopdf.lz4", 0); + + c.bench_function("xfxq2_cheby", |b| { + b.iter(|| { + pdf.xfxq2( + std::hint::black_box(2), + std::hint::black_box(&[1e-2, 5e-1, 10.0]), + ) + }) + }); +} + +fn xfxq2_cheby_batch(c: &mut Criterion) { + let pdf = PDF::load("MAP22_grids_FF_Km_N3LL.neopdf.lz4", 0); + + c.bench_function("xfxq2_cheby_batch", |b| { + b.iter(|| { + pdf.xfxq2_cheby_batch( + std::hint::black_box(2), + std::hint::black_box(&[&[1e-2, 5e-1, 10.0]]), + ) + }) + }); +} + fn xfxq2s(c: &mut Criterion) { let pdf = PDF::load("NNPDF40_nnlo_as_01180", 0); @@ -45,5 +71,12 @@ fn xfxq2_members(c: &mut Criterion) { }); } -criterion_group!(benches, xfxq2, xfxq2s, xfxq2_members); +criterion_group!( + benches, + xfxq2, + xfxq2s, + xfxq2_members, + xfxq2_cheby, + xfxq2_cheby_batch +); criterion_main!(benches);