Pass the computation of a Operator slice as a callback

Author: Radonirinaunimi

examples/cpp/evolve-grid-identity.cpp

@@ -24,27 +24,36 @@ std::vector<std::size_t> unravel_index(std::size_t flat_index, const std::vector
     return coords;
 }
 
-std::vector<double> generate_fake_ekos(
-    std::vector<int> pids_in,
-    std::vector<double> x_in,
-    std::vector<int> pids_out,
-    std::vector<double> x_out
+extern "C" void generate_fake_ekos(
+    const int* pids_in,
+    const double* x_in,
+    const int* pids_out,
+    const double* x_out,
+    double* eko_buffer,
+    pineappl_conv_type conv_type,
+    std::size_t pids_in_len,
+    std::size_t x_in_len,
+    std::size_t pids_out_len,
+    std::size_t x_out_len
 ) {
-    std::size_t flat_len = x_out.size() * x_in.size() * pids_out.size() * pids_in.size();
-    std::vector<double> ops(flat_len);
-
+    // Ignore unused variables
+    (void)pids_in;
+    (void)x_in;
+    (void)pids_out;
+    (void)x_out;
+    (void) conv_type;
+
+    std::size_t flat_len = pids_in_len * x_in_len * pids_out_len * x_out_len;
     // NOTE: The EKO has to have as shape: (pids_in, x_in, pids_out, x_out)
-    std::vector<std::size_t> shape = {pids_in.size(), x_in.size(), pids_out.size(), x_out.size()};
+    std::vector<std::size_t> shape = {pids_in_len, x_in_len, pids_out_len, x_out_len};
     for (std::size_t i = 0; i != flat_len; i++) {
         std::vector<std::size_t> coords = unravel_index(i, shape);
 
         double delta_ik = (coords[0] == coords[2]) ? 1.0 : 0.0;
         double delta_jl = (coords[1] == coords[3]) ? 1.0 : 0.0;
 
-        ops[i] = delta_ik * delta_jl;
+        eko_buffer[i] = delta_ik * delta_jl;
     }
-
-    return ops;
 }
 
 void print_results(std::vector<double> dxsec_grid, std::vector<double> dxsec_fktable) {
@@ -152,18 +161,6 @@ int main() {
     // std::vector<int> pids_out = {-6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 21, 22};
     std::vector<int> pids_out = pids_in;
 
-    // The Evolution Operator is a vector with length `N_conv * N_Q2_slices * Σ product(OP shape)`
-    std::vector<double> op_slices;
-    std::size_t flat_len = x_in.size() * x_in.size() * pids_in.size() * pids_out.size();
-    for (std::size_t _i = 0; _i != conv_types.size(); _i++) {
-        for (std::size_t j = 0; j != fac1.size(); j++) {
-            std::vector<double> eko = generate_fake_ekos(pids_in, x_in, pids_out, x_in);
-            for (std::size_t k = 0; k != flat_len; k++) {
-                op_slices.push_back(eko[k]);
-            }
-        }
-    }
-
     // Construct the values of alphas table
     std::vector<double> alphas_table;
     for (double q2 : ren1) {
@@ -189,7 +186,7 @@ int main() {
     //     - `ren1`: values of the renormalization scales
     //     - `alphas_table`: values of alphas for each renormalization scales
     pineappl_fk_table* fktable = pineappl_grid_evolve(grid, opinfo_slices.data(),
-        max_orders.data(), op_slices.data(), x_in.data(),
+        max_orders.data(), generate_fake_ekos, x_in.data(),
         x_in.data(), pids_in.data(), pids_out.data(),
         tensor_shape.data(), xi.data(), ren1.data(), alphas_table.data());
 

examples/cpp/evolve-grid.cpp

@@ -82,16 +82,36 @@ std::vector<std::size_t> unravel_index(std::size_t flat_index, const std::vector
     return coords;
 }
 
-std::vector<double> generate_fake_ekos(std::string filename) {
-    std::ifstream input_file(filename);
-    std::vector<double> ops;
+extern "C" void generate_fake_ekos(
+    const int* pids_in,
+    const double* x_in,
+    const int* pids_out,
+    const double* x_out,
+    double* eko_buffer,
+    pineappl_conv_type conv_type,
+    std::size_t pids_in_len,
+    std::size_t x_in_len,
+    std::size_t pids_out_len,
+    std::size_t x_out_len
+) {
+    // Ignore unused variables
+    (void)pids_in;
+    (void)x_in;
+    (void)pids_out;
+    (void)x_out;
+    (void) conv_type;
+    (void)pids_in_len;
+    (void)x_in_len;
+    (void)pids_out_len;
+    (void)x_out_len;
+
+    std::ifstream input_file("../../test-data/EKO_LHCB_WP_7TEV.txt");
     double weight_value;
 
+    std::size_t count = 0;
     while (input_file >> weight_value) {
-        ops.push_back(weight_value);
+        eko_buffer[count++] = weight_value;
     }
-
-    return ops;
 }
 
 void print_results(std::vector<double> dxsec_grid, std::vector<double> dxsec_fktable) {
@@ -126,7 +146,6 @@ void print_results(std::vector<double> dxsec_grid, std::vector<double> dxsec_fkt
 int main() {
     // TODO: How to get a Grid that can be evolved??
     std::string filename = "../../test-data/LHCB_WP_7TEV_opt.pineappl.lz4";
-    std::string ekoname = "../../test-data/EKO_LHCB_WP_7TEV.txt";
 
     // disable LHAPDF banners to guarantee deterministic output
     LHAPDF::setVerbosity(0);
@@ -200,18 +219,6 @@ int main() {
     std::vector<int> pids_in = PIDS;
     std::vector<int> pids_out = PIDS;
 
-    // The Evolution Operator is a vector with length `N_conv * N_Q2_slices * Σ product(OP shape)`
-    std::vector<double> op_slices;
-    std::size_t flat_len = XGRID.size() * XGRID.size() * pids_in.size() * pids_out.size();
-    for (std::size_t _i = 0; _i != conv_types.size(); _i++) {
-        for (std::size_t j = 0; j != fac1.size(); j++) {
-            std::vector<double> eko = generate_fake_ekos(ekoname);
-            for (std::size_t k = 0; k != flat_len; k++) {
-                op_slices.push_back(eko[k]);
-            }
-        }
-    }
-
     // Construct the values of alphas table
     std::vector<double> alphas_table;
     for (double q2 : ren1) {
@@ -237,7 +244,7 @@ int main() {
     //     - `ren1`: values of the renormalization scales
     //     - `alphas_table`: values of alphas for each renormalization scales
     pineappl_fk_table* fktable = pineappl_grid_evolve(grid, opinfo_slices.data(),
-        max_orders.data(), op_slices.data(), XGRID.data(),
+        max_orders.data(), generate_fake_ekos, XGRID.data(),
         XGRID.data(), pids_in.data(), pids_out.data(),
         tensor_shape.data(), xi.data(), ren1.data(), alphas_table.data());
 

pineappl_capi/src/lib.rs

@@ -2132,9 +2132,20 @@ pub unsafe extern "C" fn pineappl_grid_evolve_info(
 #[no_mangle]
 pub unsafe extern "C" fn pineappl_grid_evolve(
     grid: *mut Grid,
-    op_info: *mut OperatorInfo,
+    operator_infos: *mut OperatorInfo,
     max_orders: *const u8,
-    operators: *mut f64,
+    operators: extern "C" fn(
+        *const i32,
+        *const f64,
+        *const i32,
+        *const f64,
+        *mut f64,
+        ConvType,
+        usize,
+        usize,
+        usize,
+        usize,
+    ),
     x_in: *const f64,
     x_out: *const f64,
     pids_in: *const i32,
@@ -2167,46 +2178,50 @@ pub unsafe extern "C" fn pineappl_grid_evolve(
         .map(pineappl::convolutions::Conv::conv_type)
         .collect();
 
-    let op_info = unsafe {
-        slice::from_raw_parts(op_info, conv_types.len() * evolve_info.fac1.len())
+    let operator_infos = unsafe {
+        slice::from_raw_parts(operator_infos, conv_types.len() * evolve_info.fac1.len())
             .chunks_exact(evolve_info.fac1.len())
     };
-
     let total_shape: usize = eko_shape.iter().product();
-    let operators = unsafe {
-        slice::from_raw_parts(
-            operators,
-            conv_types.len() * evolve_info.fac1.len() * total_shape,
-        )
-        .chunks_exact(evolve_info.fac1.len() * total_shape)
-    };
 
-    let slices = op_info
-        .zip(operators)
-        .map(|(op_infos, op_vals)| {
-            op_infos
-                .iter()
-                .zip(op_vals.chunks_exact(total_shape))
-                .map(|(op_info, values)| {
-                    let operator_slice_info = OperatorSliceInfo {
-                        pid_basis: op_info.pid_basis,
-                        fac0: op_info.fac0,
-                        pids0: pids_out.to_vec(),
-                        x0: x_out.to_vec(),
-                        fac1: op_info.fac1,
-                        pids1: pids_in.to_vec(),
-                        x1: x_in.to_vec(),
-                        conv_type: op_info.conv_type,
-                    };
-
-                    let array = Array4::from_shape_vec(
-                        Ix4(eko_shape[0], eko_shape[1], eko_shape[2], eko_shape[3]),
-                        values.to_vec(),
-                    )
-                    .unwrap();
-
-                    Ok::<_, std::io::Error>((operator_slice_info, CowArray::from(array)))
-                })
+    let slices = operator_infos
+        .map(|op_infos| {
+            op_infos.iter().map(|op_info| {
+                let operator_slice_info = OperatorSliceInfo {
+                    pid_basis: op_info.pid_basis,
+                    fac0: op_info.fac0,
+                    pids0: pids_out.to_vec(),
+                    x0: x_out.to_vec(),
+                    fac1: op_info.fac1,
+                    pids1: pids_in.to_vec(),
+                    x1: x_in.to_vec(),
+                    conv_type: op_info.conv_type,
+                };
+
+                // Let Rust manages the memory by allocating a flat vector of zeros as a placeholder
+                // to hold the evolution operator data.
+                let mut eko_slice = vec![0.0; total_shape];
+                operators(
+                    pids_in.as_ptr(),
+                    x_in.as_ptr(),
+                    pids_out.as_ptr(),
+                    x_out.as_ptr(),
+                    eko_slice.as_mut_ptr(),
+                    op_info.conv_type,
+                    pids_in.len(),
+                    x_in.len(),
+                    pids_out.len(),
+                    x_out.len(),
+                );
+
+                let array = Array4::from_shape_vec(
+                    Ix4(eko_shape[0], eko_shape[1], eko_shape[2], eko_shape[3]),
+                    eko_slice,
+                )
+                .unwrap();
+
+                Ok::<_, std::io::Error>((operator_slice_info, CowArray::from(array)))
+            })
         })
         .collect();