diff --git a/examples/cpp/Makefile b/examples/cpp/Makefile index 7f44b4f6..aef4689b 100644 --- a/examples/cpp/Makefile +++ b/examples/cpp/Makefile @@ -13,6 +13,7 @@ PROGRAMS = \ advanced-filling \ convolve-grid-deprecated \ convolve-grid \ + get-subgrids \ deprecated \ display-channels-deprecated \ display-channels \ @@ -47,6 +48,9 @@ convolve-grid: convolve-grid.cpp deprecated: deprecated.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ +get-subgrids: get-subgrids.cpp + $(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@ + display-channels-deprecated: display-channels-deprecated.cpp $(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@ diff --git a/examples/cpp/advanced-convolution.cpp b/examples/cpp/advanced-convolution.cpp index e16f9e30..16f67eed 100644 --- a/examples/cpp/advanced-convolution.cpp +++ b/examples/cpp/advanced-convolution.cpp @@ -123,5 +123,14 @@ int main(int argc, char* argv[]) { << normalizations.at(bin) << '\n'; } + // modify the particle ID representation of the Grid + pineappl_pid_basis ref_pid_repr_basis = pineappl_grid_pid_basis(grid); + assert(ref_pid_repr_basis == PINEAPPL_PID_BASIS_EVOL); + + pineappl_grid_rotate_pid_basis(grid, PINEAPPL_PID_BASIS_PDG); + + pineappl_pid_basis mod_pid_repr_basis = pineappl_grid_pid_basis(grid); + assert(mod_pid_repr_basis == PINEAPPL_PID_BASIS_PDG); + pineappl_grid_delete(grid); } diff --git a/examples/cpp/advanced-filling.cpp b/examples/cpp/advanced-filling.cpp index 01a1fbfe..8c7a4a72 100644 --- a/examples/cpp/advanced-filling.cpp +++ b/examples/cpp/advanced-filling.cpp @@ -158,6 +158,17 @@ int main() { } //-----------------------------------------------------------------------// + // Check that the Grid contains an empty subgrid at (b, o, c)=(0, 0, 0) + auto subgrid_dim = pineappl_grid_kinematics_len(grid); + std::vector subgrid_shape(subgrid_dim); + std::vector zero_vector(subgrid_dim, 0); + pineappl_grid_subgrid_shape(grid, 0, 0, 0, subgrid_shape.data()); + assert(subgrid_shape == zero_vector); + + // Check that a call to an empty subgrid does not panic + std::vector subgrid_array(0); + pineappl_grid_subgrid_array(grid, 0, 0, 0, subgrid_array.data()); + pineappl_grid_write(grid, "advanced-filling.pineappl.lz4"); // release memory diff --git a/examples/cpp/display-channels.cpp b/examples/cpp/display-channels.cpp index c22da961..636fd328 100644 --- a/examples/cpp/display-channels.cpp +++ b/examples/cpp/display-channels.cpp @@ -22,6 +22,9 @@ int main(int argc, char* argv[]) { // read the grid from a file auto* grid = pineappl_grid_read(filename.c_str()); + // How many convolutions are there? + auto n_conv = pineappl_grid_convolutions_len(grid); + // extract all channels auto* channels = pineappl_grid_channels(grid); @@ -36,7 +39,7 @@ int main(int argc, char* argv[]) { auto combinations = pineappl_channels_combinations(channels, channel); std::vector factors(combinations); - std::vector pids(2 * combinations); + std::vector pids(n_conv * combinations); // read out the channel with index given by `channel`, writing the particle identifiers into // `pids` and the corresponding factors into `factors` diff --git a/examples/cpp/get-subgrids.cpp b/examples/cpp/get-subgrids.cpp new file mode 100644 index 00000000..b37f0db4 --- /dev/null +++ b/examples/cpp/get-subgrids.cpp @@ -0,0 +1,95 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +std::vector unravel_index(std::size_t flat_index, const std::vector& shape) { + std::size_t ndim = shape.size(); + std::vector coords(ndim); + + for (int i = ndim - 1; i >= 0; --i) { + coords[i] = flat_index % shape[i]; + flat_index /= shape[i]; + } + + return coords; +} + +std::string coords_to_string(const std::vector& coords) { + std::ostringstream osstream; + + osstream << "("; + for (std::size_t i = 0; i < coords.size(); ++i) { + osstream << coords[i]; + if (i != coords.size() - 1) { + osstream << ", "; + } + } + osstream << ")"; + + return osstream.str(); +} + + +int main() { + std::string filename = "drell-yan-rap-ll.pineappl.lz4"; + + // read the grid from a file + auto* grid = pineappl_grid_read(filename.c_str()); + + // Determine the number of bins and the index of order and channel + std::size_t n_bins = pineappl_grid_bin_count(grid); + std::size_t order = 0; + std::size_t channel = 0; + + // Get the dimension of the subgrids + std::size_t subgrid_dim = pineappl_grid_kinematics_len(grid); + + std::cout << std::right << std::setw(10) << "bin" << std::setw(10) << "sg idx" + << std::setw(6 * subgrid_dim) << "sg coordinates" << std::setw(16) + << "sg value" << "\n"; + std::cout << std::right << std::setw(10) << "---" << std::setw(10) << "------" + << std::setw(6 * subgrid_dim) << "--------------" << std::setw(16) + << "------------" << "\n"; + + for (std::size_t b = 0; b < n_bins; ++b) { + std::vector subgrid_shape(subgrid_dim); + pineappl_grid_subgrid_shape(grid, b, order, channel, subgrid_shape.data()); + + // Check if the subgrid is not empty + if (subgrid_shape[0] != 0) { + std::size_t flat_shape = std::accumulate(subgrid_shape.begin(), + subgrid_shape.end(), 1, std::multiplies()); + std::vector subgrid_array(flat_shape); + + pineappl_grid_subgrid_array(grid, b, order, channel, subgrid_array.data()); + + for (std::size_t index = 0; index < subgrid_array.size(); ++index) { + if (subgrid_array[index] != 0) { + // Unravel the index to recover the standard coordinates + std::vector coords = unravel_index(index, subgrid_shape); + + std::cout << std::right << std::setw(10) << b << std::setw(10) + << index << std::setw(6 * subgrid_dim) << coords_to_string(coords) + << std::setw(16) << subgrid_array[index] << "\n"; + + // Compare to some reference value + if (b==0 && index==41020) { + assert(subgrid_array[index] == -4.936156925096015e-07); + } + break; + } + } + } + } + + pineappl_grid_delete(grid); +} diff --git a/examples/cpp/output b/examples/cpp/output index e7262e30..d5ad2d4d 100644 --- a/examples/cpp/output +++ b/examples/cpp/output @@ -146,6 +146,32 @@ idx left right dsig/dx dx 21 2.1 2.2 3.507294e-02 0.1 22 2.2 2.3 1.976542e-02 0.1 23 2.3 2.4 5.590552e-03 0.1 + bin sg idx sg coordinates sg value + --- ------ -------------- ------------ + 0 41020 (16, 20, 20) -4.93616e-07 + 1 41020 (16, 20, 20) -4.5499e-08 + 2 40971 (16, 19, 21) -4.10008e-08 + 3 40971 (16, 19, 21) -2.78597e-07 + 4 40971 (16, 19, 21) -2.01924e-07 + 5 40922 (16, 18, 22) -3.39322e-10 + 6 40922 (16, 18, 22) -8.43358e-08 + 7 40922 (16, 18, 22) -2.99247e-07 + 8 40922 (16, 18, 22) -1.5117e-07 + 9 40872 (16, 17, 22) -9.95013e-11 + 10 40873 (16, 17, 23) -1.30578e-07 + 11 40873 (16, 17, 23) -2.55869e-07 + 12 40823 (16, 16, 23) -1.0823e-09 + 13 40823 (16, 16, 23) -4.22686e-09 + 14 40824 (16, 16, 24) -1.68149e-07 + 15 40774 (16, 15, 24) -7.58702e-10 + 16 40774 (16, 15, 24) -2.41887e-08 + 17 40774 (16, 15, 24) -8.83139e-09 + 18 40725 (16, 14, 25) -1.14932e-09 + 19 40725 (16, 14, 25) -3.03857e-08 + 20 40725 (16, 14, 25) -2.66414e-08 + 21 40675 (16, 13, 25) -3.60243e-10 + 22 40676 (16, 13, 26) -1.22785e-08 + 23 40626 (16, 12, 26) -7.04493e-11 idx p-p c#0 l#0 p-p~ c#0 l# p-d c#0 l#0 p-d dx --- ------------ ----------- ------------- ------------ ------ 0 5.263109e-01 5.263109e-01 5.263109e-01 5.263109e-01 0.1 diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index 0c5d9587..f93c544d 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -171,7 +171,12 @@ impl From> for PackedArray { } /// Converts a `multi_index` into a flat index. -fn ravel_multi_index(multi_index: &[usize], shape: &[usize]) -> usize { +/// +/// # Panics +/// +/// TODO +#[must_use] +pub fn ravel_multi_index(multi_index: &[usize], shape: &[usize]) -> usize { assert_eq!(multi_index.len(), shape.len()); multi_index diff --git a/pineappl/src/subgrid.rs b/pineappl/src/subgrid.rs index c115f1f9..799bb51a 100644 --- a/pineappl/src/subgrid.rs +++ b/pineappl/src/subgrid.rs @@ -33,7 +33,7 @@ impl Subgrid for EmptySubgridV1 { Vec::new() } - fn shape(&mut self) -> &[usize] { + fn shape(&self) -> &[usize] { panic!("EmptySubgridV1 doesn't have a shape"); } @@ -166,7 +166,7 @@ impl Subgrid for ImportSubgridV1 { Box::new(self.array.indexed_iter()) } - fn shape(&mut self) -> &[usize] { + fn shape(&self) -> &[usize] { self.array.shape() } @@ -275,7 +275,7 @@ impl Subgrid for InterpSubgridV1 { self.array.is_empty() } - fn shape(&mut self) -> &[usize] { + fn shape(&self) -> &[usize] { self.array.shape() } @@ -471,7 +471,7 @@ pub trait Subgrid { fn scale(&mut self, factor: f64); /// Return the shape of the subgrid - fn shape(&mut self) -> &[usize]; + fn shape(&self) -> &[usize]; /// Assume that the convolution functions for indices `a` and `b` for this grid are the same /// and use this to optimize the size of the grid. diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs index 20b792ef..c28ee961 100644 --- a/pineappl_capi/src/lib.rs +++ b/pineappl_capi/src/lib.rs @@ -60,7 +60,9 @@ use pineappl::boc::{Bin, BinsWithFillLimits, Channel, Kinematics, Order, ScaleFu use pineappl::convolutions::{Conv, ConvType, ConvolutionCache}; use pineappl::grid::{Grid, GridOptFlags}; use pineappl::interpolation::{Interp, InterpMeth, Map, ReweightMeth}; +use pineappl::packed_array::ravel_multi_index; use pineappl::pids::PidBasis; +use pineappl::subgrid::Subgrid; use std::collections::HashMap; use std::ffi::{CStr, CString}; use std::fs::File; @@ -895,6 +897,32 @@ pub unsafe extern "C" fn pineappl_grid_split_lumi(grid: *mut Grid) { grid.split_channels(); } +/// Change the particle ID basis of a given Grid. +/// +/// # Safety +/// +/// If `grid` does not point to a valid `Grid` object (for example when `grid` is the null pointer) +/// or if `pid_basis` does not refer to a correct basis, then this function is not safe to call. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_rotate_pid_basis(grid: *mut Grid, pid_basis: PidBasis) { + let grid = unsafe { &mut *grid }; + + grid.rotate_pid_basis(pid_basis); +} + +/// Get the particle ID basis of a Grid. +/// +/// # 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_pid_basis(grid: *mut Grid) -> PidBasis { + let grid = unsafe { &mut *grid }; + + *grid.pid_basis() +} + /// Optimizes the grid representation for space efficiency. /// /// # Safety @@ -1732,15 +1760,15 @@ pub unsafe extern "C" fn pineappl_grid_split_channels(grid: *mut Grid) { grid.split_channels(); } -/// Similar to `pineappl_lumi_entry` but for luminosity channels that involve 3 or more partons, ie. -/// in the case of multiple convolutions. +/// Read out the channel with index `entry` of the given `channels`. The PDG ids and factors will +/// be copied into `pdg_ids` and `factors`. /// /// # Safety /// -/// The parameter `channels` must point to a valid `Channels` object created by `pineappl_channels_new` or -/// `pineappl_grid_channels`. The parameter `factors` must point to an array as long as the size -/// returned by `pineappl_channels_combinations` and `pdg_ids` must point to an array that is twice as -/// long. +/// The parameter `channels` must point to a valid [`Channels`] object created by +/// [`pineappl_channels_new`] or [`pineappl_grid_channels`]. The parameter `factors` must point to +/// an array as long as the result of [`pineappl_channels_combinations`] and `pdg_ids` must +/// point to an array as long as the result multiplied with the number of convolutions. #[no_mangle] pub unsafe extern "C" fn pineappl_channels_entry( channels: *const Channels, @@ -1750,7 +1778,10 @@ pub unsafe extern "C" fn pineappl_channels_entry( ) { let channels = unsafe { &*channels }; let entry = channels.0[entry].entry(); - let pdg_ids = unsafe { slice::from_raw_parts_mut(pdg_ids, 3 * entry.len()) }; + // if the channel has no entries we assume no convolutions, which is OK we don't copy anything + // in this case + let convolutions = entry.get(0).map_or(0, |x| x.0.len()); + let pdg_ids = unsafe { slice::from_raw_parts_mut(pdg_ids, convolutions * entry.len()) }; let factors = unsafe { slice::from_raw_parts_mut(factors, entry.len()) }; entry @@ -1867,3 +1898,90 @@ pub unsafe extern "C" fn pineappl_grid_convolve( mu_scales, )); } + +/// Get the number of convolutions for a given Grid. +/// +/// # 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_convolutions_len(grid: *mut Grid) -> usize { + let grid = unsafe { &mut *grid }; + + grid.convolutions().len() +} + +/// Get the number of different kinematics for a given Grid. +/// +/// # 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_kinematics_len(grid: *mut Grid) -> usize { + let grid = unsafe { &mut *grid }; + + grid.kinematics().len() +} + +/// Get the shape of a subgrid for a given bin, channel, and order. +/// +/// # 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. Additionally, the pointer that specifies the shape of the +/// subgrid has to be an array whose size must be as given by `pineappl_grid_kinematics_len`. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_subgrid_shape( + grid: *const Grid, + bin: usize, + order: usize, + channel: usize, + shape: *mut usize, +) { + let grid = unsafe { &*grid }; + let subgrid = &grid.subgrids()[[order, bin, channel]]; + let subgrid_shape = if subgrid.is_empty() { + // avoid calling `Subgrid::shape()` for empty grids, which may panic + let subgrid_dim = grid.kinematics().len(); + &vec![0; subgrid_dim] + } else { + subgrid.shape() + }; + let shape = unsafe { slice::from_raw_parts_mut(shape, grid.kinematics().len()) }; + + shape.copy_from_slice(&subgrid_shape); +} + +/// Get the subgrid for a given bin, channel, and order +/// +/// # 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. Additionally, the pointer that specifies the size of the subgrid +/// when flattened must be an array; its size must be computed by multiplying the shape dimension as +/// given by `pineappl_grid_subgrid_shape`. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_subgrid_array( + grid: *const Grid, + bin: usize, + order: usize, + channel: usize, + subgrid_array: *mut f64, +) { + let grid = unsafe { &*grid }; + let subgrid = &grid.subgrids()[[order, bin, channel]]; + + // avoid calling `Subgrid::shape()` for empty grids, which may panic + if !subgrid.is_empty() { + let shape = subgrid.shape(); + let subgrid_array = + unsafe { slice::from_raw_parts_mut(subgrid_array, shape.iter().product()) }; + + for (index, value) in subgrid.indexed_iter() { + let ravel_index = ravel_multi_index(index.as_slice(), &shape); + subgrid_array[ravel_index] = value; + } + } +}