Add fastNLO exporter

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

@@ -12,6 +12,9 @@ use std::process::ExitCode;
 #[cfg(feature = "applgrid")]
 mod applgrid;
 
+#[cfg(feature = "fastnlo")]
+mod fastnlo;
+
 #[cfg(feature = "applgrid")]
 fn convert_into_applgrid(
     output: &Path,
@@ -42,6 +45,36 @@ fn convert_into_applgrid(
     ))
 }
 
+#[cfg(feature = "fastnlo")]
+fn convert_into_fastnlo(
+    output: &Path,
+    grid: &Grid,
+    conv_funs: &mut [Pdf],
+    _: usize,
+    discard_non_matching_scales: bool,
+) -> Result<(&'static str, Vec<f64>, usize, Vec<bool>)> {
+    // TODO: check also scale-varied results
+
+    let (mut fastnlo, order_mask) =
+        fastnlo::convert_into_fastnlo(grid, output, discard_non_matching_scales)?;
+    let results = fastnlo::convolve_fastnlo(fastnlo.pin_mut(), conv_funs);
+
+    Ok(("fastNLO", results, 1, order_mask))
+}
+
+#[cfg(not(feature = "fastnlo"))]
+fn convert_into_fastnlo(
+    _: &Path,
+    _: &Grid,
+    _: &mut [Pdf],
+    _: usize,
+    _: bool,
+) -> Result<(&'static str, Vec<f64>, usize, Vec<bool>)> {
+    Err(anyhow!(
+        "you need to install `pineappl` with feature `fastnlo`"
+    ))
+}
+
 fn convert_into_grid(
     output: &Path,
     grid: &Grid,
@@ -58,13 +91,27 @@ fn convert_into_grid(
                 scales,
                 discard_non_matching_scales,
             );
+        } else if extension == "tab"
+            || (extension == "gz"
+                && output
+                    .with_extension("")
+                    .extension()
+                    .map_or(false, |ext| ext == "tab"))
+        {
+            return convert_into_fastnlo(
+                output,
+                grid,
+                conv_funs,
+                scales,
+                discard_non_matching_scales,
+            );
         }
     }
 
     Err(anyhow!("could not detect file format"))
 }
 
-/// Converts PineAPPL grids to APPLgrid files.
+/// Converts PineAPPL grids to APPLgrid/fastNLO files.
 #[derive(Parser)]
 pub struct Opts {
     /// Path to the input grid.

pineappl_cli/src/export/fastnlo.rs

@@ -0,0 +1,113 @@
+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::{self, fastNLOLHAPDF};
+use std::path::Path;
+use std::pin::Pin;
+
+pub fn convert_into_fastnlo(
+    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!()
+}
+
+pub fn convolve_fastnlo(_grid: Pin<&mut fastNLOLHAPDF>, conv_funs: &mut [Pdf]) -> Vec<f64> {
+    // TODO: add support for convolving an APPLgrid with two functions
+    assert_eq!(conv_funs.len(), 1);
+
+    todo!()
+}

pineappl_cli/tests/export.rs

@@ -3,7 +3,7 @@ use assert_cmd::Command;
 #[cfg(feature = "applgrid")]
 use assert_fs::NamedTempFile;
 
-const HELP_STR: &str = "Converts PineAPPL grids to APPLgrid files
+const HELP_STR: &str = "Converts PineAPPL grids to APPLgrid/fastNLO files
 
 Usage: pineappl export [OPTIONS] <INPUT> <OUTPUT> <CONV_FUNS>
 

pineappl_cli/tests/main.rs

@@ -10,7 +10,7 @@ Commands:
   convolve  Convolutes a PineAPPL grid with a PDF set
   diff      Compares the numerical content of two grids with each other
   evolve    Evolve a grid with an evolution kernel operator to an FK table
-  export    Converts PineAPPL grids to APPLgrid files
+  export    Converts PineAPPL grids to APPLgrid/fastNLO files
   help      Display a manpage for selected subcommands
   import    Converts APPLgrid/fastNLO/FastKernel files to PineAPPL grids
   merge     Merges one or more PineAPPL grids together

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++11") // 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));