Add code to initialize an new fastNLO table

pineappl_cli/Cargo.toml

@@ -54,6 +54,6 @@ rustc-args = [ "--cfg feature=\"docs-only\"" ]
 [features]
 applgrid = ["dep:cxx", "dep:pineappl_applgrid"]
 evolve = ["dep:base64", "dep:either", "dep:tar", "dep:lz4_flex", "dep:ndarray-npy", "dep:serde", "dep:serde_yaml"]
-fastnlo = ["dep:pineappl_fastnlo"]
+fastnlo = ["dep:cxx", "dep:pineappl_fastnlo"]
 fktable = ["dep:flate2",  "dep:tar"]
 static = ["lhapdf/static", "pineappl/static", "pineappl_applgrid?/static", "pineappl_fastnlo?/static"]

pineappl_cli/src/export/fastnlo.rs

@@ -1,16 +1,107 @@
-use anyhow::Result;
+use anyhow::{bail, Result};
 use cxx::UniquePtr;
+use float_cmp::assert_approx_eq;
 use lhapdf::Pdf;
+use pineappl::boc::Order;
 use pineappl::grid::Grid;
-use pineappl_fastnlo::ffi::fastNLOLHAPDF;
+use pineappl_fastnlo::ffi::{self, fastNLOLHAPDF};
 use std::path::Path;
 use std::pin::Pin;
 
 pub fn convert_into_fastnlo(
-    _grid: &Grid,
+    grid: &Grid,
     _output: &Path,
     _discard_non_matching_scales: bool,
 ) -> Result<(UniquePtr<fastNLOLHAPDF>, Vec<bool>)> {
+    let bin_info = grid.bin_info();
+    let dim = bin_info.dimensions();
+
+    if dim > 3 {
+        bail!(
+            "grid has {} dimensions, but fastNLO only supports up to three-dimensional distributions",
+            dim
+        );
+    }
+
+    let bin_limits = bin_info.limits();
+    let left_bin_limits: Vec<Vec<_>> = bin_limits
+        .iter()
+        .map(|limits| limits.iter().map(|&(left, _)| left).collect())
+        .collect();
+    let right_bin_limits: Vec<Vec<_>> = bin_limits
+        .iter()
+        .map(|limits| limits.iter().map(|&(_, right)| right).collect())
+        .collect();
+    let normalizations = bin_info.normalizations();
+
+    let order_mask = Order::create_mask(grid.orders(), 3, 0, false);
+    let orders_with_mask: Vec<_> = grid
+        .orders()
+        .iter()
+        .cloned()
+        .zip(order_mask.iter().copied())
+        .collect();
+    let lo_alphas = orders_with_mask
+        .iter()
+        .filter_map(|&(Order { alphas, .. }, keep)| keep.then_some(alphas))
+        .min()
+        // UNWRAP: this will fail for `Grid` with no orders, but this shouldn't happen
+        .unwrap();
+    //let loops = orders_with_mask
+    //    .iter()
+    //    .filter_map(|&(Order { alphas, .. }, keep)| keep.then_some(alphas))
+    //    .max()
+    //    .unwrap()
+    //    - lo_alphas;
+
+    let convolutions: Vec<i32> = grid
+        .convolutions()
+        .iter()
+        .filter_map(|conv| conv.pid())
+        .collect();
+
+    let channels: Vec<Vec<_>> = grid
+        .channels()
+        .iter()
+        .map(|channel| {
+            channel
+                .entry()
+                .iter()
+                .map(|&(a, b, factor)| {
+                    assert_approx_eq!(f64, factor, 1.0, ulps = 4);
+                    ffi::pair_int_int {
+                        first: a,
+                        second: b,
+                    }
+                })
+                .collect()
+        })
+        .collect();
+
+    //for (fnlo_order, order) in order_mask
+    //    .iter()
+    //    .enumerate()
+    //    .filter_map(|(index, keep)| keep.then_some(index))
+    //    .enumerate()
+    //{}
+
+    let _fastnlo = ffi::make_fastnlo_create(
+        // UNWRAP: negative numbers and overflow should not happen
+        lo_alphas.try_into().unwrap(),
+        &left_bin_limits,
+        &right_bin_limits,
+        &normalizations,
+        // TODO: calculate channels for each order separately
+        // UNWRAP: negative numbers and overflow should not happen
+        channels.len().try_into().unwrap(),
+        // UNWRAP: negative numbers and overflow should not happen
+        channels.len().try_into().unwrap(),
+        // UNWRAP: negative numbers and overflow should not happen
+        channels.len().try_into().unwrap(),
+        &convolutions,
+        &channels,
+    );
+
     todo!()
 }
 

pineappl_fastnlo/build.rs

@@ -65,7 +65,7 @@ fn main() {
         .file("src/fastnlo.cpp")
         .include(fnlo_include_path.trim())
         .includes(lhapdf_include_paths)
-        .std("c++17") // apparently not supported by MSVC, but fastNLO probably can't be compiled on Windows
+        .std("c++17")
         .compile("fnlo-bridge");
 
     println!("cargo:rerun-if-changed=src/lib.rs");

pineappl_fastnlo/src/fastnlo.cpp

@@ -1,6 +1,7 @@
 #include "pineappl_fastnlo/src/fastnlo.hpp"
 
 #include <algorithm>
+#include <cassert>
 #include <iterator>
 #include <string>
 
@@ -65,6 +66,156 @@ std::unique_ptr<fastNLOLHAPDF> make_fastnlo_lhapdf_with_name_file_set(
     return std::unique_ptr<fastNLOLHAPDF>(new fastNLOLHAPDF(arg0, arg1, PDFSet));
 }
 
+std::unique_ptr<fastNLOCreate> make_fastnlo_create(
+    int alphas_lo,
+    rust::Slice<rust::Vec<double> const> left_bin_limits,
+    rust::Slice<rust::Vec<double> const> right_bin_limits,
+    rust::Slice<double const> normalizations,
+    int lo_channels,
+    int nlo_channels,
+    int nnlo_channels,
+    rust::Slice<int const> convolutions,
+    rust::Slice<rust::Vec<pair_int_int> const> channels
+) {
+    assert(left_bin_limits.size() == right_bin_limits.size());
+    auto const bins = left_bin_limits.size();
+    assert(bins == normalizations.size());
+    assert(bins > 0);
+    auto const dimensions = left_bin_limits.at(0).size();
+    assert(dimensions > 0);
+    assert(convolutions.size() <= 2);
+    assert(convolutions.size() >= 1);
+
+    std::vector<std::vector<double>> bin_limits(dimensions);
+
+    // TODO: check if this is the right ordering
+    for (std::size_t i = 0; i != dimensions; ++i) {
+        assert(left_bin_limits.at(i).size() == dimensions);
+        assert(right_bin_limits.at(i).size() == dimensions);
+
+        //bin_limits.at(i).resize(2 * limits);
+
+        //for (std::size_t j = 0; j != limits; ++j) {
+        //    bin_limits.at(i).at(2 * j + 0) = left_bin_limits.at(j).at(i);
+        //    bin_limits.at(i).at(2 * j + 1) = right_bin_limits.at(j).at(i);
+        //}
+        bin_limits.at(i).resize(bins + 1);
+        bin_limits.at(i).at(0) = left_bin_limits.at(0).front();
+
+        for (std::size_t j = 0; j != bins; ++j) {
+            bin_limits.at(i).at(j + 1) = right_bin_limits.at(j).at(i);
+        }
+    }
+
+    fastNLO::GeneratorConstants gconst;
+    // TODO: add PineAPPL's version number
+    gconst.Name = "PineAPPL-fastNLO interface";
+
+    fastNLO::ProcessConstants pconst;
+    pconst.LeadingOrder = alphas_lo;
+    pconst.NPDF = convolutions.size();
+    pconst.NSubProcessesLO = lo_channels;
+    pconst.NSubProcessesNLO = nlo_channels;
+    pconst.NSubProcessesNNLO = nnlo_channels;
+
+    if (convolutions.size() == 1) {
+        pconst.IPDFdef1 = 2;
+    } else {
+        pconst.IPDFdef1 = 3;
+    }
+
+    // TODO: is this the correct value to set the linear combinations ourselves?
+    pconst.IPDFdef2 = 0;
+    pconst.IPDFdef3LO = 2;
+    pconst.IPDFdef3NLO = 0;
+    pconst.IPDFdef3NNLO = 0;
+
+    if (convolutions.size() == 1) {
+        // TODO: not yet implemented
+        assert(false);
+    } else {
+        pconst.NPDFDim = 2;
+    }
+
+    std::vector<std::vector<std::pair<int, int>>> linear_combinations(channels.size());
+    for (std::size_t i = 0; i != channels.size(); ++i) {
+        std::vector<std::pair<int, int>> entries(channels.at(i).size());
+        for (std::size_t j = 0; j != channels.at(i).size(); ++j) {
+            auto const first = channels.at(i).at(j).first;
+            auto const second = channels.at(i).at(j).second;
+            entries.at(j) = std::make_pair(first, second);
+        }
+        linear_combinations.at(i) = entries;
+    }
+    pconst.PDFCoeffLO = linear_combinations;
+
+    fastNLO::ScenarioConstants sconst;
+    sconst.DifferentialDimension = dimensions;
+    sconst.DimensionIsDifferential = std::vector<int>(dimensions, 0);
+    sconst.CalculateBinSize = false;
+    sconst.BinSize = std::vector<double>(normalizations.begin(), normalizations.end());
+
+    switch (sconst.DifferentialDimension) {
+        case 1:
+            sconst.SingleDifferentialBinning = bin_limits.at(0);
+            break;
+
+        case 2:
+            sconst.DoubleDifferentialBinning = bin_limits;
+            break;
+
+        case 3:
+            sconst.TripleDifferentialBinning = bin_limits;
+            break;
+
+        default:
+            // ASSERT: there are no or too many dimensions, which fastNLO doesn't support
+            assert(false);
+    }
+    sconst.FlexibleScaleTable = true;
+
+    if (convolutions.size() == 1) {
+        sconst.PDF1 = convolutions.at(0);
+        // TODO: do we leave PDF2 unchanged (set to 'proton') for DIS?
+    } else {
+        sconst.PDF1 = convolutions.at(0);
+        sconst.PDF2 = convolutions.at(1);
+    }
+
+    sconst.ReadBinningFromSteering = true;
+    sconst.IgnoreWarmupBinningCheck = true;
+    sconst.X_NNodeCounting = "NodesPerBin";
+    sconst.Mu1_NNodeCounting = "NodesPerBin";
+    sconst.Mu2_NNodeCounting = "NodesPerBin";
+
+    fastNLO::WarmupConstants wconst(sconst);
+    wconst.Binning = bin_limits;
+
+    // these values are probably irrelevant but must nevertheless given
+    wconst.Values.resize(bins, std::vector<double>{
+        // bin index
+        0,
+        // x-min
+        2e-7,
+        // x-max
+        1.0,
+        // scale1-min
+        10.0,
+        // scale1-max
+        100.0,
+        // scale2-min
+        10.0,
+        // scale2-max
+        100.0
+    });
+    for (std::size_t i = 0; i != wconst.Values.size(); ++i) {
+        wconst.Values.at(i).at(0) = static_cast <double> (i);
+    }
+    // wconst.headerValues = ;
+
+    return std::unique_ptr<fastNLOCreate>(new fastNLOCreate(gconst, pconst, sconst, wconst));
+}
+
 rust::Vec<double> GetCrossSection(fastNLOReader& reader, bool lNorm)
 {
     return std_vector_to_rust_vec(reader.GetCrossSection(lNorm));

pineappl_fastnlo/src/fastnlo.hpp

@@ -8,6 +8,7 @@
 #include <fastnlotk/fastNLOCoeffAddBase.h>
 #include <fastnlotk/fastNLOCoeffAddFix.h>
 #include <fastnlotk/fastNLOCoeffAddFlex.h>
+#include <fastnlotk/fastNLOCreate.h>
 #include <fastnlotk/fastNLOLHAPDF.h>
 #include <fastnlotk/fastNLOPDFLinearCombinations.h>
 #include <fastnlotk/fastNLOReader.h>
@@ -21,6 +22,18 @@ std::unique_ptr<fastNLOLHAPDF> make_fastnlo_lhapdf_with_name_file_set(
     int PDFSet
 );
 
+std::unique_ptr<fastNLOCreate> make_fastnlo_create(
+    int alphas_lo,
+    rust::Slice<rust::Vec<double> const> left_bin_limits,
+    rust::Slice<rust::Vec<double> const> right_bin_limits,
+    rust::Slice<double const> normalizations,
+    int lo_channels,
+    int nlo_channels,
+    int nnlo_channels,
+    rust::Slice<int const> convolutions,
+    rust::Slice<rust::Vec<pair_int_int> const> channels
+);
+
 rust::Vec<double> CalcPDFLinearCombination(
     fastNLOPDFLinearCombinations const& lc,
     fastNLOCoeffAddBase const& base,

pineappl_fastnlo/src/lib.rs

@@ -134,6 +134,12 @@ pub mod ffi {
         fn GetPDFPDG(&self, _: i32) -> i32;
     }
 
+    unsafe extern "C++" {
+        include!("fastnlotk/fastNLOCreate.h");
+
+        type fastNLOCreate;
+    }
+
     unsafe extern "C++" {
         include!("pineappl_fastnlo/src/fastnlo.hpp");
 
@@ -189,6 +195,18 @@ pub mod ffi {
             _: &str,
             _: i32,
         ) -> UniquePtr<fastNLOLHAPDF>;
+
+        fn make_fastnlo_create(
+            alphas_lo: i32,
+            left_bin_limits: &[Vec<f64>],
+            right_bin_limits: &[Vec<f64>],
+            normalizations: &[f64],
+            lo_channels: i32,
+            nlo_channels: i32,
+            nnlo_channels: i32,
+            convolutions: &[i32],
+            channels: &[Vec<pair_int_int>],
+        ) -> UniquePtr<fastNLOCreate>;
     }
 }