Accommodate for two different EKOs in Evolution

.github/workflows/rust.yml

@@ -25,7 +25,7 @@ jobs:
         uses: actions/cache@v4
         with:
           path: test-data
-          key: test-data-v11
+          key: test-data-v12
       - name: Download test data
         if: steps.cache-test-data.outputs.cache-hit != 'true'
         run: |
@@ -55,6 +55,9 @@ jobs:
           curl -s -C - -O 'https://ploughshare.web.cern.ch/ploughshare/db/applfast/applfast-h1-dijets-appl-arxiv-0010054/grids/applfast-h1-dijets-appl-arxiv-0010054-xsec000.appl'
           curl -s -C - -O 'https://ploughshare.web.cern.ch/ploughshare/db/applfast/applfast-h1-incjets-fnlo-arxiv-0706.3722/grids/applfast-h1-incjets-fnlo-arxiv-0706.3722-xsec000.tab.gz'
           curl -s -C - -O 'https://ploughshare.web.cern.ch/ploughshare/db/atlas/atlas-atlas-wpm-arxiv-1109.5141/grids/atlas-atlas-wpm-arxiv-1109.5141-xsec001.appl'
+          curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/STAR_WMWP_510GEV_WM-AL-POL.pineappl.lz4'
+          curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/STAR_WMWP_510GEV_WM-AL-POL_PolPDF.tar'
+          curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/STAR_WMWP_510GEV_WM-AL-POL_UnpolPDF.tar'
 
       - name: Set RUSTDOCFLAGS
         run: |

CHANGELOG.md

@@ -7,6 +7,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- added new method `Grid::evolve_with_slice_iter2` which is able to perform
+  evolutions with two different EKOs
+
 ### Fixed
 
 - fixed CI to build CAPI and CLI

maintainer/generate-coverage.sh

@@ -27,6 +27,9 @@ wget --no-verbose --no-clobber -P test-data 'https://ploughshare.web.cern.ch/plo
 wget --no-verbose --no-clobber -P test-data 'https://ploughshare.web.cern.ch/ploughshare/db/applfast/applfast-h1-dijets-appl-arxiv-0010054/grids/applfast-h1-dijets-appl-arxiv-0010054-xsec000.appl'
 wget --no-verbose --no-clobber -P test-data 'https://ploughshare.web.cern.ch/ploughshare/db/applfast/applfast-h1-incjets-fnlo-arxiv-0706.3722/grids/applfast-h1-incjets-fnlo-arxiv-0706.3722-xsec000.tab.gz'
 wget --no-verbose --no-clobber -P test-data 'https://ploughshare.web.cern.ch/ploughshare/db/atlas/atlas-atlas-wpm-arxiv-1109.5141/grids/atlas-atlas-wpm-arxiv-1109.5141-xsec001.appl'
+wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/STAR_WMWP_510GEV_WM-AL-POL.pineappl.lz4'
+wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/STAR_WMWP_510GEV_WM-AL-POL_PolPDF.tar'
+wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/STAR_WMWP_510GEV_WM-AL-POL_UnpolPDF.tar'
 
 # we compile with different flags and don't want to destroy the other target directory
 export CARGO_TARGET_DIR="$(mktemp -d)"

maintainer/pineappl-ci/script.sh

@@ -42,7 +42,7 @@ ldconfig
 cd ..
 
 # install PDF sets
-for pdf in NNPDF31_nlo_as_0118_luxqed NNPDF40_nnlo_as_01180 NNPDF40_nlo_as_01180; do
+for pdf in NNPDF31_nlo_as_0118_luxqed NNPDF40_nnlo_as_01180 NNPDF40_nlo_as_01180 NNPDF40_nlo_pch_as_01180; do
     curl "https://lhapdfsets.web.cern.ch/current/${pdf}.tar.gz" | tar xzf - -C /usr/local/share/LHAPDF
 done
 

pineappl/src/evolution.rs

@@ -269,7 +269,7 @@ fn operator_slices(
 type X1aX1bOp2Tuple = (Vec<Vec<f64>>, Array2<f64>);
 
 fn ndarray_from_subgrid_orders_slice(
-    info: &OperatorSliceInfo,
+    fac1: f64,
     subgrids: &ArrayView1<SubgridEnum>,
     orders: &[Order],
     order_mask: &[bool],
@@ -353,7 +353,7 @@ fn ndarray_from_subgrid_orders_slice(
         for ((ifac1, ix1, ix2), value) in subgrid.indexed_iter() {
             let Mu2 { ren, fac } = subgrid.mu2_grid()[ifac1];
 
-            if !approx_eq!(f64, xif * xif * fac, info.fac1, ulps = EVOLUTION_TOL_ULPS) {
+            if !approx_eq!(f64, xif * xif * fac, fac1, ulps = EVOLUTION_TOL_ULPS) {
                 continue;
             }
 
@@ -413,7 +413,7 @@ pub(crate) fn evolve_slice_with_one(
 
         for (subgrids_o, channel1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.channels()) {
             let (mut x1, array) = ndarray_from_subgrid_orders_slice(
-                info,
+                info.fac1,
                 &subgrids_o,
                 grid.orders(),
                 order_mask,
@@ -556,7 +556,7 @@ pub(crate) fn evolve_slice_with_two(
 
         for (subgrids_o, channel1) in subgrids_oc.axis_iter(Axis(1)).zip(grid.channels()) {
             let (x1, array) = ndarray_from_subgrid_orders_slice(
-                info,
+                info.fac1,
                 &subgrids_o,
                 grid.orders(),
                 order_mask,
@@ -633,3 +633,151 @@ pub(crate) fn evolve_slice_with_two(
             .collect(),
     ))
 }
+
+pub(crate) fn evolve_slice_with_two2(
+    grid: &Grid,
+    operators: &[ArrayView4<f64>],
+    infos: &[OperatorSliceInfo],
+    order_mask: &[bool],
+    xi: (f64, f64),
+    alphas_table: &AlphasTable,
+) -> Result<(Array3<SubgridEnum>, Vec<Channel>), GridError> {
+    let gluon_has_pid_zero = gluon_has_pid_zero(grid);
+
+    // TODO: implement matching of different scales for different EKOs
+    let mut fac1_scales: Vec<_> = infos.iter().map(|info| info.fac1).collect();
+    fac1_scales.sort_by(f64::total_cmp);
+    assert!(fac1_scales.windows(2).all(|scales| approx_eq!(
+        f64,
+        scales[0],
+        scales[1],
+        ulps = EVOLUTION_TOL_ULPS
+    )));
+    let fac1 = fac1_scales[0];
+
+    // TODO: generalize by iterating up to `n`
+    let (pid_indices, pids01): (Vec<_>, Vec<_>) = izip!(0..2, operators, infos)
+        .map(|(d, operator, info)| {
+            pid_slices(operator, info, gluon_has_pid_zero, &|pid1| {
+                grid.channels()
+                    .iter()
+                    .flat_map(Channel::entry)
+                    .any(|tuple| match d {
+                        // TODO: `Channel::entry` should return a tuple of a `Vec` and an `f64`
+                        0 => tuple.0 == pid1,
+                        1 => tuple.1 == pid1,
+                        _ => unreachable!(),
+                    })
+            })
+        })
+        .collect::<Result<Vec<_>, _>>()?
+        .into_iter()
+        .unzip();
+
+    let mut channels0: Vec<_> = pids01
+        .iter()
+        .map(|pids| pids.iter().map(|&(pid0, _)| pid0))
+        .multi_cartesian_product()
+        .collect();
+    channels0.sort_unstable();
+    channels0.dedup();
+    let channels0 = channels0;
+
+    let mut sub_fk_tables = Vec::with_capacity(grid.bin_info().bins() * channels0.len());
+
+    // TODO: generalize to `n`
+    let mut last_x1 = vec![Vec::new(); 2];
+    let mut eko_slices = vec![Vec::new(); 2];
+
+    for subgrids_oc in grid.subgrids().axis_iter(Axis(1)) {
+        assert_eq!(infos[0].x0.len(), infos[1].x0.len());
+
+        let mut tables =
+            vec![Array2::zeros((infos[0].x0.len(), infos[1].x0.len())); channels0.len()];
+
+        for (subgrids_o, channel1) in subgrids_oc.axis_iter(Axis(1)).zip(grid.channels()) {
+            let (x1, array) = ndarray_from_subgrid_orders_slice(
+                fac1,
+                &subgrids_o,
+                grid.orders(),
+                order_mask,
+                xi,
+                alphas_table,
+            )?;
+
+            for (last_x1, x1, pid_indices, slices, operator, info) in izip!(
+                &mut last_x1,
+                x1,
+                &pid_indices,
+                &mut eko_slices,
+                operators,
+                infos
+            ) {
+                if (last_x1.len() != x1.len())
+                    || last_x1
+                        .iter()
+                        .zip(x1.iter())
+                        .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS))
+                {
+                    *slices = operator_slices(operator, info, pid_indices, &x1)?;
+                    *last_x1 = x1;
+                }
+            }
+
+            let mut tmp = Array2::zeros((last_x1[0].len(), infos[1].x0.len()));
+
+            for (pids1, factor) in channel1
+                .entry()
+                .iter()
+                .map(|&(pida1, pidb1, factor)| ([pida1, pidb1], factor))
+            {
+                for (fk_table, ops) in
+                    channels0
+                        .iter()
+                        .zip(tables.iter_mut())
+                        .filter_map(|(pids0, fk_table)| {
+                            izip!(pids0, &pids1, &pids01, &eko_slices)
+                                .map(|(&pid0, &pid1, pids, slices)| {
+                                    pids.iter().zip(slices).find_map(|(&(p0, p1), op)| {
+                                        ((p0 == pid0) && (p1 == pid1)).then_some(op)
+                                    })
+                                })
+                                // TODO: avoid using `collect`
+                                .collect::<Option<Vec<_>>>()
+                                .map(|ops| (fk_table, ops))
+                        })
+                {
+                    // tmp = array * ops[1]^T
+                    linalg::general_mat_mul(1.0, &array, &ops[1].t(), 0.0, &mut tmp);
+                    // fk_table += factor * ops[0] * tmp
+                    linalg::general_mat_mul(factor, ops[0], &tmp, 1.0, fk_table);
+                }
+            }
+        }
+
+        sub_fk_tables.extend(tables.into_iter().map(|table| {
+            ImportOnlySubgridV2::new(
+                SparseArray3::from_ndarray(table.insert_axis(Axis(0)).view(), 0, 1),
+                vec![Mu2 {
+                    // TODO: FK tables don't depend on the renormalization scale
+                    //ren: -1.0,
+                    ren: infos[0].fac0,
+                    fac: infos[0].fac0,
+                }],
+                infos[0].x0.clone(),
+                infos[1].x0.clone(),
+            )
+            .into()
+        }));
+    }
+
+    Ok((
+        Array1::from_iter(sub_fk_tables)
+            .into_shape((1, grid.bin_info().bins(), channels0.len()))
+            .unwrap(),
+        channels0
+            .iter()
+            .map(|c| channel![c[0], c[1], 1.0])
+            .collect(),
+    ))
+}

pineappl/src/grid.rs

@@ -1489,6 +1489,145 @@ impl Grid {
         Ok(FkTable::try_from(grid).unwrap_or_else(|_| unreachable!()))
     }
 
+    /// Converts this `Grid` into an [`FkTable`] using `slices` that must iterate over a [`Result`]
+    /// of tuples of an [`OperatorSliceInfo`] and the corresponding sliced operator. The parameter
+    /// `order_mask` can be used to include or exclude orders from this operation, and must
+    /// correspond to the ordering given by [`Grid::orders`]. Orders that are not given are
+    /// enabled, and in particular if `order_mask` is empty all orders are activated.
+    ///
+    /// # Errors
+    ///
+    /// Returns a [`GridError::EvolutionFailure`] if either the `operator` or its `info` is
+    /// incompatible with this `Grid`. Returns a [`GridError::Other`] if the iterator from `slices`
+    /// return an error.
+    pub fn evolve_with_slice_iter2<'a, E: Into<anyhow::Error>>(
+        &self,
+        slices_a: impl IntoIterator<Item = Result<(OperatorSliceInfo, CowArray<'a, f64, Ix4>), E>>,
+        slices_b: impl IntoIterator<Item = Result<(OperatorSliceInfo, CowArray<'a, f64, Ix4>), E>>,
+        order_mask: &[bool],
+        xi: (f64, f64),
+        alphas_table: &AlphasTable,
+    ) -> Result<FkTable, GridError> {
+        use super::evolution::EVOLVE_INFO_TOL_ULPS;
+        use itertools::izip;
+
+        let mut lhs: Option<Self> = None;
+        let mut fac1 = Vec::new();
+
+        // TODO: simplify the ugly repetition below by offloading some ops into fn
+        for (result_a, result_b) in izip!(slices_a, slices_b) {
+            // Operate on `slices_a`
+            let (info_a, operator_a) = result_a.map_err(|err| GridError::Other(err.into()))?;
+
+            let op_info_dim_a = (
+                info_a.pids1.len(),
+                info_a.x1.len(),
+                info_a.pids0.len(),
+                info_a.x0.len(),
+            );
+
+            if operator_a.dim() != op_info_dim_a {
+                return Err(GridError::EvolutionFailure(format!(
+                    "operator information {:?} does not match the operator's dimensions: {:?}",
+                    op_info_dim_a,
+                    operator_a.dim(),
+                )));
+            }
+
+            // Operate on `slices_b`
+            let (info_b, operator_b) = result_b.map_err(|err| GridError::Other(err.into()))?;
+
+            let op_info_dim_b = (
+                info_b.pids1.len(),
+                info_b.x1.len(),
+                info_b.pids0.len(),
+                info_b.x0.len(),
+            );
+
+            if operator_b.dim() != op_info_dim_b {
+                return Err(GridError::EvolutionFailure(format!(
+                    "operator information {:?} does not match the operator's dimensions: {:?}",
+                    op_info_dim_b,
+                    operator_b.dim(),
+                )));
+            }
+
+            let views = [operator_a.view(), operator_b.view()];
+            let infos = [info_a, info_b];
+
+            let (subgrids, channels) = if self.convolutions()[0] != Convolution::None
+                && self.convolutions()[1] != Convolution::None
+            {
+                evolution::evolve_slice_with_two2(
+                    self,
+                    &views,
+                    &infos,
+                    order_mask,
+                    xi,
+                    alphas_table,
+                )
+            } else {
+                evolution::evolve_slice_with_one(
+                    self,
+                    &views[0],
+                    &infos[1],
+                    order_mask,
+                    xi,
+                    alphas_table,
+                )
+            }?;
+
+            let mut rhs = Self {
+                subgrids,
+                channels,
+                bin_limits: self.bin_limits.clone(),
+                orders: vec![Order::new(0, 0, 0, 0)],
+                subgrid_params: SubgridParams::default(),
+                more_members: self.more_members.clone(),
+            };
+
+            assert_eq!(infos[0].pid_basis, infos[1].pid_basis);
+
+            // TODO: use a new constructor to set this information
+            rhs.set_pid_basis(infos[0].pid_basis);
+
+            if let Some(lhs) = &mut lhs {
+                lhs.merge(rhs)?;
+            } else {
+                lhs = Some(rhs);
+            }
+
+            // NOTE: The following should be shared by the 2 EKOs(?)
+            fac1.push(infos[0].fac1);
+        }
+
+        // UNWRAP: if we can't compare two numbers there's a bug
+        fac1.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!()));
+
+        // make sure we've evolved all slices
+        if let Some(muf2) = self
+            .evolve_info(order_mask)
+            .fac1
+            .into_iter()
+            .map(|mu2| xi.1 * xi.1 * mu2)
+            .find(|&grid_mu2| {
+                !fac1
+                    .iter()
+                    .any(|&eko_mu2| approx_eq!(f64, grid_mu2, eko_mu2, ulps = EVOLVE_INFO_TOL_ULPS))
+            })
+        {
+            return Err(GridError::EvolutionFailure(format!(
+                "no operator for muf2 = {muf2} found in {fac1:?}"
+            )));
+        }
+
+        // TODO: convert this unwrap into error
+        let grid = lhs.unwrap();
+
+        // UNWRAP: merging evolved slices should be a proper FkTable again
+        Ok(FkTable::try_from(grid).unwrap_or_else(|_| unreachable!()))
+    }
+
     /// Deletes bins with the corresponding `bin_indices`. Repeated indices and indices larger or
     /// equal the bin length are ignored.
     pub fn delete_bins(&mut self, bin_indices: &[usize]) {