diff --git a/CHANGELOG.md b/CHANGELOG.md index 97402030..2198cbf5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- added a `fix_convolution` method to integrate out one of the convolution + dimensions in a grid by convolving it with a non-perturbative function + ## [1.3.0] - 06/12/2025 ### Added diff --git a/examples/cpp/Makefile b/examples/cpp/Makefile index f5e07a99..031cb608 100644 --- a/examples/cpp/Makefile +++ b/examples/cpp/Makefile @@ -22,7 +22,8 @@ PROGRAMS = \ display-orders \ display-orders-deprecated \ merge-grids \ - modify-grid + modify-grid \ + fix-convolution all: $(PROGRAMS) @@ -47,6 +48,9 @@ convolve-grid-deprecated: convolve-grid-deprecated.cpp convolve-grid: convolve-grid.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ +fix-convolution: fix-convolution.cpp + $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ + deprecated: deprecated.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ diff --git a/examples/cpp/fix-convolution.cpp b/examples/cpp/fix-convolution.cpp new file mode 100644 index 00000000..6edeb01d --- /dev/null +++ b/examples/cpp/fix-convolution.cpp @@ -0,0 +1,101 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +void print_results(const std::vector& r1, const std::vector& r2) { + std::cout << std::scientific << std::setprecision(6); + std::cout << "Bin | Original Grid | Fixed Grid | Rel. Diff." << std::endl; + std::cout << "--- | --------------- | --------------- | -----------------" << std::endl; + for (size_t i = 0; i < r1.size(); ++i) { + double rel_diff = 0.0; + if (r1[i] != 0.0) { + rel_diff = std::abs((r1[i] - r2[i]) / r1[i]); + } + std::cout << std::setw(3) << i << " | " + << std::setw(15) << r1[i] << " | " + << std::setw(15) << r2[i] << " | " + << std::setw(17) << rel_diff << std::endl; + assert(rel_diff < 1e-9); + } +} + +int main() { + // Grid with 3 convolutions (2x pol PDF, 1x unpol FF) + const std::string filename = "../../test-data/SIHP-PP-POLARIZED-STAR-NLO.pineappl.lz4"; + const std::string pol_pdf_set = "NNPDFpol11_100"; + const std::string ff_set = "MAPFF10NLOPIsum"; + + LHAPDF::setVerbosity(0); + auto pol_pdf = std::unique_ptr(LHAPDF::mkPDF(pol_pdf_set, 0)); + auto ff = std::unique_ptr(LHAPDF::mkPDF(ff_set, 0)); + + pineappl_grid* grid = pineappl_grid_read(filename.c_str()); + assert(grid != nullptr); + + auto xfx = [](int32_t id, double x, double q2, void* pdf) { + return static_cast (pdf)->xfxQ2(id, x, q2); + }; + auto alphas = [](double q2, void* pdf) { + return static_cast (pdf)->alphasQ2(q2); + }; + + // Convolve original grid + size_t bins = pineappl_grid_bin_count(grid); + std::vector results_orig(bins); + std::vector pdfs_orig_vec = { pol_pdf.get(), pol_pdf.get(), ff.get() }; + void** pdfs_orig = reinterpret_cast(pdfs_orig_vec.data()); + std::vector mu_scales = { 1.0, 1.0, 1.0 }; + pineappl_grid_convolve( + grid, + xfx, + alphas, + pdfs_orig, + pol_pdf.get(), + nullptr, + nullptr, + nullptr, + 1, + mu_scales.data(), + results_orig.data() + ); + + // Fix the third convolution (fragmentation function) + pineappl_grid* grid_fixed = pineappl_grid_fix_convolution(grid, 2, xfx, ff.get(), 1.0); + assert(grid_fixed != nullptr); + + // Convolve the new grid + std::vector results_fixed(bins); + std::vector pdfs_fixed_vec = { pol_pdf.get(), pol_pdf.get() }; + void** pdfs_fixed = reinterpret_cast(pdfs_fixed_vec.data()); + pineappl_grid_convolve( + grid_fixed, + xfx, + alphas, + pdfs_fixed, + pol_pdf.get(), + nullptr, + nullptr, + nullptr, + 1, + mu_scales.data(), + results_fixed.data() + ); + + print_results(results_orig, results_fixed); + + std::cout << "\nSuccess: results from original and fixed grid match." << std::endl; + + pineappl_grid_delete(grid); + pineappl_grid_delete(grid_fixed); + + return 0; +} diff --git a/examples/cpp/fix-convolution.output b/examples/cpp/fix-convolution.output new file mode 100644 index 00000000..c8deb031 --- /dev/null +++ b/examples/cpp/fix-convolution.output @@ -0,0 +1,10 @@ +Bin | Original Grid | Fixed Grid | Rel. Diff. +--- | --------------- | --------------- | ----------------- + 0 | 2.260512e+03 | 2.260512e+03 | 1.207021e-15 + 1 | 1.036130e+03 | 1.036130e+03 | 1.097226e-15 + 2 | 4.894751e+02 | 4.894751e+02 | 1.045182e-15 + 3 | 2.402394e+02 | 2.402394e+02 | 2.011198e-15 + 4 | 1.246446e+02 | 1.246446e+02 | 7.980768e-16 + 5 | 5.268035e+01 | 5.268035e+01 | 8.092688e-16 + +Success: results from original and fixed grid match. diff --git a/examples/object-oriented-cpp/PineAPPL.hpp b/examples/object-oriented-cpp/PineAPPL.hpp index d7fe1de4..b517b28c 100644 --- a/examples/object-oriented-cpp/PineAPPL.hpp +++ b/examples/object-oriented-cpp/PineAPPL.hpp @@ -13,6 +13,7 @@ #include #include #include +#include /** @brief Object oriented interface to PineAPPL.*/ namespace PineAPPL { @@ -123,8 +124,7 @@ struct Grid { /** @brief Underlying raw object. */ pineappl_grid *raw; - /** @brief Constructor (protected to avoid direct instantiation). */ - protected: + /** @brief Constructor. */ Grid(pineappl_grid *grid) : raw(grid) {} /** @brief Deleted copy/move semantics. */ @@ -320,6 +320,31 @@ struct Grid { raw_mu_scales.data(), results.data()); return results; } + + /** + * @brief Fix one of the convolutions in the Grid and return a new Grid with + * lower dimension. + * + * @param conv_idx index of the convolution to fix + * @param pdf the PDF set to use for the convolution + * @param xi the scale factor (xif or xia) + * @return a new grid with one less convolution + */ + std::unique_ptr fix_convolution(std::size_t conv_idx, LHAPDF::PDF *pdf, + double xi) const { + auto xfx = [](std::int32_t id, double x, double q2, void *pdf_ptr) { + return static_cast(pdf_ptr)->xfxQ2(id, x, q2); + }; + + pineappl_grid *new_raw_grid = + pineappl_grid_fix_convolution(this->raw, conv_idx, xfx, pdf, xi); + + if (new_raw_grid == nullptr) { + return nullptr; + } + + return std::unique_ptr(new Grid(new_raw_grid)); + } }; } // namespace PineAPPL diff --git a/pineappl/src/convolutions.rs b/pineappl/src/convolutions.rs index b60f53ea..4031f0ad 100644 --- a/pineappl/src/convolutions.rs +++ b/pineappl/src/convolutions.rs @@ -319,6 +319,12 @@ impl ConvType { pub const fn is_pdf(&self) -> bool { matches!(self, Self::UnpolPDF | Self::PolPDF) } + + /// TODO + #[must_use] + pub const fn is_ff(&self) -> bool { + matches!(self, Self::UnpolFF | Self::PolFF) + } } /// Data type that indentifies different types of convolutions. diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index af45d6f0..efd5bf30 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -6,6 +6,7 @@ use super::error::{Error, Result}; use super::evolution::{self, AlphasTable, EvolveInfo, OperatorSliceInfo}; use super::fk_table::FkTable; use super::interpolation::Interp; +use super::packed_array::PackedArray; use super::pids::PidBasis; use super::reference::Reference; use super::subgrid::{ @@ -1437,6 +1438,229 @@ impl Grid { self.channels = new_channels; } + + /// Fix one of the convolutions in the Grid and return a new Grid with lower convolution dimension. + /// + /// This function integrates out one of the convolution dimensions of the grid by convolving it + /// with the provided function `xfx`. + /// + /// The `conv_idx` parameter specifies which convolution to fix. The `xfx` function provides + /// the values of the parton distribution function or fragmentation function for a given parton + /// ID, `x` value, and scale `mu2`. The `xi` parameter is a scale factor for the factorization + /// or fragmentation scale. + /// + /// # Special handling of fragmentation functions + /// + /// If the convolution being fixed is a fragmentation function, the dependency on the + /// fragmentation scale is removed from the grid. This has a direct impact on the perturbative + /// orders (`Order`). Specifically, the `logxia` of each `Order` is set to zero. + /// + /// As a result, multiple original orders might collapse into a single new order. When this + /// happens, the corresponding subgrids are merged together, ensuring that the total + /// contribution is preserved. The final grid is then optimized to remove any empty or + /// duplicate structures. + /// + /// # Panics + /// + /// This function panics if internal invariants are violated, which typically indicates a bug in + /// the library. + /// + /// # Errors + /// + /// Returns an error if `conv_idx` is out of bounds or if the grid has only one convolution, + /// as the last convolution cannot be fixed. + pub fn fix_convolution( + &self, + conv_idx: usize, + xfx: &mut dyn FnMut(i32, f64, f64) -> f64, + xi: f64, + ) -> Result { + if self.convolutions.len() <= 1 { + return Err(Error::General( + "cannot fix the last convolution".to_string(), + )); + } + + if conv_idx >= self.convolutions.len() { + return Err(Error::General(format!( + "convolution index {} out of bounds (max: {})", + conv_idx, + self.convolutions.len() - 1 + ))); + } + + let mut new_convolutions = self.convolutions.clone(); + new_convolutions.remove(conv_idx); + + let mut new_kinematics = self.kinematics.clone(); + let mut new_interps = self.interps.clone(); + let kin_pos = new_kinematics + .iter() + .position(|k| *k == Kinematics::X(conv_idx)) + .unwrap(); + new_kinematics.remove(kin_pos); + new_interps.remove(kin_pos); + + for kin in &mut new_kinematics { + if let Kinematics::X(i) = kin { + if *i > conv_idx { + *i -= 1; + } + } + } + + let mut new_channel_map: BTreeMap, Vec<(usize, i32, f64)>> = BTreeMap::new(); + for (ichan, chan) in self.channels().iter().enumerate() { + for (pids, factor) in chan.entry() { + let mut new_pids = pids.clone(); + let fixed_pid = new_pids.remove(conv_idx); + + new_channel_map + .entry(new_pids) + .or_default() + .push((ichan, fixed_pid, *factor)); + } + } + + let new_channels: Vec = new_channel_map + .keys() + .map(|pids| Channel::new(vec![(pids.clone(), 1.0)])) + .collect(); + let new_channel_pids: Vec<_> = new_channel_map.keys().cloned().collect(); + + let conv_to_fix = &self.convolutions[conv_idx]; + let (new_orders, order_map) = { + let mut unique_orders = Vec::new(); + + let other_conv_idx_opt = if self.convolutions.len() == 2 { + Some(1 - conv_idx) + } else { + None + }; + + let map: Vec = self + .orders + .iter() + .map(|o| { + let mut new_o = o.clone(); + + if conv_to_fix.conv_type().is_ff() { + new_o.logxia = 0; + } else if let Some(other_conv_idx) = other_conv_idx_opt { + if conv_to_fix.conv_type().is_pdf() + && self.convolutions[other_conv_idx].conv_type().is_ff() + { + new_o.logxif = 0; + } + } + + unique_orders + .iter() + .position(|uo| uo == &new_o) + .unwrap_or_else(|| { + unique_orders.push(new_o); + unique_orders.len() - 1 + }) + }) + .collect(); + (unique_orders, map) + }; + + let mut new_subgrids: Array3 = Array3::from_shape_simple_fn( + (new_orders.len(), self.bwfl.len(), new_channels.len()), + || EmptySubgridV1.into(), + ); + + for (inew_chan, new_pids) in new_channel_pids.iter().enumerate() { + let origins = &new_channel_map[new_pids]; + + (0..self.orders().len()).for_each(|iord| { + for ibin in 0..self.bwfl().bins().len() { + let mut intermediate_sg: Option = None; + + for &(ichan_orig, pid_fixed, factor) in origins { + let sg_orig = &self.subgrids[[iord, ibin, ichan_orig]]; + if sg_orig.is_empty() { + continue; + } + + let mut new_node_values = sg_orig.node_values(); + new_node_values.remove(kin_pos); + + let mut sg_new_array = + PackedArray::new(new_node_values.iter().map(Vec::len).collect()); + + let scale_form = if conv_to_fix.conv_type().is_pdf() { + &self.scales.fac + } else { + &self.scales.frg + }; + + for (mut idxs_orig, val_orig) in sg_orig.indexed_iter() { + let x_val = idxs_orig.remove(kin_pos); + + let sg_orig_node_values = sg_orig.node_values(); + let self_kinematics = self.kinematics(); + + let scale_dims: Vec<_> = sg_orig_node_values + .iter() + .enumerate() + .filter(|(i, _)| { + matches!(self_kinematics.get(*i), Some(Kinematics::Scale(_))) + }) + .map(|(_, v)| v.len()) + .collect(); + let mu2_nodes_calc = + scale_form.calc(&sg_orig_node_values, self_kinematics); + let mu2_idx = scale_form.idx(&idxs_orig, &scale_dims); + let mu2_val = mu2_nodes_calc[mu2_idx] * xi * xi; + + let x = sg_orig_node_values[kin_pos][x_val]; + let pdf_val = xfx(pid_fixed, x, mu2_val) / x; + let final_val = val_orig * factor * pdf_val; + + sg_new_array[idxs_orig.as_slice()] += final_val; + } + + let sg_contrib: SubgridEnum = + ImportSubgridV1::new(sg_new_array, new_node_values).into(); + + if let Some(ref mut isg) = intermediate_sg { + isg.merge(&sg_contrib, None); + } else { + intermediate_sg = Some(sg_contrib); + } + } + + if let Some(sg) = intermediate_sg { + let new_iord = order_map[iord]; + new_subgrids[[new_iord, ibin, inew_chan]].merge(&sg, None); + } + } + }); + } + + let mut new_grid = Self::new( + self.bwfl.clone(), + new_orders, + new_channels, + *self.pid_basis(), + new_convolutions, + new_interps, + new_kinematics, + self.scales.clone(), + ); + new_grid.subgrids = new_subgrids; + new_grid.metadata = self.metadata.clone(); + + new_grid.optimize_using( + GridOptFlags::STRIP_EMPTY_ORDERS + | GridOptFlags::STRIP_EMPTY_CHANNELS + | GridOptFlags::MERGE_SAME_CHANNELS, + ); + + Ok(new_grid) + } } #[cfg(test)] diff --git a/pineappl/tests/drell_yan_lo.rs b/pineappl/tests/drell_yan_lo.rs index 2bcb5f25..760ba0f4 100644 --- a/pineappl/tests/drell_yan_lo.rs +++ b/pineappl/tests/drell_yan_lo.rs @@ -366,26 +366,51 @@ fn perform_grid_tests( assert_approx_eq!(f64, *result, *reference, ulps = 8); } - // TEST 5b: `convolve` with `ConvolutionCache::with_two` - let mut xfx1 = |id, x, q2| pdf.xfx_q2(id, x, q2); - let mut xfx2 = |id, x, q2| pdf.xfx_q2(id, x, q2); - let mut alphas2 = |_| 0.0; - let mut convolution_cache2 = ConvolutionCache::new( - vec![ - Conv::new(ConvType::UnpolPDF, 2212), - Conv::new(ConvType::UnpolPDF, 2212), - ], - vec![&mut xfx1, &mut xfx2], - &mut alphas2, + // TEST 6a: Fix conv_idx = 0 (first hadron) + let mut xfx_fixed_0 = |id, x, q2| pdf.xfx_q2(id, x, q2); + let grid_fixed_0 = grid.fix_convolution(0, &mut xfx_fixed_0, 1.0)?; + + let mut xfx_convolve_0 = |id, x, q2| pdf.xfx_q2(id, x, q2); + let mut alphas_convolve_0 = |_| 0.0; + let mut convolution_cache_one_0 = ConvolutionCache::new( + vec![Conv::new(ConvType::UnpolPDF, 2212)], + vec![&mut xfx_convolve_0], + &mut alphas_convolve_0, + ); + let bins_fixed_0 = grid_fixed_0.convolve( + &mut convolution_cache_one_0, + &[], + &[], + &[], + &[(1.0, 1.0, 1.0)], ); - let bins2 = grid.convolve(&mut convolution_cache2, &[], &[], &[], &[(1.0, 1.0, 1.0)]); - for (result, reference) in bins2.iter().zip(reference.iter()) { - assert_approx_eq!(f64, *result, *reference, ulps = 16); + for (result, original_val) in bins_fixed_0.iter().zip(bins.iter()) { + assert_approx_eq!(f64, *result, *original_val, ulps = 16); } - mem::drop(convolution_cache2); - mem::drop(bins2); + // TEST 6b: Fix conv_idx = 1 (second hadron) + let mut xfx_fixed_1 = |id, x, q2| pdf.xfx_q2(id, x, q2); + let grid_fixed_1 = grid.fix_convolution(1, &mut xfx_fixed_1, 1.0)?; + + let mut xfx_convolve_1 = |id, x, q2| pdf.xfx_q2(id, x, q2); + let mut alphas_convolve_1 = |_| 0.0; + let mut convolution_cache_one_1 = ConvolutionCache::new( + vec![Conv::new(ConvType::UnpolPDF, 2212)], + vec![&mut xfx_convolve_1], + &mut alphas_convolve_1, + ); + let bins_fixed_1 = grid_fixed_1.convolve( + &mut convolution_cache_one_1, + &[], + &[], + &[], + &[(1.0, 1.0, 1.0)], + ); + + for (result, original_val) in bins_fixed_1.iter().zip(bins.iter()) { + assert_approx_eq!(f64, *result, *original_val, ulps = 16); + } // TEST 7a: `optimize_using` - tests `symmetrize` for each subgrid type grid.optimize_using(GridOptFlags::SYMMETRIZE_CHANNELS); diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs index cdab06e6..563efdd0 100644 --- a/pineappl_capi/src/lib.rs +++ b/pineappl_capi/src/lib.rs @@ -1962,6 +1962,32 @@ pub unsafe extern "C" fn pineappl_grid_convolve( )); } +/// Fix one of the convolutions in the Grid and return a new Grid with lower dimension. +/// +/// # Safety +/// +/// TODO +/// +/// # Panics +/// +/// TODO +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_fix_convolution( + grid: *const Grid, + conv_idx: usize, + xfx: extern "C" fn(pdg_id: i32, x: f64, q2: f64, state: *mut c_void) -> f64, + state: *mut c_void, + xi: f64, +) -> Box { + let grid = unsafe { &*grid }; + let mut xfx_closure = |id, x, q2| xfx(id, x, q2, state); + + Box::new( + grid.fix_convolution(conv_idx, &mut xfx_closure, xi) + .unwrap(), + ) +} + /// Get the type of convolutions for this Grid. /// /// # Safety diff --git a/pineappl_cli/src/helpers.rs b/pineappl_cli/src/helpers.rs index 0d85c55c..b3569358 100644 --- a/pineappl_cli/src/helpers.rs +++ b/pineappl_cli/src/helpers.rs @@ -346,7 +346,7 @@ pub fn convolve_scales( .normalizations() .into_iter() .enumerate() - .filter(|(index, _)| (bins.is_empty() || bins.contains(index))) + .filter(|(index, _)| bins.is_empty() || bins.contains(index)) .flat_map(|(_, norm)| iter::repeat(norm).take(scales.len())), ) .for_each(|(value, norm)| *value *= norm); diff --git a/pineappl_cli/src/write.rs b/pineappl_cli/src/write.rs index 57dfdea5..5da5d6ce 100644 --- a/pineappl_cli/src/write.rs +++ b/pineappl_cli/src/write.rs @@ -1,6 +1,6 @@ use super::helpers; use super::{GlobalConfiguration, Subcommand}; -use anyhow::Result; +use anyhow::{Context, Result}; use clap::builder::{PossibleValuesParser, TypedValueParser}; use clap::{ value_parser, Arg, ArgAction, ArgMatches, Args, Command, Error, FromArgMatches, Parser, @@ -36,6 +36,7 @@ enum OpsArg { DeleteKey(String), DeleteOrders(Vec>), DivBinNormDims(Vec), + FixConvolution((usize, String)), MergeBins(Vec>), MergeChannelFactors(bool), MulBinNorm(f64), @@ -62,12 +63,14 @@ struct MoreArgs { impl FromArgMatches for MoreArgs { fn from_arg_matches(matches: &ArgMatches) -> Result { let mut matches = matches.clone(); - let mut args = Vec::new(); + let mut args: Vec> = Vec::new(); let ids: Vec<_> = matches.ids().map(|id| id.as_str().to_owned()).collect(); for id in ids { let indices: Vec<_> = matches.indices_of(&id).unwrap().collect(); - args.resize(indices.iter().max().unwrap() + 1, None); + if args.len() <= *indices.iter().max().unwrap() { + args.resize(indices.iter().max().unwrap() + 1, None); + } match id.as_str() { "cc" => { @@ -215,6 +218,20 @@ impl FromArgMatches for MoreArgs { }); } } + "fix_convolution" => { + for (index, arg) in indices.into_iter().zip( + matches + .remove_occurrences(&id) + .unwrap() + .map(Iterator::collect::>), + ) { + assert_eq!(arg.len(), 1); + let s = &arg[0]; + let (conv_idx_str, pdf) = s.split_once(':').unwrap(); + let conv_idx = conv_idx_str.parse::().unwrap(); + args[index] = Some(OpsArg::FixConvolution((conv_idx, pdf.to_string()))) + } + } "rotate_pid_basis" => { for (index, arg) in indices.into_iter().zip( matches @@ -340,6 +357,14 @@ impl Args for MoreArgs { .value_name("DIM1,...") .value_parser(value_parser!(usize)), ) + .arg( + Arg::new("fix_convolution") + .action(ArgAction::Append) + .help("Fix one of the convolutions with a PDF set") + .long("fix-convolution") + .num_args(1) + .value_name("IDX:CONV_FUN"), + ) .arg( Arg::new("merge_bins") .action(ArgAction::Append) @@ -571,6 +596,13 @@ impl Subcommand for Opts { // UNWRAP: this cannot fail because we only modify the normalizations .unwrap(); } + OpsArg::FixConvolution((conv_idx, pdf_set)) => { + // TODO: account for the variation of scale + let pdf = lhapdf::Pdf::with_setname_and_member(pdf_set, 0) + .with_context(|| format!("Unable to load PDF set '{}'", pdf_set))?; + let mut xfx = |id, x, q2| pdf.xfx_q2(id, x, q2); + grid = grid.fix_convolution(*conv_idx, &mut xfx, 1.0)?; + } OpsArg::MergeBins(ranges) => { // TODO: sort after increasing start indices for range in ranges.iter().rev().cloned() { diff --git a/pineappl_cli/tests/write.rs b/pineappl_cli/tests/write.rs index f15e5689..9fa725bc 100644 --- a/pineappl_cli/tests/write.rs +++ b/pineappl_cli/tests/write.rs @@ -19,6 +19,7 @@ Options: --delete-orders Delete orders with the specified indices --delete-key Delete an internal key-value pair --div-bin-norm-dims Divide each bin normalizations by the bin lengths for the given dimensions + --fix-convolution Fix one of the convolutions with a PDF set --merge-bins Merge specific bins together --merge-channel-factors[=] Merge channel factors into the grid [possible values: true, false] --mul-bin-norm Multiply all bin normalizations with the given factor @@ -444,6 +445,17 @@ const DELETE_ORDERS_STR: &str = "o order 1 O(as^1 a^2 lf^1) "; +const THREE_CONVOLUTIONS_STR: &str = "b pT dsig/dpT (pol) + [GeV] [pb/GeV] +-+-----------------+------------------+-------------- +0 5.108395099639893 6.045444965362549 2.2605116e3 +1 6.045444965362549 6.982494831085205 1.0361301e3 +2 6.982494831085205 7.992245197296143 4.8947508e2 +3 7.992245197296143 8.960753917694092 2.4023939e2 +4 8.960753917694092 9.929026126861572 1.2464463e2 +5 9.929026126861572 11.660773754119873 5.2680349e1 +"; + #[test] fn help() { Command::cargo_bin("pineappl") @@ -1138,3 +1150,31 @@ fn delete_orders() { .success() .stdout(DELETE_ORDERS_STR); } + +#[test] +fn fix_convolution() { + let output = NamedTempFile::new("fix_convolution.pineappl.lz4").unwrap(); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "write", + "--fix-convolution=2:MAPFF10NLOPIsum", + "../test-data/SIHP-PP-POLARIZED-STAR-NLO.pineappl.lz4", + output.path().to_str().unwrap(), + ]) + .assert() + .success() + .stdout(""); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "convolve", + output.path().to_str().unwrap(), + "NNPDFpol11_100+p", + ]) + .assert() + .success() + .stdout(THREE_CONVOLUTIONS_STR); +} diff --git a/pineappl_py/src/grid.rs b/pineappl_py/src/grid.rs index b248afa5..34941203 100644 --- a/pineappl_py/src/grid.rs +++ b/pineappl_py/src/grid.rs @@ -434,6 +434,37 @@ impl PyGrid { .into_pyarray(py) } + /// Fix one of the convolutions in the Grid and return a new Grid with lower dimension. + /// + /// # Panics + /// + /// TODO + /// + /// Parameters + /// ---------- + /// ``conv_idx``: usize + /// index of the convolution (zero-based) + /// ``xfxs`` : callable + /// lhapdf-like callable with arguments `pid, x, Q2` returning x*pdf + /// ``xi``: float + #[must_use] + #[pyo3(signature = (conv_idx, xfx, xi = 1.0))] + pub fn fix_convolution( + &self, + conv_idx: usize, + xfx: Py, + xi: f64, + py: Python<'_>, + ) -> Self { + let mut xfx = move |id: i32, x: f64, q2: f64| { + xfx.call1(py, (id, x, q2)).unwrap().extract(py).unwrap() + }; + + Self { + grid: self.grid.fix_convolution(conv_idx, &mut xfx, xi).unwrap(), + } + } + /// Collect information for convolution with an evolution operator. /// /// # Panics diff --git a/pineappl_py/tests/test_grid.py b/pineappl_py/tests/test_grid.py index a3306aef..ae12ba9a 100644 --- a/pineappl_py/tests/test_grid.py +++ b/pineappl_py/tests/test_grid.py @@ -376,7 +376,10 @@ def test_incosistent_convolutions( xfxs=[pdf.polarized_pdf], # Requires ONE single PDF alphas=pdf.alphasQ, ) - assert "couldn't match Conv { conv_type: UnpolPDF, pid: 2212 } with a convolution function from cache [Conv { conv_type: PolPDF, pid: 2212 }]" == str(err_func.value) + assert ( + "couldn't match Conv { conv_type: UnpolPDF, pid: 2212 } with a convolution function from cache [Conv { conv_type: PolPDF, pid: 2212 }]" + == str(err_func.value) + ) @pytest.mark.parametrize("params,expected", TESTING_SPECS) def test_toy_convolution(self, fake_grids, params, expected): @@ -567,6 +570,110 @@ def test_three_convolutions_with_ff( expected_results, ) + def test_fix_convolution( + self, + pdf, + download_objects, + gridname: str = "SIHP-PP-POLARIZED-STAR-NLO.pineappl.lz4", + ): + expected_results = [ + -3.90292729e09, + +3.43682719e11, + -3.58390524e10, + -4.66855347e10, + -2.15171695e09, + +1.57010877e10, + ] # Numbers computed using `v1.0.0a2` + + grid = download_objects(f"{gridname}") + g = Grid.read(grid) + + # Fix the convolution for Fragmentation Function part + fix_g_conv = g.fix_convolution(conv_idx=2, xfx=pdf.ff_set) + + np.testing.assert_allclose( + fix_g_conv.convolve( + pdg_convs=fix_g_conv.convolutions, + xfxs=[pdf.polarized_pdf, pdf.polarized_pdf, pdf.ff_set], + alphas=pdf.alphasQ, + ), + expected_results, + ) + + def test_fix_convolution_dis(self, fake_grids, pdf): + channels = [Channel([([i], 0.1)]) for i in range(-5, 5)] + + g = fake_grids.grid_with_generic_convolution( + nb_convolutions=1, + channels=channels, + orders=ORDERS, + convolutions=[CONVOBJECT], + ) + + with pytest.raises(BaseException) as err_func: + g.fix_convolution(conv_idx=0, xfx=pdf.unpolarized_pdf) + assert "cannot fix the last convolution" in str(err_func.value) + + def test_fix_convolution_wrong_index(self, fake_grids, pdf): + channels = [Channel([([i, -i], 0.1)]) for i in range(-5, 5)] + + g = fake_grids.grid_with_generic_convolution( + nb_convolutions=2, + channels=channels, + orders=ORDERS, + convolutions=[CONVOBJECT, CONVOBJECT], + ) + + with pytest.raises(BaseException) as err_func: + g.fix_convolution(conv_idx=4, xfx=pdf.unpolarized_pdf) + assert "convolution index 4 out of bounds (max: 1)" in str(err_func.value) + + def test_fix_convolution_logxia(self, fake_grids, pdf): + rndgen = Generator(PCG64(seed=1234)) + binning = [1e-2, 1e-1, 0.5, 1] + + channels = [Channel([([i, -i, i], 0.1)]) for i in range(-5, 5)] + orders = [Order(2, 0, 0, 0, 0), Order(2, 0, 0, 0, 1)] + + convbools = [(False, False), (False, False), (False, True)] + hpids = [2212, 2212, 211] + convtypes = [ConvType(polarized=p, time_like=t) for p, t in convbools] + convolutions = [ + Conv(convolution_types=c, pid=p) for c, p in zip(convtypes, hpids) + ] + + g = fake_grids.grid_with_generic_convolution( + nb_convolutions=len(convolutions), + channels=channels, + orders=orders, + convolutions=convolutions, + bins=binning, + ) + + # Fill with non-empty subgrids + _q2grid = np.geomspace(1e3, 1e5, 5) + _xgrid = np.geomspace(1e-5, 1, 4) + comb_nodes = [_q2grid] + [_xgrid for _ in convolutions] + ntuples = [np.array(list(kins)) for kins in itertools.product(*comb_nodes)] + obs = [rndgen.uniform(binning[0], binning[-1]) for _ in ntuples] + for pto in range(len(ORDERS)): + for channel_id in range(len(channels)): + g.fill_array( + order=pto, + observables=obs, + channel=channel_id, + ntuples=ntuples, + weights=np.repeat(1, len(obs)), + ) + + # Fix the Fragmentation Function + g_fix = g.fix_convolution(conv_idx=2, xfx=pdf.ff_set) + orders_fix = g_fix.orders() + + # Check that the orders have been merged + assert len(orders_fix) == 1 + assert orders_fix[0].as_tuple() == (2, 0, 0, 0, 0) + def test_many_convolutions(self, fake_grids, pdf, nb_convolutions: int = 3): """Test for fun many convolutions.""" expected_results = [