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/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<std::size_t> subgrid_shape(subgrid_dim);
+    std::vector<std::size_t> 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<double> subgrid_array(0);
+    pineappl_grid_subgrid_array(grid, 0, 0, 0, subgrid_array.data());
+
     pineappl_grid_write(grid, "advanced-filling.pineappl.lz4");
 
     // release memory

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<double> factors(combinations);
-        std::vector<int> pids(2 * combinations);
+        std::vector<int> 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`

examples/cpp/get-subgrids.cpp

@@ -0,0 +1,95 @@
+#include <pineappl_capi.h>
+
+#include <cassert>
+#include <cstddef>
+#include <iomanip>
+#include <ios>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <numeric>
+
+
+std::vector<std::size_t> unravel_index(std::size_t flat_index, const std::vector<std::size_t>& shape) {
+    std::size_t ndim = shape.size();
+    std::vector<std::size_t> 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<std::size_t>& 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<std::size_t> 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::size_t>());
+            std::vector<double> 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<std::size_t> 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);
+}

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

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/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.