From 73bba9d32034c22e62bc0adde2b71d2380687582 Mon Sep 17 00:00:00 2001 From: Tanjona Rabemananjara Date: Sat, 8 Mar 2025 14:32:35 +0100 Subject: [PATCH] Start polishing codes in the main `pineappl` crate --- pineappl/src/boc.rs | 2 +- pineappl/src/convolutions.rs | 21 +++++----- pineappl/src/evolution.rs | 6 ++- pineappl/src/fk_table.rs | 8 ++-- pineappl/src/grid.rs | 60 ++++++++++++++--------------- pineappl/src/interpolation.rs | 72 ++++++++++++++++++++++++----------- 6 files changed, 98 insertions(+), 71 deletions(-) diff --git a/pineappl/src/boc.rs b/pineappl/src/boc.rs index adf81ba3..dd47dc70 100644 --- a/pineappl/src/boc.rs +++ b/pineappl/src/boc.rs @@ -16,7 +16,7 @@ use std::str::FromStr; /// TODO #[repr(C)] -#[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize)] +#[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize, Hash)] pub enum Kinematics { /// TODO Scale(usize), diff --git a/pineappl/src/convolutions.rs b/pineappl/src/convolutions.rs index 331f7602..63e9c033 100644 --- a/pineappl/src/convolutions.rs +++ b/pineappl/src/convolutions.rs @@ -31,7 +31,7 @@ pub struct ConvolutionCache<'a> { } impl<'a> ConvolutionCache<'a> { - /// TODO + /// Constructor. pub fn new( convolutions: Vec, xfx: Vec<&'a mut dyn FnMut(i32, f64, f64) -> f64>, @@ -160,7 +160,8 @@ impl<'a> ConvolutionCache<'a> { } } -/// TODO +/// A cache for evaluating the full convolutions. Methods like [`Grid::convolve`] and +/// [`FkTable::convolve`] calls this `struct`. pub struct GridConvCache<'a, 'b> { cache: &'b mut ConvolutionCache<'a>, perm: Vec<(usize, bool)>, @@ -171,7 +172,7 @@ pub struct GridConvCache<'a, 'b> { } impl GridConvCache<'_, '_> { - /// TODO + /// Perform the convolution with the PDF(s) for a given channel entry of a subgrid. pub fn as_fx_prod(&mut self, pdg_ids: &[i32], as_order: u8, indices: &[usize]) -> f64 { // TODO: here we assume that // - indices[0] is the (squared) factorization scale, @@ -280,22 +281,22 @@ impl GridConvCache<'_, '_> { } } -/// TODO +/// Defines the various types of convolutions. #[repr(C)] #[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize)] pub enum ConvType { - /// Unpolarized parton distribution function. + /// Unpolarized initial-state hadron. UnpolPDF, - /// Polarized parton distribution function. + /// Polarized initial-state hadron. PolPDF, - /// Unpolarized fragmentation function. + /// Unpolarized final-state hadron. UnpolFF, - /// Polarized fragmentation function. + /// Polarized final-state hadron. PolFF, } impl ConvType { - /// TODO + /// Constructor. #[must_use] pub const fn new(polarized: bool, time_like: bool) -> Self { match (polarized, time_like) { @@ -306,7 +307,7 @@ impl ConvType { } } - /// TODO + /// Check whether the current convolution is of type initial-state hadron. #[must_use] pub const fn is_pdf(&self) -> bool { matches!(self, Self::UnpolPDF | Self::PolPDF) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 4b54d0d5..f62e70d9 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -215,7 +215,11 @@ fn ndarray_from_subgrid_orders_slice( alphas_table: &AlphasTable, ) -> Result { // TODO: remove these assumptions from the following code - assert_eq!(grid.kinematics()[0], Kinematics::Scale(0)); + assert_eq!( + grid.kinematics()[0], + Kinematics::Scale(0), + "For the `Kinematics`, the scale must always go first." + ); let x_start = grid .kinematics() .iter() diff --git a/pineappl/src/fk_table.rs b/pineappl/src/fk_table.rs index 933677c8..cd7fd331 100644 --- a/pineappl/src/fk_table.rs +++ b/pineappl/src/fk_table.rs @@ -107,10 +107,6 @@ impl FkTable { &self.grid } - // TODO: when trying to convert the following function to `const` as per clippy's suggestion, - // the compiler errors out with: 'the destructor for this type cannot be evaluated in constant - // functions' - /// Converts the `FkTable` back to a [`Grid`]. #[must_use] pub fn into_grid(self) -> Grid { @@ -122,7 +118,9 @@ impl FkTable { /// /// # Panics /// - /// TODO + /// This function may panic in two cases: first when there is a mismatch between the `xgrid` + /// and the subgrid nodes, and second when the defined kinematic is not part of the `Kinematic` + /// object. #[must_use] pub fn table(&self) -> ArrayD { let x_grid = self.x_grid(); diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 659e9279..e118aaeb 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -19,16 +19,13 @@ use itertools::Itertools; use lz4_flex::frame::{FrameDecoder, FrameEncoder}; use ndarray::{s, Array2, Array3, ArrayView3, ArrayViewMut3, Axis, CowArray, Dimension, Ix4, Zip}; use serde::{Deserialize, Serialize}; -use std::collections::BTreeMap; +use std::collections::{BTreeMap, HashMap}; use std::io::{BufRead, BufReader, BufWriter, Read, Write}; use std::ops::{Bound, RangeBounds}; use std::{iter, mem}; const BIN_AXIS: Axis = Axis(1); -// const ORDER_AXIS: Axis = Axis(0); -// const CHANNEL_AXIS: Axis = Axis(2); - #[derive(Clone, Deserialize, Serialize)] struct Mmv4; @@ -151,12 +148,13 @@ impl Grid { } } - /// TODO - pub fn reference(&self) -> &Reference { + /// Get the Monte Carlo reference for the given grid. + #[must_use] + pub const fn reference(&self) -> &Reference { &self.reference } - /// TODO + /// Set the Monte Carlo reference for the given grid. pub fn set_reference(&mut self, reference: Reference) { // TODO: check that the number of bins and channels is consistent between the grid and // `reference` @@ -192,10 +190,6 @@ impl Grid { /// and `channel_mask`. A variation of the scales is performed using the factors in `xi`; the /// first factor varies the renormalization scale, the second the factorization scale. Note /// that for the variation to be trusted all non-zero log-grids must be contained. - /// - /// # Panics - /// - /// TODO pub fn convolve( &self, cache: &mut ConvolutionCache, @@ -279,10 +273,6 @@ impl Grid { /// Fills the grid with an ntuple for the given `order`, `observable`, and `channel`. The /// parameter `ntuple` must contain the variables specified by the `kinematics` parameter in /// the constructor [`Grid::new`] in the same order. - /// - /// # Panics - /// - /// TODO pub fn fill( &mut self, order: usize, @@ -438,12 +428,28 @@ impl Grid { Ok(()) } + /// Check if the two `Kinematics` objects are the same even in the case they are ordered + /// differently. + fn are_kinematics_equal(kin_a: &[Kinematics], kin_b: &[Kinematics]) -> bool { + let mut count_a = HashMap::new(); + let mut count_b = HashMap::new(); + + for item in kin_a { + *count_a.entry(item).or_insert(0) += 1; + } + for item in kin_b { + *count_b.entry(item).or_insert(0) += 1; + } + + count_a == count_b + } + /// Merge non-empty `Subgrid`s contained in `other` into `self`. Subgrids with the same bin /// limits are summed and subgrids with non-overlapping bin limits create new bins in `self`. /// /// # Errors /// - /// If `self` and `other` in have different convolutions, PID bases, kinematics, + /// Raises an error if `self` and `other` have different convolutions, PID bases, kinematics, /// interpolations, or scales an error is returned. If the bin limits of `self` and `other` /// are different and if the bin limits of `other` cannot be merged with `self` an error is /// returned. @@ -454,10 +460,6 @@ impl Grid { if self.pid_basis() != other.pid_basis() { return Err(Error::General("PID bases do not match".to_owned())); } - // TODO: relax check if kinematic variables are permutations of each other - if self.kinematics() != other.kinematics() { - return Err(Error::General("kinematics do not match".to_owned())); - } // TODO: relax check if subgrid types don't use interpolation if self.interpolations() != other.interpolations() { return Err(Error::General("interpolations do not match".to_owned())); @@ -465,6 +467,9 @@ impl Grid { if self.scales() != other.scales() { return Err(Error::General("scales do not match".to_owned())); } + if !Self::are_kinematics_equal(self.kinematics(), other.kinematics()) { + return Err(Error::General("kinematics do not match".to_owned())); + } let mut new_orders: Vec = Vec::new(); let mut new_bins = 0; @@ -653,10 +658,6 @@ impl Grid { /// Scales each subgrid by a factor which is the product of the given values `alphas`, `alpha`, /// `logxir`, and `logxif`, each raised to the corresponding powers for each subgrid. In /// addition, every subgrid is scaled by a factor `global` independently of its order. - /// - /// # Panics - /// - /// TODO pub fn scale_by_order( &mut self, alphas: f64, @@ -719,11 +720,12 @@ impl Grid { self.subgrids.view_mut() } - /// TODO + /// Set the bin with filled limits for the grid. This is used to redefine the bin specifications + /// and normalizations. /// /// # Errors /// - /// TODO + /// Raises an error if the length of the bins in the grid and in the redefinition are different. pub fn set_bwfl(&mut self, bwfl: BinsWithFillLimits) -> Result<()> { let bins = bwfl.len(); let grid_bins = self.bwfl().len(); @@ -739,7 +741,7 @@ impl Grid { Ok(()) } - /// TODO + /// Get the bin specifications for this grid. #[must_use] pub const fn bwfl(&self) -> &BinsWithFillLimits { &self.bwfl @@ -992,10 +994,6 @@ impl Grid { } /// Return the metadata of this grid. - /// - /// # Panics - /// - /// TODO #[must_use] pub fn metadata_mut(&mut self) -> &mut BTreeMap { &mut self.metadata diff --git a/pineappl/src/interpolation.rs b/pineappl/src/interpolation.rs index 02211279..067398d6 100644 --- a/pineappl/src/interpolation.rs +++ b/pineappl/src/interpolation.rs @@ -60,35 +60,35 @@ fn lagrange_weights(i: usize, n: usize, u: f64) -> f64 { product / convert::f64_from_usize(factorials) } -/// TODO +/// Represents the method used to reweight the grid. #[repr(C)] #[derive(Clone, Copy, Deserialize, Eq, PartialEq, Serialize)] pub enum ReweightMeth { - /// TODO + /// Reweight the x grid according the `APPLGrid` definition. ApplGridX, - /// TODO + /// No reweighting is applied. NoReweight, } -/// TODO +/// Defines the mapping methods. #[repr(C)] #[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize)] pub enum Map { - /// TODO + /// Uses the `APPLGrid` method [`applgrid::fq20`] for the mapping. ApplGridF2, - /// TODO + /// Uses the `APPLGrid` method [`applgrid::fx2`] for the mapping. ApplGridH0, } -/// TODO +/// Specifies the interpolation method. #[repr(C)] #[derive(Clone, Copy, Deserialize, Eq, PartialEq, Serialize)] pub enum InterpMeth { - /// TODO + /// Uses Lagrange interpolation method. Lagrange, } -/// TODO +/// Defines the interpolation configurations. #[derive(Clone, Deserialize, PartialEq, Serialize)] pub struct Interp { min: f64, @@ -101,7 +101,7 @@ pub struct Interp { } impl Interp { - /// TODO + /// Constructor. /// /// # Panics /// @@ -157,7 +157,7 @@ impl Interp { convert::f64_from_usize(index).mul_add(self.deltay(), self.min) } - /// TODO + /// Reweights a particular x value. #[must_use] pub fn reweight(&self, x: f64) -> f64 { match self.reweight { @@ -166,7 +166,12 @@ impl Interp { } } - /// TODO + /// Interpolates the given x-value and returns the corresponding node index and fractional position. + /// + /// The function maps `x` to `y` using `map_x_to_y` and checks if it falls within the interpolation + /// range. If `y` is outside the range `[min, max]`, `None` is returned. If there is only one node, + /// the function returns `(0, 0.0)`. Otherwise, it calculates the index of the closest node and + /// determines the fractional value `y_fraction` within the interpolation interval. #[must_use] pub fn interpolate(&self, x: f64) -> Option<(usize, f64)> { let y = self.map_x_to_y(x); @@ -191,7 +196,10 @@ impl Interp { } } - /// TODO + /// Reweights a given node given some fractional value. + /// + /// The function applies weights based on the interpolation method set in `interp_meth`. + /// Currently, it supports Lagrange interpolation. #[must_use] pub fn node_weights(&self, fraction: f64) -> ArrayVec { (0..=self.order) @@ -201,13 +209,16 @@ impl Interp { .collect() } - /// TODO + /// Returns the interpolation order. #[must_use] pub const fn order(&self) -> usize { self.order } - /// TODO + /// Returns the values of the interpolation nodes. + /// + /// If there is only one node, the function returns the maps of `min` directly. Otherwise, + /// it returns the maps of the `y`'s values. #[must_use] pub fn node_values(&self) -> Vec { if self.nodes == 1 { @@ -233,43 +244,43 @@ impl Interp { } } - /// TODO + /// Returns the number of nodes. #[must_use] pub const fn nodes(&self) -> usize { self.nodes } - /// TODO + /// Returns the minimum value between the different mappings. #[must_use] pub fn min(&self) -> f64 { self.map_y_to_x(self.min).min(self.map_y_to_x(self.max)) } - /// TODO + /// Returns the maximum value between the different mappings. #[must_use] pub fn max(&self) -> f64 { self.map_y_to_x(self.min).max(self.map_y_to_x(self.max)) } - /// TODO + /// Returns the mapping method used. #[must_use] pub const fn map(&self) -> Map { self.map } - /// TODO + /// Returns the interpolation method used. #[must_use] pub const fn interp_meth(&self) -> InterpMeth { self.interp_meth } - /// TODO + /// Returns the reweighting method used. #[must_use] pub const fn reweight_meth(&self) -> ReweightMeth { self.reweight } - /// TODO + /// Returns the interpolation object within a specified range. #[must_use] pub fn sub_interp(&self, range: Range) -> Self { Self { @@ -284,7 +295,22 @@ impl Interp { } } -/// TODO +/// Performs multi-dimensional interpolation and updates the provided array. +/// +/// # Parameters +/// - `interps`: A slice of `Interp` objects, each representing an interpolation configurations. +/// - `ntuple`: A slice of floating-point values corresponding to interpolation objects. +/// - `weight`: A weight factor applied to the interpolated values. +/// - `array`: A mutable reference to a `PackedArray` where interpolated values will be stored. +/// +/// # Returns +/// - `true` if interpolation is successful and values are updated. +/// - `false` if interpolation fails (e.g., if weight is zero or interpolation points are invalid). +/// +/// # Notes +/// - Ensures that the number of interpolation variables matches the number of provided values. +/// - Uses Cartesian product expansion to compute weighted interpolated values. +/// - Currently allocates memory during iteration; future optimization may be needed. pub fn interpolate( interps: &[Interp], ntuple: &[f64],