Exposes more functions to C-API for APFEL++

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 $@
 

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);
 }

examples/cpp/get-subgrids.cpp

@@ -0,0 +1,54 @@
+#include <pineappl_capi.h>
+
+#include <cstddef>
+#include <iomanip>
+#include <ios>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <numeric>
+
+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_convolutions_len(grid) + 1;
+
+    std::cout << std::right << std::setw(10) << "bin" << std::setw(10) << "sg idx"
+        << std::setw(16) << "sg value" << "\n";
+    std::cout << std::right << std::setw(10) << "---" << std::setw(10) << "------"
+        << std::setw(16) << "------------" << "\n";
+
+    for (std::size_t b = 0; b < n_bins; ++b) {
+        std::vector<std::size_t> subgrid_shape(subgrid_dim);
+        pineappl_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::size_t>());
+            std::vector<double> subgrid_array(flat_shape);
+
+            pineappl_subgrid_array(grid, b, order, channel, subgrid_array.data());
+
+            // TODO: Unravel of the index & check against some reference values
+            for (std::size_t value = 0; value < subgrid_array.size(); ++value) {
+                if (subgrid_array[value] != 0) {
+                    std::cout << std::right << std::setw(10) << b << std::setw(10)
+                        << value << std::setw(16) << subgrid_array[value] << "\n";
+                    break;
+                }
+            }
+        }
+    }
+
+    pineappl_grid_delete(grid);
+}

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 value
+       ---    ------    ------------
+         0     41020    -4.93616e-07
+         1     41020     -4.5499e-08
+         2     40971    -4.10008e-08
+         3     40971    -2.78597e-07
+         4     40971    -2.01924e-07
+         5     40922    -3.39322e-10
+         6     40922    -8.43358e-08
+         7     40922    -2.99247e-07
+         8     40922     -1.5117e-07
+         9     40872    -9.95013e-11
+        10     40873    -1.30578e-07
+        11     40873    -2.55869e-07
+        12     40823     -1.0823e-09
+        13     40823    -4.22686e-09
+        14     40824    -1.68149e-07
+        15     40774    -7.58702e-10
+        16     40774    -2.41887e-08
+        17     40774    -8.83139e-09
+        18     40725    -1.14932e-09
+        19     40725    -3.03857e-08
+        20     40725    -2.66414e-08
+        21     40675    -3.60243e-10
+        22     40676    -1.22785e-08
+        23     40626    -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

pineappl/src/packed_array.rs

@@ -171,7 +171,12 @@ impl<T: Copy + Default + PartialEq> From<ArrayViewD<'_, T>> for PackedArray<T> {
 }
 
 /// 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

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
@@ -1867,3 +1895,77 @@ 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_convolutions_len(grid: *mut Grid) -> usize {
+    let grid = unsafe { &mut *grid };
+    grid.convolutions().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 must be an array whose size must be `(n+1)` where `n` represents the number of
+/// convolutions as given by `pineappl_convolutions_len`.
+#[no_mangle]
+pub unsafe extern "C" fn pineappl_subgrid_shape(
+    grid: *mut Grid,
+    bin: usize,
+    order: usize,
+    channel: usize,
+    shape: *mut usize,
+) {
+    let grid = unsafe { &mut *grid };
+    let mut subgrid = grid.subgrids()[[order, bin, channel]].clone();
+
+    let subgrid_shape = if subgrid.is_empty() {
+        let subgrid_dim = grid.convolutions().len() + 1;
+        &vec![0; subgrid_dim]
+    } else {
+        subgrid.shape()
+    };
+
+    let shape = unsafe { slice::from_raw_parts_mut(shape, grid.convolutions().len() + 1) };
+    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_subgrid_shape`.
+#[no_mangle]
+pub unsafe extern "C" fn pineappl_subgrid_array(
+    grid: *mut Grid,
+    bin: usize,
+    order: usize,
+    channel: usize,
+    subgrid_array: *mut f64,
+) {
+    let grid = unsafe { &mut *grid };
+    let mut subgrid = grid.subgrids()[[order, bin, channel]].clone();
+
+    let flattened_shape = subgrid.shape().iter().copied().product();
+    let mut flattened_subgrid_array = vec![0.0; flattened_shape];
+
+    for (index, value) in subgrid.indexed_iter() {
+        let ravel_index = ravel_multi_index(index.as_slice(), subgrid.clone().shape());
+        flattened_subgrid_array[ravel_index] = value;
+    }
+
+    let subgrid_array = unsafe { slice::from_raw_parts_mut(subgrid_array, flattened_shape) };
+
+    subgrid_array.copy_from_slice(&flattened_subgrid_array);
+}