Refine combination of LHAPDF sets logics

.github/actions/cache-data/action.yml

@@ -10,7 +10,7 @@ runs:
       uses: actions/cache@v4
       with:
         path: neopdf-data
-        key: data-v4
+        key: data-v5
     - name: Download data if cache miss
       if: steps.cache-data.outputs.cache-hit != 'true'
       run: |

docs/cli-tutorials.md

@@ -118,9 +118,9 @@ neopdf compute alphas_q2 --pdf-name NNPDF40_nnlo_as_01180 --member 0 --q2 10
 
 ---
 
-## Converting LHAPDF into NeoPDF Format
+## Converting LHAPDF Sets into NeoPDF Format
 
-### Convert LHAPDF set to NeoPDF format
+### Convert standard LHAPDF to NeoPDF format
 
 To convert a standard LHAPDF set into the new, compressed `NeoPDF` format:
 
@@ -135,7 +135,7 @@ neopdf write convert --pdf-name NNPDF40_nnlo_as_01180 --output NNPDF40_nnlo_as_0
 To combine several LHAPDF nuclear PDF sets into a single `NeoPDF` grid:
 
 ```bash
-neopdf write combine --pdf-names nNNPDF30_nlo_as_0118_p nNNPDF30_nlo_as_0118_A12 nNNPDF30_nlo_as_0118_A40 --output nNNPDF30_nlo_as_0118.neopdf.lz4
+neopdf write combine-npdfs --pdf-names nNNPDF30_nlo_as_0118_p nNNPDF30_nlo_as_0118_A12 nNNPDF30_nlo_as_0118_A40 --output nNNPDF30_nlo_as_0118.neopdf.lz4
 ```
 
 - `--pdf-names`: List of PDF set names to combine.
@@ -170,7 +170,7 @@ nNNPDF30_nlo_as_0118_A208_Z82
 the command now becomes:
 
 ```bash
-neopdf write combine --names-file nuclearpdfs.txt --output nNNPDF30_nlo_as_0118.neopdf.lz4
+neopdf write combine-npdfs --names-file nuclearpdfs.txt --output nNNPDF30_nlo_as_0118.neopdf.lz4
 ```
 
 This will generate a `nNNPDF30_nlo_as_0118.neopdf.lz4` PDF grid that also contains the $A$
@@ -179,3 +179,50 @@ dependence. One can now check the value of $xf(A, x, Q^2)$ for the Iron $^{56}_{
 ```bash
 neopdf compute xfx_q2 --pdf-name nNNPDF30_nlo_as_0118.neopdf.lz4 --member 0 --pid 21 56 1e-3 10.0
 ```
+
+### Combine multiple  $\alpha_s$ LHAPDF sets
+
+To combine multiple $\alpha_s$ LHAPDF sets into a single NeoPDF grid with an explicit dependence
+on $\alpha_s$, the procedure is the same as when combining multiple nuclear PDFs with the option
+`combine-npdfs` replaced with `combine-alphas`. That is, given an input file with the names of the
+PDF sets:
+
+```yaml title="alphaspdfs.txt"
+NNPDF40_nnlo_as_01160
+NNPDF40_nnlo_as_01170
+NNPDF40_nnlo_as_01175
+NNPDF40_nnlo_as_01185
+NNPDF40_nnlo_as_01190
+```
+
+The command to run is:
+
+```bash
+neopdf write combine-alphas --names-file alphaspdfs.txt --output NNPDF40_nnlo.neopdf.lz4
+```
+
+Note that the names of the PDF sets can be passed via the command line using the option
+`--pdf--name`. We can then interpolate the $\alpha_s = 0.1180$ value:
+
+```bash
+neopdf compute xfx_q2 --pdf-name NNPDF40_nnlo.neopdf.lz4 --member 0 --pid 21 0.1180 1e-3 10.0
+```
+
+!!! danger "Warning"
+
+    Given that LHAPDF doesn't provide specific attributes to extract the values of $A$ and $\alpha_s$,
+    their values are inferred from the set name. The extraction of these values are therefore set
+    dependent and we try to support as many sets as possible while keeping the implementation modular:
+
+    ```rust
+    // Regexes to extract A from the PDF set name
+    let re_nnpdf = Regex::new(r"_A(\d+)").unwrap();
+    let re_ncteq = Regex::new(r"_(\d+)_(\d+)$").unwrap();
+    let re_epps = Regex::new(r"[a-zA-Z]+(\d+)$").unwrap();
+
+    // Regexes to extract alpha_s from the PDF set name
+    let re_nnpdf_ct = Regex::new(r"_as_(\d+)_?").unwrap();
+    let re_msht = Regex::new(r"_as(\d+)").unwrap();
+    ```
+
+    This approach is not fullproof and some PDF sets might not be supported at all.

maintainer/download-data.sh

@@ -9,6 +9,11 @@ LHAPDF_SETS=(
   nNNPDF30_nlo_as_0118_A4_Z2
   nNNPDF30_nlo_as_0118_A6_Z3
   nNNPDF30_nlo_as_0118_A9_Z4
+  NNPDF40_nnlo_as_01160
+  NNPDF40_nnlo_as_01170
+  NNPDF40_nnlo_as_01175
+  NNPDF40_nnlo_as_01185
+  NNPDF40_nnlo_as_01190
 )
 
 NEOPDF_SETS=(

neopdf/src/converter.rs

@@ -1,6 +1,6 @@
 //! This module provides utilities for converting LHAPDF sets to the NeoPDF format and for
 //! combining multiple nuclear PDF sets into a single NeoPDF file.
-//!
+//_!
 //! Main functions:
 //! - `convert_lhapdf`: Converts an LHAPDF set to NeoPDF format and writes it to disk.
 //! - `combine_lhapdf_npdfs`: Combines several nuclear PDF sets (with different nucleon
@@ -45,7 +45,8 @@ pub fn convert_lhapdf<P: AsRef<std::path::Path>>(
     Ok(())
 }
 
-/// Combines a list of nuclear PDF sets (differing in nucleon number A) into a single NeoPDF file with explicit A dependence.
+/// Combines a list of nuclear PDF sets (differing in nucleon number A) into a single NeoPDF
+/// file with explicit A dependence.
 ///
 /// # Arguments
 /// * `pdf_names` - List of PDF set names (each with a different A).
@@ -57,20 +58,23 @@ pub fn combine_lhapdf_npdfs<P: AsRef<std::path::Path>>(
     pdf_names: &[&str],
     output_path: P,
 ) -> Result<(), Box<dyn std::error::Error>> {
-    // let pdfs = [ "nNNPDF30_nlo_as_0118_p", ... ];
-    // combine_nuclear_pdfs_with_a_dependence(&pdfs, "output.neopdf.lz4")?;
     if pdf_names.is_empty() {
         return Err("No PDF set names provided".into());
     }
 
-    // Regex to extract A from the PDF set name (e.g., _A40_)
-    let re_a = Regex::new(r"_A(\d+)").unwrap();
+    // Regexes to extract A from the PDF set name
+    let re_nnpdf = Regex::new(r"_A(\d+)").unwrap();
+    let re_ncteq = Regex::new(r"_(\d+)_(\d+)$").unwrap();
+    let re_epps = Regex::new(r"[a-zA-Z]+(\d+)$").unwrap();
     let mut a_values = Vec::new();
     let mut all_members: Vec<Vec<(MetaData, GridArray)>> = Vec::new();
 
     for &pdf_name in pdf_names {
-        // Extract A
-        let a = if let Some(cap) = re_a.captures(pdf_name) {
+        let a = if let Some(cap) = re_nnpdf.captures(pdf_name) {
+            cap[1].parse::<f64>().unwrap()
+        } else if let Some(cap) = re_ncteq.captures(pdf_name) {
+            cap[1].parse::<f64>().unwrap()
+        } else if let Some(cap) = re_epps.captures(pdf_name) {
             cap[1].parse::<f64>().unwrap()
         } else if pdf_name.ends_with("_p") {
             1.0 // proton
@@ -100,10 +104,8 @@ pub fn combine_lhapdf_npdfs<P: AsRef<std::path::Path>>(
     meta.interpolator_type = InterpolatorType::LogTricubic;
 
     for member_idx in 0..num_members {
-        // For each set, get the GridArray for this member
         let member_arrays: Vec<&GridArray> = all_members.iter().map(|v| &v[member_idx].1).collect();
 
-        // Assume all have the same pids and subgrid structure
         let pids = member_arrays[0].pids.clone();
         let num_subgrids = member_arrays[0].subgrids.len();
         if !member_arrays
@@ -164,3 +166,121 @@ pub fn combine_lhapdf_npdfs<P: AsRef<std::path::Path>>(
     GridArrayCollection::compress(&combined_grids, &meta, output_path)?;
     Ok(())
 }
+
+/// Combines a list of PDF sets (differing in alpha_s) into a single NeoPDF file with explicit
+/// `alpha_s` dependence.
+///
+/// # Arguments
+/// * `pdf_names` - List of PDF set names (each with a different alpha_s).
+/// * `output_path` - Output NeoPDF file path.
+///
+/// # Errors
+/// Returns an error if loading or writing fails, or if the sets are not compatible.
+pub fn combine_lhapdf_alphas<P: AsRef<std::path::Path>>(
+    pdf_names: &[&str],
+    output_path: P,
+) -> Result<(), Box<dyn std::error::Error>> {
+    if pdf_names.is_empty() {
+        return Err("No PDF set names provided".into());
+    }
+
+    // Regexes to extract alpha_s from the PDF set name
+    let re_nnpdf_ct = Regex::new(r"_as_(\d+)_?").unwrap();
+    let re_msht = Regex::new(r"_as(\d+)").unwrap();
+    let mut alphas_values = Vec::new();
+    let mut all_members: Vec<Vec<(MetaData, GridArray)>> = Vec::new();
+
+    for &pdf_name in pdf_names {
+        let alphas = if let Some(cap) = re_nnpdf_ct.captures(pdf_name) {
+            cap[1].parse::<f64>().unwrap() / 10000.0
+        } else if let Some(cap) = re_msht.captures(pdf_name) {
+            cap[1].parse::<f64>().unwrap() / 1000.0
+        } else {
+            return Err(format!("Could not extract alpha_s from PDF name: {}", pdf_name).into());
+        };
+        alphas_values.push(alphas);
+        let set = LhapdfSet::new(pdf_name);
+        let members = set.members();
+        if members.is_empty() {
+            return Err(format!("No members found in set: {}", pdf_name).into());
+        }
+        all_members.push(members);
+    }
+
+    // Check all sets have the same number of members
+    let num_members = all_members[0].len();
+    if !all_members.iter().all(|v| v.len() == num_members) {
+        return Err("All sets must have the same number of members".into());
+    }
+
+    // For each member index, combine the corresponding member from each set along the alpha_s dimension
+    let mut combined_grids = Vec::with_capacity(num_members);
+    let mut meta = all_members[0][0].0.clone();
+    meta.set_desc = format!("Combined alpha_s PDFs: {}", pdf_names.join(", "));
+    meta.num_members = num_members as u32;
+    meta.interpolator_type = InterpolatorType::LogTricubic;
+
+    for member_idx in 0..num_members {
+        // For each set, get the GridArray for this member
+        let member_arrays: Vec<&GridArray> = all_members.iter().map(|v| &v[member_idx].1).collect();
+
+        // Assume all have the same pids and subgrid structure
+        let pids = member_arrays[0].pids.clone();
+        let num_subgrids = member_arrays[0].subgrids.len();
+        if !member_arrays
+            .iter()
+            .all(|ga| ga.pids == pids && ga.subgrids.len() == num_subgrids)
+        {
+            return Err("All sets must have the same flavors and subgrid structure".into());
+        }
+
+        // For each subgrid, stack along the alpha_s dimension
+        let mut combined_subgrids = Vec::with_capacity(num_subgrids);
+        for subgrid_idx in 0..num_subgrids {
+            // For each set, get the subgrid
+            let subgrids: Vec<&SubGrid> = member_arrays
+                .iter()
+                .map(|ga| &ga.subgrids[subgrid_idx])
+                .collect();
+
+            // Check x, q2, nucleons shapes match
+            let xs = &subgrids[0].xs;
+            let q2s = &subgrids[0].q2s;
+            let kts = &subgrids[0].kts;
+            let nucleons = &subgrids[0].nucleons;
+            if !subgrids.iter().all(|sg| {
+                sg.xs == *xs && sg.q2s == *q2s && sg.nucleons == *nucleons && sg.kts == *kts
+            }) {
+                return Err("All sets must have the same x, q2, kT, and nucleons grids".into());
+            }
+
+            // Concatenate along the alphas axis to get [..., alphas=pdf_names.len(), ...]
+            let grid_views: Vec<_> = subgrids.iter().map(|sg| sg.grid.view()).collect();
+            let concatenated = concatenate(Axis(1), &grid_views.to_vec())?;
+            let alphas = Array1::from(alphas_values.clone());
+            let new_subgrid = SubGrid {
+                xs: xs.clone(),
+                q2s: q2s.clone(),
+                kts: kts.clone(),
+                grid: concatenated,
+                nucleons: nucleons.clone(),
+                alphas,
+                nucleons_range: subgrids[0].nucleons_range,
+                alphas_range: subgrids[0].alphas_range,
+                kt_range: subgrids[0].kt_range,
+                x_range: subgrids[0].x_range,
+                q2_range: subgrids[0].q2_range,
+            };
+            combined_subgrids.push(new_subgrid);
+        }
+        let combined_grid = GridArray {
+            pids: pids.clone(),
+            subgrids: combined_subgrids,
+        };
+        combined_grids.push(combined_grid);
+    }
+
+    let combined_grids: Vec<&GridArray> = combined_grids.iter().collect();
+    GridArrayCollection::compress(&combined_grids, &meta, output_path)?;
+    Ok(())
+}

neopdf_cli/src/converter.rs

@@ -1,5 +1,5 @@
 //! CLI logic for `NeoPDF` conversion utilities.
-//!
+//_!
 //! This module defines the command-line interface for converting LHAPDF sets to `NeoPDF` format
 //! and for combining multiple nuclear PDFs into a single `NeoPDF` file with explicit A dependence.
 
@@ -30,7 +30,7 @@ pub enum Commands {
         output: String,
     },
     /// Combine multiple nuclear PDFs into a single `NeoPDF` with A dependence.
-    Combine {
+    CombineNpdfs {
         /// List of PDF set names (each with a different A).
         #[arg(
             short = 'n',
@@ -45,6 +45,22 @@ pub enum Commands {
         #[arg(short, long)]
         output: String,
     },
+    /// Combine multiple PDFs with different `alpha_s` values into a single `NeoPDF`.
+    CombineAlphas {
+        /// List of PDF set names (each with a different `alpha_s`).
+        #[arg(
+            short = 'n',
+            long = "pdf-names",
+            required_unless_present = "names_file"
+        )]
+        pdf_names: Option<Vec<String>>,
+        /// Path to a file containing PDF set names, one per line.
+        #[arg(short = 'f', long = "names-file", conflicts_with = "pdf_names")]
+        names_file: Option<String>,
+        /// Output path for the combined `NeoPDF` file.
+        #[arg(short, long)]
+        output: String,
+    },
 }
 
 // Loads PDF names from either command line arguments or a file.
@@ -101,7 +117,7 @@ pub fn run_cli(cli: Cli) -> Result<(), Box<dyn std::error::Error>> {
         Commands::Convert { pdf_name, output } => {
             converter::convert_lhapdf(pdf_name, output)?;
         }
-        Commands::Combine {
+        Commands::CombineNpdfs {
             pdf_names,
             names_file,
             output,
@@ -110,6 +126,15 @@ pub fn run_cli(cli: Cli) -> Result<(), Box<dyn std::error::Error>> {
             let names_str: Vec<&str> = names.iter().map(String::as_str).collect();
             converter::combine_lhapdf_npdfs(&names_str, output)?;
         }
+        Commands::CombineAlphas {
+            pdf_names,
+            names_file,
+            output,
+        } => {
+            let names = load_pdf_names(pdf_names.as_deref(), names_file.as_deref())?;
+            let names_str: Vec<&str> = names.iter().map(String::as_str).collect();
+            converter::combine_lhapdf_alphas(&names_str, output)?;
+        }
     }
     Ok(())
 }