Allow to specify multiple convolution functions with CLI

CHANGELOG.md

@@ -39,6 +39,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - changed member `lumi_id_types` of `OperatorInfo` and `OperatorSliceInfo` to
   `pid_basis`, which is now of type `PidBasis`
 - renamed module `pineappl::lumi` to `pineappl::convolutions`
+- renamed switch `--pdf` to `--conv-fun` in the subcommand `uncert`. This
+  switch now optionally accepts a list of indices, which determines the
+  corresponding convolution function (PDF/FF), for which the uncertainty should
+  calculated
 
 ### Removed
 

pineappl_cli/src/uncert.rs

@@ -1,22 +1,28 @@
-use super::helpers::{self, ConvoluteMode};
+use super::helpers::{self, ConvFuns, ConvoluteMode};
 use super::{GlobalConfiguration, Subcommand};
-use anyhow::Result;
+use anyhow::{Error, Result};
 use clap::builder::{PossibleValuesParser, TypedValueParser};
 use clap::{Args, Parser, ValueHint};
 use prettytable::{cell, Row};
 use rayon::{prelude::*, ThreadPoolBuilder};
 use std::num::NonZeroUsize;
 use std::path::PathBuf;
 use std::process::ExitCode;
-use std::slice;
 use std::thread;
 
 #[derive(Args)]
 #[group(multiple = true, required = true)]
 struct Group {
-    /// Calculate the PDF uncertainties.
-    #[arg(long)]
-    pdf: bool,
+    /// Calculate convolution function uncertainties.
+    #[arg(
+        default_missing_value = "0",
+        num_args = 0..=1,
+        long,
+        require_equals = true,
+        value_delimiter = ',',
+        value_name = "IDX"
+    )]
+    conv_fun: Vec<usize>,
     /// Show absolute numbers of the scale-varied results.
     #[arg(
         default_missing_value = "7",
@@ -49,18 +55,17 @@ struct Group {
     scale_env: Option<u16>,
 }
 
-/// Calculates scale and PDF uncertainties.
+/// Calculates scale and convolution function uncertainties.
 #[derive(Parser)]
 pub struct Opts {
     /// Path to the input grid.
     #[arg(value_hint = ValueHint::FilePath)]
     input: PathBuf,
-    /// LHAPDF id or name of the PDF set.
-    #[arg(value_parser = helpers::parse_pdfset)]
-    pdfset: String,
+    /// LHAPDF ID(s) or name(s) of the PDF(s)/FF(s).
+    conv_funs: ConvFuns,
     #[command(flatten)]
     group: Group,
-    /// Confidence level in per cent, for PDF uncertainties.
+    /// Confidence level in per cent, for convolution function uncertainties.
     #[arg(default_value_t = lhapdf::CL_1_SIGMA, long)]
     cl: f64,
     /// Show integrated numbers (without bin widths) instead of differential ones.
@@ -89,7 +94,7 @@ pub struct Opts {
 impl Subcommand for Opts {
     fn run(&self, cfg: &GlobalConfiguration) -> Result<ExitCode> {
         let grid = helpers::read_grid(&self.input)?;
-        let (set, _) = helpers::create_pdfset(&self.pdfset)?;
+        let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?;
 
         let limits = helpers::convolve_limits(
             &grid,
@@ -101,34 +106,53 @@ impl Subcommand for Opts {
             },
         );
 
-        let pdf_results = if self.group.pdf {
-            ThreadPoolBuilder::new()
-                .num_threads(self.threads)
-                .build_global()
-                .unwrap();
+        ThreadPoolBuilder::new()
+            .num_threads(self.threads)
+            .build_global()
+            .unwrap();
+
+        let conv_fun_results: Vec<Vec<_>> = self
+            .group
+            .conv_fun
+            .iter()
+            .map(|&index| {
+                let set = conv_funs[index].set();
+                let results: Vec<_> = set
+                    .mk_pdfs()?
+                    .into_par_iter()
+                    .map(|fun| {
+                        // TODO: do not create objects that are getting overwritten in any case
+                        let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?;
+                        conv_funs[index] = fun;
 
-            set.mk_pdfs()?
-                .into_par_iter()
-                .flat_map(|mut pdf| {
-                    helpers::convolve(
-                        &grid,
-                        slice::from_mut(&mut pdf),
-                        &self.orders,
-                        &[],
-                        &[],
-                        1,
-                        if self.integrated {
-                            ConvoluteMode::Integrated
-                        } else {
-                            ConvoluteMode::Normal
-                        },
-                        cfg,
-                    )
-                })
-                .collect()
-        } else {
-            vec![]
-        };
+                        Ok::<_, Error>(helpers::convolve(
+                            &grid,
+                            &mut conv_funs,
+                            &self.orders,
+                            &[],
+                            &[],
+                            1,
+                            if self.integrated {
+                                ConvoluteMode::Integrated
+                            } else {
+                                ConvoluteMode::Normal
+                            },
+                            cfg,
+                        ))
+                    })
+                    .collect::<Result<_, _>>()?;
+
+                // transpose results
+                let results: Vec<Vec<_>> = (0..results[0].len())
+                    .map(|bin| (0..results.len()).map(|pdf| results[pdf][bin]).collect())
+                    .collect();
+
+                results
+                    .into_iter()
+                    .map(|values| Ok(set.uncertainty(&values, self.cl, false)?))
+                    .collect::<Result<_, _>>()
+            })
+            .collect::<Result<_, Error>>()?;
         let scales_max = self
             .group
             .scale_env
@@ -140,7 +164,7 @@ impl Subcommand for Opts {
             .unwrap_or(1);
         let scale_results = helpers::convolve(
             &grid,
-            slice::from_mut(&mut helpers::create_pdf(&self.pdfset)?),
+            &mut conv_funs,
             &self.orders,
             &[],
             &[],
@@ -163,9 +187,10 @@ impl Subcommand for Opts {
         }
         title.add_cell(cell!(c->format!("{y_label}\n[{y_unit}]")));
 
-        if self.group.pdf {
-            title.add_cell(cell!(c->"PDF central"));
-            title.add_cell(cell!(c->"PDF\n[%]").with_hspan(2));
+        for &index in &self.group.conv_fun {
+            title.add_cell(
+                cell!(c->format!("{}", self.conv_funs.lhapdf_names[index])).with_hspan(3),
+            );
         }
 
         if let Some(scales) = self.group.scale_abs {
@@ -190,24 +215,6 @@ impl Subcommand for Opts {
             .zip(scale_results.chunks_exact(scales_max))
             .enumerate()
         {
-            let (pdf_cen, pdf_neg, pdf_pos) = if self.group.pdf {
-                let values: Vec<_> = pdf_results
-                    .iter()
-                    .skip(bin)
-                    .step_by(limits.len())
-                    .copied()
-                    .collect();
-                let uncertainty = set.uncertainty(&values, self.cl, false)?;
-
-                (
-                    uncertainty.central,
-                    -100.0 * uncertainty.errminus / uncertainty.central,
-                    100.0 * uncertainty.errplus / uncertainty.central,
-                )
-            } else {
-                (0.0, 0.0, 0.0)
-            };
-
             let row = table.add_empty_row();
             row.add_cell(cell!(r->format!("{bin}")));
             for (left, right) in left_right_limits {
@@ -217,10 +224,10 @@ impl Subcommand for Opts {
 
             row.add_cell(cell!(r->format!("{:.*e}", self.digits_abs, scale_res[0])));
 
-            if self.group.pdf {
-                row.add_cell(cell!(r->format!("{:.*e}", self.digits_abs, pdf_cen)));
-                row.add_cell(cell!(r->format!("{:.*}", self.digits_rel, pdf_neg)));
-                row.add_cell(cell!(r->format!("{:.*}", self.digits_rel, pdf_pos)));
+            for uncertainty in conv_fun_results.iter().map(|results| &results[bin]) {
+                row.add_cell(cell!(r->format!("{:.*e}", self.digits_abs, uncertainty.central)));
+                row.add_cell(cell!(r->format!("{:.*}", self.digits_rel, -100.0 * uncertainty.errminus / uncertainty.central)));
+                row.add_cell(cell!(r->format!("{:.*}", self.digits_rel, 100.0 * uncertainty.errplus / uncertainty.central)));
             }
 
             if let Some(scales) = self.group.scale_abs {