From 9d5cbb5ccc159ed7b42e190e945bc4fc942b2e63 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Thu, 24 Jul 2025 14:16:29 +0200 Subject: [PATCH 1/2] Add an interface to set subgrid int the C-API --- examples/cpp/Makefile | 4 + examples/cpp/set-subgrids.cpp | 218 +++++++++++++++++++++++++++++++ examples/cpp/set-subgrids.output | 4 + pineappl_capi/src/lib.rs | 103 ++++++++++++++- 4 files changed, 327 insertions(+), 2 deletions(-) create mode 100644 examples/cpp/set-subgrids.cpp create mode 100644 examples/cpp/set-subgrids.output diff --git a/examples/cpp/Makefile b/examples/cpp/Makefile index f5e07a99..3c02602d 100644 --- a/examples/cpp/Makefile +++ b/examples/cpp/Makefile @@ -14,6 +14,7 @@ PROGRAMS = \ convolve-grid-deprecated \ convolve-grid \ get-subgrids \ + set-subgrids \ evolve-grid \ evolve-grid-identity \ deprecated \ @@ -53,6 +54,9 @@ deprecated: deprecated.cpp get-subgrids: get-subgrids.cpp $(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@ +set-subgrids: set-subgrids.cpp + $(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@ + evolve-grid: evolve-grid.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ diff --git a/examples/cpp/set-subgrids.cpp b/examples/cpp/set-subgrids.cpp new file mode 100644 index 00000000..a3e1583b --- /dev/null +++ b/examples/cpp/set-subgrids.cpp @@ -0,0 +1,218 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using KinematicsTuple = std::tuple; + +struct KinInterpolation { + std::vector x_interp; + std::vector z_interp; +}; + +template +std::vector geomspace(T start, T stop, int num, bool endpoint = false) { + std::vector result(num); + + if (num == 1) { + result[0] = start; + return result; + } + + T log_start = std::log(start); + T log_stop = std::log(stop); + T step = (log_stop - log_start) / (endpoint ? (num - 1) : num); + + for (int i = 0; i < num; ++i) { + result[i] = std::exp(log_start + i * step); + } + + return result; +} + +KinInterpolation kinematics_interpolation_points() { + std::vector x = geomspace(1e-5, 1.0, 50); + std::vector z = geomspace(1e-5, 1.0, 50); + + return {x, x}; +} + +std::vector generate_subgrid_arrays( + std::mt19937& rng, + std::vector x, + std::vector z +) { + std::size_t kin_length = x.size() * z.size(); + std::vector subgrid(kin_length); + + // NOTE: `subgrid` is a flatten matrix whose layout was `[q2=1][x][z]`. + // The order of the kinematics shoud match the kinematics declaration. + for (std::size_t i = 0; i != kin_length; i++) { + subgrid[i] = std::generate_canonical(rng); + } + + return subgrid; +} + +void fill_grid(pineappl_grid* grid, std::vector kin_tuples) { + auto rng = std::mt19937(); + + auto* channels = pineappl_grid_channels(grid); + std::size_t n_bins = pineappl_grid_bin_count(grid); + std::size_t n_orders = pineappl_grid_order_count(grid); + std::size_t n_channels = pineappl_channels_count(channels); + + // Get the kinematics + KinInterpolation kins = kinematics_interpolation_points(); + std::vector x_interp = kins.x_interp; + std::vector z_interp = kins.z_interp; + + // Extract the shape of the subgrid - Q2 always passed as an array of ONE element + std::vector subgrid_shape = {1, x_interp.size(), z_interp.size()}; + + for (std::size_t b = 0; b != n_bins; b++) { + for (std::size_t o = 0; o != n_orders; o++) { + for (std::size_t c = 0; c != n_channels; c++) { + // Construct the node values of {Q2, x_inter, z_interp} + // NOTE: Pay attention to the order, it should match the kinematics declaration + // and how the subgrid was constructed (see `generate_subgrid_array`). + std::vector node_values = { std::get<0>(kin_tuples[b]) }; + node_values.insert(node_values.end(), x_interp.begin(), x_interp.end()); + node_values.insert(node_values.end(), z_interp.begin(), z_interp.end()); + + // Mock the subgrids for a given bin, order, and channel + std::vector subgrid_arrays = generate_subgrid_arrays(rng, x_interp, z_interp); + + // set the subgrids + pineappl_grid_set_subgrid( + grid, + b, o, c, // kinematics index + node_values.data(), + subgrid_arrays.data(), + subgrid_shape.data() + ); + } + } + } +} + +int main() { + // --- + // Create all channels + + std::size_t nb_convolutions = 2; + auto* channels = pineappl_channels_new(nb_convolutions); + + int32_t pids1[] = { 21, 21 }; + double factors1[] = { 1.0 }; + // define the channel #0 + pineappl_channels_add(channels, 1, pids1, factors1); + + // create another channel this channel is the down-type-antidown-type quark channel; here we + int32_t pids2[] = { 1, -1, 3, -3, 5, -5 }; + // define the channel #1 + pineappl_channels_add(channels, 3, pids2, nullptr); + + // --- + // Specify the perturbative orders that will be filled into the grid + std::vector orders = { + 1, 0, 0, 0, 0, // order #0: LO QCD + 2, 0, 0, 0, 0, // order #1: NLO QCD + }; + + // --- + // Specify the bin limits + + // In SIDIS, a bin is defined as a tuple (Q2, x, z) values (3D). + std::vector kin_obs = {std::make_tuple(1e3, 1e-5, 1e-2), std::make_tuple(1e4, 1e-2, 1e-3)}; + // We are going to define some placeholder 1D bins that we'll overwrite later. + std::vector bins; + for (std::size_t i = 0; i < kin_obs.size() + 1; ++i) { + bins.push_back(static_cast(i)); + } + + // --- + // Construct the objects that are needed to fill the Grid + + pineappl_pid_basis pid_basis = PINEAPPL_PID_BASIS_EVOL; + pineappl_conv convs[] = { + { PINEAPPL_CONV_TYPE_UNPOL_PDF, 2212 }, + { PINEAPPL_CONV_TYPE_UNPOL_FF, 211 }, // Assumes Pion + }; + + // Define the kinematics required for this process. In the following example we have ONE single + // scale and two momentum fractions (corresponding to the two initial- and final-state hadrons). + // The format of the kinematics is: { type, value }. + pineappl_kinematics scales = { PINEAPPL_KINEMATICS_SCALE, 0 }; + pineappl_kinematics x1 = { PINEAPPL_KINEMATICS_X, 0 }; + pineappl_kinematics x2 = { PINEAPPL_KINEMATICS_X, 1 }; + pineappl_kinematics kinematics[3] = { scales, x1, x2 }; + + // Define the specificities of the interpolations for each of the kinematic variables. + pineappl_reweight_meth scales_reweight = PINEAPPL_REWEIGHT_METH_NO_REWEIGHT; // Reweighting method + pineappl_reweight_meth moment_reweight = PINEAPPL_REWEIGHT_METH_APPL_GRID_X; + pineappl_map scales_mapping = PINEAPPL_MAP_APPL_GRID_H0; // Mapping method + pineappl_map moment_mapping = PINEAPPL_MAP_APPL_GRID_F2; + pineappl_interp_meth interpolation_meth = PINEAPPL_INTERP_METH_LAGRANGE; + pineappl_interp interpolations[3] = { + { 1e2, 1e8, 40, 3, scales_reweight, scales_mapping, interpolation_meth }, // Interpolation fo `scales` + { 2e-7, 1.0, 50, 3, moment_reweight, moment_mapping, interpolation_meth }, // Interpolation fo `x1` + { 2e-7, 1.0, 50, 3, moment_reweight, moment_mapping, interpolation_meth }, // Interpolation fo `x2` + }; + + // Define the unphysical scale objects + pineappl_scale_func_form scale_mu = { PINEAPPL_SCALE_FUNC_FORM_SCALE, 0 }; + pineappl_scale_func_form mu_scales[3] = { scale_mu, scale_mu, scale_mu }; + + // --- + // Create the grid using the previously set information about orders, bins and channels + + auto* grid = pineappl_grid_new2(bins.size() - 1, bins.data(), orders.size() / 5, orders.data(), + channels, pid_basis, convs, 3, interpolations, kinematics, mu_scales); + + pineappl_channels_delete(channels); + + // --- + // Fill the grid with phase-space points + fill_grid(grid, kin_obs); + + // --- + // NOTE: We now have to remap the bins to 3D + + // We need to flatten the array of `KinematicsTuple` and define the normalizations + std::vector flat_kin_obs; + for (const auto& kin : kin_obs) { + flat_kin_obs.push_back(std::get<0>(kin)); + flat_kin_obs.push_back(std::get<1>(kin)); + flat_kin_obs.push_back(std::get<2>(kin)); + } + std::vector normalizations(bins.size() - 1, 1.0); + pineappl_grid_set_bwfl( + grid, + flat_kin_obs.data(), // lower bin limits + flat_kin_obs.data(), // upper bin limits + bins.size() - 1, // number of bins + 3, // Dimension of the bins (Q2, x, z) + normalizations.data() + ); + pineappl_grid_optimize(grid); + + // --- + // Write the grid to disk - the filename can be anything ... + std::string filename = "sidis-toygrid.pineappl.lz4"; + pineappl_grid_write(grid, filename.c_str()); + + // destroy the object + pineappl_grid_delete(grid); + + std::cout << "Generated " << filename << " containing a toy SIDIS.\n\n" + "Try running the following command to check the bins:\n" + " - pineappl read --bins " << filename << "\n"; +} diff --git a/examples/cpp/set-subgrids.output b/examples/cpp/set-subgrids.output new file mode 100644 index 00000000..5bc57c99 --- /dev/null +++ b/examples/cpp/set-subgrids.output @@ -0,0 +1,4 @@ +Generated sidis-toygrid.pineappl.lz4 containing a toy SIDIS. + +Try running the following command to check the bins: + - pineappl read --bins sidis-toygrid.pineappl.lz4 diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs index d28d213e..96a69c92 100644 --- a/pineappl_capi/src/lib.rs +++ b/pineappl_capi/src/lib.rs @@ -63,9 +63,10 @@ use pineappl::evolution::{AlphasTable, OperatorSliceInfo}; use pineappl::fk_table::{FkAssumptions, FkTable}; use pineappl::grid::{Grid, GridOptFlags}; use pineappl::interpolation::{Interp as InterpMain, InterpMeth, Map, ReweightMeth}; -use pineappl::packed_array::ravel_multi_index; +use pineappl::packed_array::PackedArray; +use pineappl::packed_array::{ravel_multi_index, unravel_index}; use pineappl::pids::PidBasis; -use pineappl::subgrid::Subgrid; +use pineappl::subgrid::{ImportSubgridV1, Subgrid}; use std::collections::HashMap; use std::ffi::{CStr, CString}; use std::fs::File; @@ -2065,6 +2066,104 @@ pub unsafe extern "C" fn pineappl_grid_subgrid_array( } } +/// Set the subgrid of a Grid for a given bin, order, and channel. +/// +/// # Safety +/// +/// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer, +/// this function is not safe to call. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_set_subgrid( + grid: *mut Grid, + bin: usize, + order: usize, + channel: usize, + node_values: *mut f64, + subgrid_array: *mut f64, + subgrid_shape: *mut usize, +) { + let grid = unsafe { &mut *grid }; + let num_kins = grid.kinematics().len(); + + let subgrid_shape = unsafe { slice::from_raw_parts(subgrid_shape, num_kins) }; + let subgrid_array = + unsafe { slice::from_raw_parts(subgrid_array, subgrid_shape.iter().product()) }; + + let node_values: Vec> = { + let mut offset = 0; + subgrid_shape + .iter() + .map(|&dim_size| { + let dim_nodes = + unsafe { slice::from_raw_parts(node_values.add(offset), dim_size) }.to_vec(); + offset += dim_size; + dim_nodes + }) + .collect() + }; + + let mut sparse_array: PackedArray = + PackedArray::new(node_values.iter().map(Vec::len).collect()); + + for (index, value) in subgrid_array + .iter() + .enumerate() + .filter(|(_, value)| **value != 0.0) + { + let index_unravel = unravel_index(index, subgrid_shape); + sparse_array[index_unravel.as_slice()] = *value; + } + + let subgrid = ImportSubgridV1::new(sparse_array, node_values); + grid.subgrids_mut()[[order, bin, channel]] = subgrid.into(); +} + +/// Redefine the bin representation of the Grid; generalization of `pineappl_grid_set_remapper`. +/// +/// # Panics +/// +/// TODO +/// +/// # Safety +/// +/// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer, +/// this function is not safe to call. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_set_bwfl( + grid: *mut Grid, + bins_lower: *mut f64, + bins_upper: *mut f64, + n_bins: usize, + bin_dim: usize, + normalizations: *mut f64, +) { + let grid = unsafe { &mut *grid }; + + let bins_lower = unsafe { slice::from_raw_parts(bins_lower, bin_dim * n_bins) }; + let bins_upper = unsafe { slice::from_raw_parts(bins_upper, bin_dim * n_bins) }; + let normalizations = unsafe { slice::from_raw_parts(normalizations, n_bins) }; + + let limits: Vec> = bins_lower + .chunks(bin_dim) + .zip(bins_upper.chunks(bin_dim)) + .map(|(bl_chunk, bu_chunk)| { + bl_chunk + .iter() + .zip(bu_chunk.iter()) + .map(|(&a, &b)| (a, b)) + .collect() + }) + .collect(); + + grid.set_bwfl( + BinsWithFillLimits::from_limits_and_normalizations(limits, normalizations.to_vec()) + // UNWRAP: error handling in the CAPI is to abort + .unwrap(), + ) + // UNWRAP: error handling in the CAPI is to abort + .unwrap(); +} + /// Get the shape of the objects represented in the evolve info. /// /// # Safety From 1468365b83d913d06559358ea8430a5c85b5549d Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Thu, 24 Jul 2025 14:27:46 +0200 Subject: [PATCH 2/2] Fix memory leak --- examples/cpp/set-subgrids.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/cpp/set-subgrids.cpp b/examples/cpp/set-subgrids.cpp index a3e1583b..e629a9b5 100644 --- a/examples/cpp/set-subgrids.cpp +++ b/examples/cpp/set-subgrids.cpp @@ -101,6 +101,9 @@ void fill_grid(pineappl_grid* grid, std::vector kin_tuples) { } } } + + // Remove channels object from memory + pineappl_channels_delete(channels); } int main() {