From 53acd7f9b930561741e9740480fca0d820068c7c Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 10 Oct 2023 13:34:01 +0200 Subject: [PATCH 01/40] Add first candidate for `Grid::evolve_slice` --- pineappl/src/evolution.rs | 331 +++++++++++++++++++++++++++++++++++++- pineappl/src/grid.rs | 57 ++++++- 2 files changed, 385 insertions(+), 3 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index c93b5603..dfa12c76 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -8,8 +8,9 @@ use super::sparse_array3::SparseArray3; use super::subgrid::{Mu2, Subgrid, SubgridEnum}; use float_cmp::approx_eq; use itertools::Itertools; -use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView5, Axis}; +use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView4, ArrayView5, Axis, Slice}; use std::iter; +use std::slice; /// Number of ULPS used to de-duplicate grid values in [`Grid::evolve_info`]. pub(crate) const EVOLVE_INFO_TOL_ULPS: i64 = 64; @@ -88,6 +89,62 @@ pub struct OperatorInfo { pub lumi_id_types: String, } +/// Information about the evolution kernel operator slice (EKO) passed to [`Grid::evolve_slice`] as +/// `operator`, which is used to convert a [`Grid`] into an [`FkTable`]. The dimensions of the EKO +/// must correspond to the values given in [`fac1`], [`pids0`], [`x0`], [`pids1`] and [`x1`], +/// exactly in this order. Members with a `1` are defined at the squared factorization scale given +/// as [`fac1`] (often called process scale) and are found in the [`Grid`] that +/// [`Grid::evolve_slice`] is called with. Members with a `0` are defined at the squared +/// factorization scale [`fac0`] (often called fitting scale or starting scale) and are found in +/// the [`FkTable`] resulting from [`Grid::evolve`]. +/// +/// The EKO slice may convert a `Grid` from a basis given by the particle identifiers [`pids1`] to +/// a possibly different basis given by [`pids0`]. This basis must also be identified using +/// [`lumi_id_types`], which tells [`FkTable::convolute`] how to perform a convolution. The members +/// [`ren1`] and [`alphas`] must be the strong couplings given at the respective renormalization +/// scales. Finally, [`xir`] and [`xif`] can be used to vary the renormalization and factorization +/// scales, respectively, around their central values. +/// +/// [`FkTable::convolute`]: super::fk_table::FkTable::convolute +/// [`FkTable`]: super::fk_table::FkTable +/// [`alphas`]: Self::alphas +/// [`fac0`]: Self::fac0 +/// [`fac1`]: Self::fac1 +/// [`lumi_id_types`]: Self::lumi_id_types +/// [`pids0`]: Self::pids0 +/// [`pids1`]: Self::pids1 +/// [`ren1`]: Self::ren1 +/// [`x0`]: Self::x0 +/// [`x1`]: Self::x1 +/// [`xif`]: Self::xif +/// [`xir`]: Self::xir +pub struct OperatorSliceInfo { + /// Squared factorization scale of the `FkTable`. + pub fac0: f64, + /// Particle identifiers of the `FkTable`. + pub pids0: Vec, + /// `x`-grid coordinates of the `FkTable` + pub x0: Vec, + /// Squared factorization scale of the slice of `Grid` that should be evolved. + pub fac1: f64, + /// Particle identifiers of the `Grid`. If the `Grid` contains more particle identifiers than + /// given here, the contributions of them are silently ignored. + pub pids1: Vec, + /// `x`-grid coordinates of the `Grid`. + pub x1: Vec, + + /// Renormalization scales of the `Grid`. + pub ren1: Vec, + /// Strong couplings corresponding to the order given in [`ren1`](Self::ren1). + pub alphas: Vec, + /// Multiplicative factor for the central renormalization scale. + pub xir: f64, + /// Multiplicative factor for the central factorization scale. + pub xif: f64, + /// Identifier of the particle basis for the `FkTable`. + pub lumi_id_types: String, +} + fn gluon_has_pid_zero(grid: &Grid) -> bool { grid.key_values() .and_then(|key_values| key_values.get("lumi_id_types")) @@ -368,6 +425,127 @@ pub(crate) fn ndarray_from_subgrid_orders( Ok((fac1, x1_a, x1_b, array)) } +type X1aX1bOp3Tuple = (Vec, Vec, Array3); + +fn ndarray_from_subgrid_orders_slice( + info: &OperatorSliceInfo, + subgrids: &ArrayView1, + orders: &[Order], + order_mask: &[bool], +) -> Result { + // TODO: skip empty subgrids + + let fac1 = info.xif * info.xif * info.fac1; + let mut x1_a: Vec<_> = subgrids + .iter() + .enumerate() + .filter(|(index, _)| order_mask.get(*index).copied().unwrap_or(true)) + .flat_map(|(_, subgrid)| subgrid.x1_grid().into_owned()) + .collect(); + let mut x1_b: Vec<_> = subgrids + .iter() + .enumerate() + .filter(|(index, _)| order_mask.get(*index).copied().unwrap_or(true)) + .flat_map(|(_, subgrid)| subgrid.x2_grid().into_owned()) + .collect(); + + x1_a.sort_by(f64::total_cmp); + x1_a.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLUTION_TOL_ULPS)); + x1_b.sort_by(f64::total_cmp); + x1_b.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLUTION_TOL_ULPS)); + + let mut array = Array3::::zeros((1, x1_a.len(), x1_b.len())); + + // add subgrids for different orders, but the same bin and lumi, using the right + // couplings + for (subgrid, order) in subgrids + .iter() + .zip(orders.iter()) + .zip(order_mask.iter().chain(iter::repeat(&true))) + .filter_map(|((subgrid, order), &enabled)| { + (enabled && !subgrid.is_empty()).then_some((subgrid, order)) + }) + { + let mut logs = 1.0; + + if order.logxir > 0 { + if approx_eq!(f64, info.xir, 1.0, ulps = 4) { + continue; + } + + logs *= (info.xir * info.xir).ln(); + } + + if order.logxif > 0 { + if approx_eq!(f64, info.xif, 1.0, ulps = 4) { + continue; + } + + logs *= (info.xif * info.xif).ln(); + } + + // TODO: use `try_collect` once stabilized + let xa_indices: Vec<_> = subgrid + .x1_grid() + .iter() + .map(|&xa| { + x1_a.iter() + .position(|&x1a| approx_eq!(f64, x1a, xa, ulps = EVOLUTION_TOL_ULPS)) + .ok_or_else(|| { + GridError::EvolutionFailure(format!("no operator for x1 = {xa} found")) + }) + }) + .collect::>()?; + let xb_indices: Vec<_> = subgrid + .x2_grid() + .iter() + .map(|&xb| { + x1_b.iter() + .position(|&x1b| approx_eq!(f64, x1b, xb, ulps = EVOLUTION_TOL_ULPS)) + .ok_or_else(|| { + GridError::EvolutionFailure(format!("no operator for x1 = {xb} found")) + }) + }) + .collect::>()?; + + for ((ifac1, ix1, ix2), value) in subgrid.indexed_iter() { + let Mu2 { ren, fac } = subgrid.mu2_grid()[ifac1]; + + if !approx_eq!( + f64, + info.xif * info.xif * fac, + fac1, + ulps = EVOLUTION_TOL_ULPS + ) { + continue; + } + + let mur2 = info.xir * info.xir * ren; + + let als = if order.alphas == 0 { + 1.0 + } else if let Some(alphas) = + info.ren1 + .iter() + .zip(info.alphas.iter()) + .find_map(|(&ren1, &alphas)| { + approx_eq!(f64, ren1, mur2, ulps = EVOLUTION_TOL_ULPS).then(|| alphas) + }) + { + alphas.powi(order.alphas.try_into().unwrap()) + } else { + return Err(GridError::EvolutionFailure(format!( + "no alphas for mur2 = {mur2} found" + ))); + }; + + array[[0, xa_indices[ix1], xb_indices[ix2]]] += als * logs * value; + } + } + + Ok((x1_a, x1_b, array)) +} + pub(crate) fn evolve_with_one( grid: &Grid, operator: &ArrayView5, @@ -634,3 +812,154 @@ pub(crate) fn evolve_with_two( lumi0.iter().map(|&(a, b)| lumi_entry![a, b, 1.0]).collect(), )) } + +pub(crate) fn evolve_slice_with_two( + grid: &Grid, + operator: &ArrayView4, + info: &OperatorSliceInfo, + order_mask: &[bool], +) -> Result<(Array3, Vec), GridError> { + let op5 = operator.insert_axis(Axis(0)); + let inf = OperatorInfo { + fac0: info.fac0, + pids0: info.pids0.clone(), + x0: info.x0.clone(), + fac1: vec![info.fac1], + pids1: info.pids1.clone(), + x1: info.x1.clone(), + ren1: info.ren1.clone(), + alphas: info.alphas.clone(), + xir: info.xir, + xif: info.xif, + lumi_id_types: info.lumi_id_types.clone(), + }; + let gluon_has_pid_zero = gluon_has_pid_zero(grid); + + let (pid_indices_a, pids_a) = pids(&op5, &inf, gluon_has_pid_zero, &|pid1| { + grid.lumi() + .iter() + .flat_map(LumiEntry::entry) + .any(|&(a, _, _)| a == pid1) + })?; + let (pid_indices_b, pids_b) = pids(&op5, &inf, gluon_has_pid_zero, &|pid1| { + grid.lumi() + .iter() + .flat_map(LumiEntry::entry) + .any(|&(_, b, _)| b == pid1) + })?; + + let lumi0 = lumi0_with_two(&pids_a, &pids_b); + let mut sub_fk_tables = Vec::with_capacity(grid.bin_info().bins() * lumi0.len()); + + let mut last_x1a = Vec::new(); + let mut last_x1b = Vec::new(); + let mut operators_a = Vec::new(); + let mut operators_b = Vec::new(); + + for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { + let mut tables = vec![Array2::zeros((info.x0.len(), info.x0.len())); lumi0.len()]; + + for (lumi1, subgrids_o) in subgrids_ol.axis_iter(Axis(1)).enumerate() { + let (x1_a, x1_b, array) = + ndarray_from_subgrid_orders_slice(info, &subgrids_o, grid.orders(), order_mask)?; + + if (last_x1a.len() != x1_a.len()) + || last_x1a + .iter() + .zip(x1_a.iter()) + .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) + { + operators_a = operators( + &op5, + &inf, + slice::from_ref(&info.fac1), + &pid_indices_a, + &x1_a, + )? + .into_iter() + .map(|mut op| { + op.slice_axis_inplace(Axis(0), Slice::from(0..1)); + op + }) + .collect(); + last_x1a = x1_a; + } + + if (last_x1b.len() != x1_b.len()) + || last_x1b + .iter() + .zip(x1_b.iter()) + .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) + { + operators_b = operators( + &op5, + &inf, + slice::from_ref(&info.fac1), + &pid_indices_b, + &x1_b, + )? + .into_iter() + .map(|mut op| { + op.slice_axis_inplace(Axis(0), Slice::from(0..1)); + op + }) + .collect(); + last_x1b = x1_b; + } + + // TODO: get rid of array-index access + for &(pida1, pidb1, factor) in grid.lumi()[lumi1].entry() { + for (fk_table, opa, opb) in + lumi0 + .iter() + .zip(tables.iter_mut()) + .filter_map(|(&(pida0, pidb0), fk_table)| { + pids_a + .iter() + .zip(operators_a.iter()) + .cartesian_product(pids_b.iter().zip(operators_b.iter())) + .find_map(|((&(pa0, pa1), opa), (&(pb0, pb1), opb))| { + (pa0 == pida0 && pa1 == pida1 && pb0 == pidb0 && pb1 == pidb1) + .then_some((opa, opb)) + }) + .map(|(opa, opb)| (fk_table, opa, opb)) + }) + { + let mut result = Array2::zeros((info.x0.len(), info.x0.len())); + + for imu2 in 0..array.dim().0 { + let opa = opa.index_axis(Axis(0), imu2); + let opb = opb.index_axis(Axis(0), imu2); + let arr = array.index_axis(Axis(0), imu2); + + result += &opa.dot(&arr.dot(&opb.t())); + } + + fk_table.scaled_add(factor, &result); + } + } + } + + 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: info.fac0, + fac: info.fac0, + }], + info.x0.clone(), + info.x0.clone(), + ) + .into() + })); + } + + Ok(( + Array1::from_iter(sub_fk_tables) + .into_shape((1, grid.bin_info().bins(), lumi0.len())) + .unwrap(), + lumi0.iter().map(|&(a, b)| lumi_entry![a, b, 1.0]).collect(), + )) +} diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 3a785fb8..023761cb 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -2,7 +2,7 @@ use super::bin::{BinInfo, BinLimits, BinRemapper}; use super::empty_subgrid::EmptySubgridV1; -use super::evolution::{self, EvolveInfo, OperatorInfo}; +use super::evolution::{self, EvolveInfo, OperatorInfo, OperatorSliceInfo}; use super::fk_table::FkTable; use super::import_only_subgrid::ImportOnlySubgridV2; use super::lagrange_subgrid::{LagrangeSparseSubgridV1, LagrangeSubgridV1, LagrangeSubgridV2}; @@ -18,7 +18,7 @@ use git_version::git_version; use indicatif::{ProgressBar, ProgressStyle}; use itertools::Itertools; use lz4_flex::frame::{FrameDecoder, FrameEncoder}; -use ndarray::{s, Array3, Array5, ArrayView5, Axis, Dimension}; +use ndarray::{s, Array3, Array5, ArrayView4, ArrayView5, Axis, Dimension}; use serde::{Deserialize, Serialize}; use std::borrow::Cow; use std::cmp::Ordering; @@ -1945,6 +1945,59 @@ impl Grid { Ok(FkTable::try_from(grid).unwrap_or_else(|_| unreachable!())) } + /// Converts this `Grid` into an [`FkTable`] using an evolution kernel operator (EKO) slice + /// given as `operator`. The dimensions and properties of this operator must be described using + /// `info`. 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`. + pub fn evolve_slice( + &self, + operator: ArrayView4, + info: &OperatorSliceInfo, + order_mask: &[bool], + ) -> Result { + let op_info_dim = ( + info.pids1.len(), + info.x1.len(), + info.pids0.len(), + info.x0.len(), + ); + + if operator.dim() != op_info_dim { + return Err(GridError::EvolutionFailure(format!( + "operator information {:?} does not match the operator's dimensions: {:?}", + op_info_dim, + operator.dim(), + ))); + } + + let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { + evolution::evolve_slice_with_two(self, &operator, info, order_mask) + } else { + todo!() + //evolution::evolve_slice_with_one(self, &operator, info, order_mask) + }?; + + let mut grid = Self { + subgrids, + lumi, + bin_limits: self.bin_limits.clone(), + orders: vec![Order::new(0, 0, 0, 0)], + subgrid_params: SubgridParams::default(), + more_members: self.more_members.clone(), + }; + + // write additional metadata + grid.set_key_value("lumi_id_types", &info.lumi_id_types); + + 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]) { From b7fa693c957455c9985724365854cbe2c9a6ce83 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 10 Oct 2023 15:17:35 +0200 Subject: [PATCH 02/40] Add support for EKO format introduced in v0.13 --- CHANGELOG.md | 4 ++ pineappl_cli/src/evolve.rs | 86 ++++++++++++++++++++++++++++++++++++ pineappl_cli/tests/evolve.rs | 47 ++++++++++++++++++++ 3 files changed, 137 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index dc9c5a15..6b222190 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- added support for new EKO format introduced with EKO v0.13 + ## [0.6.2] - 09/10/2023 ### Added diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 5c79e869..761224c2 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -22,6 +22,7 @@ mod eko { use pineappl::evolution::OperatorInfo; use pineappl::pids; use serde::Deserialize; + use std::collections::HashMap; use std::fs::File; use std::io::BufReader; use std::path::Path; @@ -58,11 +59,38 @@ mod eko { rotations: Rotations, } + #[derive(Deserialize)] + struct OperatorV1 { + mu0: f64, + } + + #[derive(Deserialize)] + struct OperatorInfoV1 { + scale: f64, + } + + const BASES_V1_DEFAULT_PIDS: [i32; 14] = [22, -6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6]; + + #[derive(Deserialize)] + struct BasesV1 { + inputgrid: Option>, + inputpids: Option>>, + targetgrid: Option>, + targetpids: Option>, + xgrid: Vec, + } + + #[derive(Deserialize)] + struct MetadataV2 { + bases: BasesV1, + } + #[derive(Deserialize)] #[serde(untagged)] enum Metadata { V0(MetadataV0), V1(MetadataV1), + V2(MetadataV2), } pub fn read(eko: &Path) -> Result<(OperatorInfo, Array5)> { @@ -72,6 +100,9 @@ mod eko { let mut operator = None; let mut operators = Vec::new(); let mut fac1 = Vec::new(); + let mut fac0 = None; + let mut fac1_order = Vec::new(); + let mut fac1_translator = HashMap::new(); let base64 = GeneralPurpose::new(&URL_SAFE, PAD); @@ -83,12 +114,18 @@ mod eko { // TODO: get rid of the unwrap match file_name.to_str().unwrap() { "metadata.yaml" => metadata = Some(serde_yaml::from_reader(file)?), + "operator.yaml" => { + let operator: OperatorV1 = serde_yaml::from_reader(file)?; + fac0 = Some(operator.mu0 * operator.mu0); + } "operators.npy.lz4" => { operator = Some(Array5::::read_npy(FrameDecoder::new( BufReader::new(file), ))?); } x if x.ends_with(".npz.lz4") => { + fac1_order.push(x.strip_suffix(".npz.lz4").unwrap().to_string()); + let name = x.strip_suffix(".npz.lz4").unwrap(); let bytes = base64.decode(name.as_bytes())?; let array: [u8; 8] = bytes.as_slice().try_into().unwrap(); @@ -103,6 +140,11 @@ mod eko { fac1.push(muf2); operators.push(operator); } + x if x.ends_with("=.yaml") => { + let name = x.strip_suffix(".yaml").unwrap().to_string(); + let op_info: OperatorInfoV1 = serde_yaml::from_reader(file)?; + fac1_translator.insert(name, op_info.scale); + } _ => {} } } @@ -113,6 +155,12 @@ mod eko { operator = Some(ndarray::stack(Axis(0), &ops).unwrap()); } + if !fac1_translator.is_empty() { + for (fac1, fac1_name) in fac1.iter_mut().zip(&fac1_order) { + *fac1 = fac1_translator[fac1_name]; + } + } + let mut info = match metadata { Some(Metadata::V0(metadata)) => OperatorInfo { fac1: metadata.q2_grid.clone(), @@ -169,6 +217,44 @@ mod eko { xif: 1.0, lumi_id_types: String::new(), }, + Some(Metadata::V2(metadata)) => OperatorInfo { + fac1, + pids0: metadata.bases.inputpids.map_or_else( + || BASES_V1_DEFAULT_PIDS.to_vec(), + |basis| { + basis + .into_iter() + .map(|factors| { + let tuples: Vec<_> = BASES_V1_DEFAULT_PIDS + .iter() + .copied() + .zip(factors.into_iter()) + .collect(); + + pids::pdg_mc_ids_to_evol(&tuples).unwrap() + }) + .collect() + }, + ), + x0: metadata + .bases + .inputgrid + .unwrap_or_else(|| metadata.bases.xgrid.clone()), + pids1: metadata + .bases + .targetpids + .unwrap_or_else(|| BASES_V1_DEFAULT_PIDS.to_vec()), + x1: metadata + .bases + .targetgrid + .unwrap_or_else(|| metadata.bases.xgrid.clone()), + fac0: fac0.unwrap(), + ren1: vec![], + alphas: vec![], + xir: 1.0, + xif: 1.0, + lumi_id_types: String::new(), + }, None => bail!("no `metadata.yaml` file found"), }; diff --git a/pineappl_cli/tests/evolve.rs b/pineappl_cli/tests/evolve.rs index a220cecc..61eb8ccf 100644 --- a/pineappl_cli/tests/evolve.rs +++ b/pineappl_cli/tests/evolve.rs @@ -67,6 +67,18 @@ const LHCB_WP_7TEV_STR: &str = "b Grid FkTable rel. diff 7 2.9272102e1 2.9268451e1 -1.2474967e-4 "; +const LHCB_WP_7TEV_V2_STR: &str = "b Grid FkTable rel. diff +-+-----------+-----------+------------- +0 7.8752127e2 7.8731064e2 -2.6745204e-4 +1 7.1872113e2 7.1853123e2 -2.6421837e-4 +2 6.2322357e2 6.2306010e2 -2.6230496e-4 +3 5.0216763e2 5.0203737e2 -2.5938800e-4 +4 3.7314506e2 3.7305090e2 -2.5233796e-4 +5 2.5302044e2 2.5295968e2 -2.4013733e-4 +6 1.1971046e2 1.1968525e2 -2.1055575e-4 +7 2.9272102e1 2.9268443e1 -1.2499435e-4 +"; + const LHCB_WP_7TEV_OPTIMIZED_STR: &str = "b etal dsig/detal [] [pb] -+----+----+----------- @@ -191,6 +203,41 @@ fn lhcb_wp_7tev() { .stdout(LHCB_WP_7TEV_OPTIMIZED_STR); } +#[test] +fn lhcb_wp_7tev_v2() { + let input = NamedTempFile::new("optimized.pineappl.lz4").unwrap(); + + // we first need to optimize the grid, to strip empty x-grid values not contained in the EKO + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "write", + "--optimize", + "../test-data/LHCB_WP_7TEV.pineappl.lz4", + input.path().to_str().unwrap(), + ]) + .assert() + .success() + .stdout(""); + + let output = NamedTempFile::new("fktable2a.lz4").unwrap(); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--silence-lhapdf", + "evolve", + input.path().to_str().unwrap(), + "../test-data/LHCB_WP_7TEV_v2.tar", + output.path().to_str().unwrap(), + "NNPDF40_nlo_as_01180", + "--orders=a2,as1a2", + ]) + .assert() + .success() + .stdout(LHCB_WP_7TEV_V2_STR); +} + #[test] fn e906nlo_bin_00() { let output = NamedTempFile::new("fktable2.lz4").unwrap(); From 11f4792a2c052192385130ff7d9a1e6bc3515ee3 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 11 Oct 2023 11:37:07 +0200 Subject: [PATCH 03/40] Fix dimensions of operator slices --- pineappl/src/evolution.rs | 31 +++---------------------------- 1 file changed, 3 insertions(+), 28 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index dfa12c76..37b2c03b 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -8,9 +8,8 @@ use super::sparse_array3::SparseArray3; use super::subgrid::{Mu2, Subgrid, SubgridEnum}; use float_cmp::approx_eq; use itertools::Itertools; -use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView4, ArrayView5, Axis, Slice}; +use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView4, ArrayView5, Axis}; use std::iter; -use std::slice; /// Number of ULPS used to de-duplicate grid values in [`Grid::evolve_info`]. pub(crate) const EVOLVE_INFO_TOL_ULPS: i64 = 64; @@ -869,19 +868,7 @@ pub(crate) fn evolve_slice_with_two( .zip(x1_a.iter()) .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) { - operators_a = operators( - &op5, - &inf, - slice::from_ref(&info.fac1), - &pid_indices_a, - &x1_a, - )? - .into_iter() - .map(|mut op| { - op.slice_axis_inplace(Axis(0), Slice::from(0..1)); - op - }) - .collect(); + operators_a = operators(&op5, &inf, &inf.fac1, &pid_indices_a, &x1_a)?; last_x1a = x1_a; } @@ -891,19 +878,7 @@ pub(crate) fn evolve_slice_with_two( .zip(x1_b.iter()) .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) { - operators_b = operators( - &op5, - &inf, - slice::from_ref(&info.fac1), - &pid_indices_b, - &x1_b, - )? - .into_iter() - .map(|mut op| { - op.slice_axis_inplace(Axis(0), Slice::from(0..1)); - op - }) - .collect(); + operators_b = operators(&op5, &inf, &inf.fac1, &pid_indices_b, &x1_b)?; last_x1b = x1_b; } From b9170bc09b82fee27c339f626af4f01e2d08f372 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 11 Oct 2023 11:39:21 +0200 Subject: [PATCH 04/40] Test more digits of evolution --- pineappl_cli/tests/evolve.rs | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/pineappl_cli/tests/evolve.rs b/pineappl_cli/tests/evolve.rs index 61eb8ccf..e69d2769 100644 --- a/pineappl_cli/tests/evolve.rs +++ b/pineappl_cli/tests/evolve.rs @@ -67,16 +67,16 @@ const LHCB_WP_7TEV_STR: &str = "b Grid FkTable rel. diff 7 2.9272102e1 2.9268451e1 -1.2474967e-4 "; -const LHCB_WP_7TEV_V2_STR: &str = "b Grid FkTable rel. diff --+-----------+-----------+------------- -0 7.8752127e2 7.8731064e2 -2.6745204e-4 -1 7.1872113e2 7.1853123e2 -2.6421837e-4 -2 6.2322357e2 6.2306010e2 -2.6230496e-4 -3 5.0216763e2 5.0203737e2 -2.5938800e-4 -4 3.7314506e2 3.7305090e2 -2.5233796e-4 -5 2.5302044e2 2.5295968e2 -2.4013733e-4 -6 1.1971046e2 1.1968525e2 -2.1055575e-4 -7 2.9272102e1 2.9268443e1 -1.2499435e-4 +const LHCB_WP_7TEV_V2_STR: &str = "b Grid FkTable rel. diff +-+--------------------+--------------------+---------------------- +0 7.8752126798068639e2 7.8731064380928558e2 -2.6745204220435248e-4 +1 7.1872113080347663e2 7.1853123147848032e2 -2.6421836906898033e-4 +2 6.2322357391848550e2 6.2306009928459093e2 -2.6230495882362259e-4 +3 5.0216762988872915e2 5.0203737363369049e2 -2.5938799573266280e-4 +4 3.7314505699003126e2 3.7305089832847733e2 -2.5233795755852384e-4 +5 2.5302044227292129e2 2.5295968261889854e2 -2.4013733229188983e-4 +6 1.1971045984774410e2 1.1968525412249538e2 -2.1055574659711862e-4 +7 2.9272102213930090e1 2.9268443366651141e1 -1.2499434622803562e-4 "; const LHCB_WP_7TEV_OPTIMIZED_STR: &str = "b etal dsig/detal @@ -227,6 +227,8 @@ fn lhcb_wp_7tev_v2() { .args([ "--silence-lhapdf", "evolve", + "--digits-abs=16", + "--digits-rel=16", input.path().to_str().unwrap(), "../test-data/LHCB_WP_7TEV_v2.tar", output.path().to_str().unwrap(), From 169168f254f9eba3374f6757dd8252ce3aec3dab Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 11 Oct 2023 16:35:15 +0200 Subject: [PATCH 05/40] Add missing test file --- .github/workflows/rust.yml | 3 ++- maintainer/generate-coverage.sh | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index d2dd409e..76f661a6 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -23,7 +23,7 @@ jobs: uses: actions/cache@v3 with: path: test-data - key: test-data-v8 + key: test-data-v9 - name: Download test data if: steps.cache-test-data.outputs.cache-hit != 'true' run: | @@ -39,6 +39,7 @@ jobs: curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCBWZMU7TEV_PI_part1.appl' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.tar' + curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV_v2.tar' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_old.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.tar' diff --git a/maintainer/generate-coverage.sh b/maintainer/generate-coverage.sh index 8f0c3bdc..bf071037 100755 --- a/maintainer/generate-coverage.sh +++ b/maintainer/generate-coverage.sh @@ -17,6 +17,7 @@ wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_old.pineappl.lz4' wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.pineappl.lz4' wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.tar' +wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_v2.tar' wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/NJetEvents_0-0-2.tab.gz' wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/NUTEV_CC_NU_FE_SIGMARED.pineappl.lz4' wget --no-verbose --no-clobber -P test-data 'https://data.nnpdf.science/pineappl/test-data/NUTEV_CC_NU_FE_SIGMARED.tar' From d01f576622ae285618f6e0a5129018f9e30c002b Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 11 Oct 2023 16:49:09 +0200 Subject: [PATCH 06/40] Fix test file name --- .github/workflows/rust.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 76f661a6..a61db8ab 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -23,7 +23,7 @@ jobs: uses: actions/cache@v3 with: path: test-data - key: test-data-v9 + key: test-data-v10 - name: Download test data if: steps.cache-test-data.outputs.cache-hit != 'true' run: | @@ -39,10 +39,10 @@ jobs: curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCBWZMU7TEV_PI_part1.appl' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.tar' - curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV_v2.tar' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_old.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.tar' + curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_v2.tar' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/NJetEvents_0-0-2.tab.gz' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/NUTEV_CC_NU_FE_SIGMARED.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/NUTEV_CC_NU_FE_SIGMARED.tar' From c15de718a4fdc7381601481ba7579af703032912 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 12 Oct 2023 10:38:16 +0200 Subject: [PATCH 07/40] Add `EkoSlices` and `EkoSlicesIter` structs --- pineappl/src/evolution.rs | 1 + pineappl_cli/src/evolve.rs | 185 ++++++++++++++++++++++++++++++++++--- 2 files changed, 174 insertions(+), 12 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 37b2c03b..695681f0 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -117,6 +117,7 @@ pub struct OperatorInfo { /// [`x1`]: Self::x1 /// [`xif`]: Self::xif /// [`xir`]: Self::xir +#[derive(Clone)] pub struct OperatorSliceInfo { /// Squared factorization scale of the `FkTable`. pub fac0: f64, diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 761224c2..337e3f16 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -10,7 +10,7 @@ use std::process::ExitCode; #[cfg(feature = "evolve")] mod eko { - use anyhow::{bail, Result}; + use anyhow::{anyhow, bail, Result}; use base64::alphabet::URL_SAFE; use base64::engine::general_purpose::PAD; use base64::engine::GeneralPurpose; @@ -19,14 +19,15 @@ mod eko { use lz4_flex::frame::FrameDecoder; use ndarray::{Array4, Array5, Axis}; use ndarray_npy::{NpzReader, ReadNpyExt}; - use pineappl::evolution::OperatorInfo; + use pineappl::evolution::{OperatorInfo, OperatorSliceInfo}; use pineappl::pids; use serde::Deserialize; use std::collections::HashMap; + use std::ffi::{OsStr, OsString}; use std::fs::File; - use std::io::BufReader; + use std::io::{self, BufReader, Cursor}; use std::path::Path; - use tar::Archive; + use tar::{Archive, Entries}; #[derive(Deserialize)] struct MetadataV0 { @@ -59,6 +60,18 @@ mod eko { rotations: Rotations, } + #[derive(Deserialize)] + #[serde(untagged)] + enum Metadata { + V0(MetadataV0), + V1(MetadataV1), + V2(MetadataV2), + } + + // --- EKO FORMAT INTRODUCED WITH EKO v0.13 --- + + const BASES_V1_DEFAULT_PIDS: [i32; 14] = [22, -6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6]; + #[derive(Deserialize)] struct OperatorV1 { mu0: f64, @@ -69,8 +82,6 @@ mod eko { scale: f64, } - const BASES_V1_DEFAULT_PIDS: [i32; 14] = [22, -6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6]; - #[derive(Deserialize)] struct BasesV1 { inputgrid: Option>, @@ -85,14 +96,164 @@ mod eko { bases: BasesV1, } - #[derive(Deserialize)] - #[serde(untagged)] - enum Metadata { - V0(MetadataV0), - V1(MetadataV1), - V2(MetadataV2), + pub struct EkoSlices { + fac1: HashMap, + info: OperatorSliceInfo, + archive: Archive, } + impl EkoSlices { + pub fn new( + eko_path: &Path, + ren1: &[f64], + alphas: &[f64], + xir: f64, + xif: f64, + ) -> Result { + let mut fac1 = HashMap::new(); + let mut metadata: Option = None; + let mut operator: Option = None; + + for entry in Archive::new(File::open(eko_path)?).entries_with_seek()? { + let entry = entry?; + let path = entry.path()?; + + if path.starts_with("./operators") && (path.extension() == Some(OsStr::new("yaml"))) + { + // TODO: use let-else when available in MSRV + let file_stem = if let Some(file_stem) = path.file_stem() { + file_stem.to_os_string() + } else { + continue; + }; + + let op_info: OperatorInfoV1 = serde_yaml::from_reader(entry)?; + fac1.insert(file_stem, op_info.scale); + } else if path.as_os_str() == OsStr::new("./metadata.yaml") { + metadata = Some(serde_yaml::from_reader(entry)?); + } else if path.as_os_str() == OsStr::new("./operator.yaml") { + operator = Some(serde_yaml::from_reader(entry)?); + } + } + + let metadata = + metadata.ok_or_else(|| anyhow!("no file 'metadata.yaml' in EKO archive found"))?; + let operator = + operator.ok_or_else(|| anyhow!("no file 'operator.yaml' in EKO archive found"))?; + + let pids0 = metadata.bases.inputpids.map_or_else( + || BASES_V1_DEFAULT_PIDS.to_vec(), + |basis| { + basis + .into_iter() + .map(|factors| { + let tuples: Vec<_> = BASES_V1_DEFAULT_PIDS + .iter() + .copied() + .zip(factors.into_iter()) + .collect(); + + // UNWRAP: we assume that an evolution basis is specified, if that's + // not the case we must make the algorithm more generic + pids::pdg_mc_ids_to_evol(&tuples).unwrap() + }) + .collect() + }, + ); + + Ok(Self { + fac1, + info: OperatorSliceInfo { + lumi_id_types: pids::determine_lumi_id_types(&pids0), + fac0: operator.mu0 * operator.mu0, + pids0, + x0: metadata + .bases + .inputgrid + .unwrap_or_else(|| metadata.bases.xgrid.clone()), + fac1: 0.0, + pids1: metadata + .bases + .targetpids + .unwrap_or_else(|| BASES_V1_DEFAULT_PIDS.to_vec()), + x1: metadata + .bases + .targetgrid + .unwrap_or_else(|| metadata.bases.xgrid.clone()), + ren1: ren1.to_vec(), + alphas: alphas.to_vec(), + xir, + xif, + }, + archive: Archive::new(File::open(eko_path)?), + }) + } + + pub fn fac1(&self) -> Vec { + self.fac1.values().copied().collect() + } + + pub fn iter_mut(&mut self) -> EkoSlicesIter { + // UNWRAP: short of changing the return type we can't propagate the error, so we must + // panic here + EkoSlicesIter { + fac1: self.fac1.clone(), + info: self.info.clone(), + entries: self.archive.entries_with_seek().unwrap(), + } + } + } + + pub struct EkoSlicesIter<'a> { + fac1: HashMap, + info: OperatorSliceInfo, + entries: Entries<'a, File>, + } + + impl<'a> Iterator for EkoSlicesIter<'a> { + type Item = Result<(OperatorSliceInfo, Array4)>; + + fn next(&mut self) -> Option { + let mut fun = || { + while let Some(entry) = self.entries.next() { + let entry = entry?; + let path = entry.path()?; + + // here we're only interested in the operators themselves + if path.starts_with("./operators") + && (path.extension() == Some(OsStr::new("lz4"))) + && (path.with_extension("").extension() == Some(OsStr::new("npz"))) + { + // TODO: use let-else when available in MSRV + let file_stem = if let Some(file_stem) = path.with_extension("").file_stem() + { + file_stem.to_os_string() + } else { + continue; + }; + + let mut reader = BufReader::new(FrameDecoder::new(BufReader::new(entry))); + let mut buffer = Vec::new(); + io::copy(&mut reader, &mut buffer)?; + let mut npz = NpzReader::new(Cursor::new(buffer))?; + let operator: Array4 = npz.by_name("operator.npy")?; + + let mut info = self.info.clone(); + info.fac1 = self.fac1.get(&file_stem).copied().ok_or_else(|| anyhow!("file '{}.yaml' not found, could not determine the operator's factorization scale", file_stem.to_string_lossy()))?; + + return Ok(Some((info, operator))); + } + } + + Ok(None) + }; + + fun().transpose() + } + } + + // --- + pub fn read(eko: &Path) -> Result<(OperatorInfo, Array5)> { let mut archive = Archive::new(File::open(eko)?); From c8e85f82e2876afdb0144c6fa0370c5111976c9a Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 12 Oct 2023 15:08:57 +0200 Subject: [PATCH 08/40] Prepare for different EKO file formats --- pineappl_cli/src/evolve.rs | 161 ++++++++++++++++++++++++------------- 1 file changed, 105 insertions(+), 56 deletions(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 337e3f16..b4a9a656 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -68,8 +68,6 @@ mod eko { V2(MetadataV2), } - // --- EKO FORMAT INTRODUCED WITH EKO v0.13 --- - const BASES_V1_DEFAULT_PIDS: [i32; 14] = [22, -6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6]; #[derive(Deserialize)] @@ -96,22 +94,55 @@ mod eko { bases: BasesV1, } - pub struct EkoSlices { - fac1: HashMap, - info: OperatorSliceInfo, - archive: Archive, + pub enum EkoSlices { + V2 { + fac1: HashMap, + info: OperatorSliceInfo, + archive: Archive, + }, } impl EkoSlices { + /// Read the EKO at `eko_path` and return the contents of the `metadata.yaml` file + /// deserialized into a [`Metadata`] object. + fn read_metadata(eko_path: &Path) -> Result { + for entry in Archive::new(File::open(eko_path)?).entries_with_seek()? { + let entry = entry?; + let path = entry.path()?; + + if path.ends_with("metadata.yaml") { + return Ok(serde_yaml::from_reader(entry)?); + } + } + + Err(anyhow!("no file 'metadata.yaml' in EKO archive found")) + } + pub fn new( eko_path: &Path, ren1: &[f64], alphas: &[f64], xir: f64, xif: f64, + ) -> Result { + let metadata = Self::read_metadata(eko_path)?; + + match metadata { + Metadata::V0(v0) => todo!(), + Metadata::V1(v1) => todo!(), + Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas, xir, xif), + } + } + + fn with_v2( + metadata: MetadataV2, + eko_path: &Path, + ren1: &[f64], + alphas: &[f64], + xir: f64, + xif: f64, ) -> Result { let mut fac1 = HashMap::new(); - let mut metadata: Option = None; let mut operator: Option = None; for entry in Archive::new(File::open(eko_path)?).entries_with_seek()? { @@ -129,15 +160,11 @@ mod eko { let op_info: OperatorInfoV1 = serde_yaml::from_reader(entry)?; fac1.insert(file_stem, op_info.scale); - } else if path.as_os_str() == OsStr::new("./metadata.yaml") { - metadata = Some(serde_yaml::from_reader(entry)?); } else if path.as_os_str() == OsStr::new("./operator.yaml") { operator = Some(serde_yaml::from_reader(entry)?); } } - let metadata = - metadata.ok_or_else(|| anyhow!("no file 'metadata.yaml' in EKO archive found"))?; let operator = operator.ok_or_else(|| anyhow!("no file 'operator.yaml' in EKO archive found"))?; @@ -161,7 +188,7 @@ mod eko { }, ); - Ok(Self { + Ok(Self::V2 { fac1, info: OperatorSliceInfo { lumi_id_types: pids::determine_lumi_id_types(&pids0), @@ -189,66 +216,88 @@ mod eko { }) } + // TODO: we could make this a Cow<'_, [f64]> pub fn fac1(&self) -> Vec { - self.fac1.values().copied().collect() + match self { + Self::V2 { fac1, .. } => fac1.values().copied().collect(), + } } pub fn iter_mut(&mut self) -> EkoSlicesIter { - // UNWRAP: short of changing the return type we can't propagate the error, so we must - // panic here - EkoSlicesIter { - fac1: self.fac1.clone(), - info: self.info.clone(), - entries: self.archive.entries_with_seek().unwrap(), + match self { + Self::V2 { + fac1, + info, + archive, + } => { + EkoSlicesIter::V2 { + fac1: fac1.clone(), + info: info.clone(), + // UNWRAP: short of changing the return type of this method we can't + // propagate the error, so we must panic here + entries: archive.entries_with_seek().unwrap(), + } + } } } } - pub struct EkoSlicesIter<'a> { - fac1: HashMap, - info: OperatorSliceInfo, - entries: Entries<'a, File>, + pub enum EkoSlicesIter<'a> { + V2 { + fac1: HashMap, + info: OperatorSliceInfo, + entries: Entries<'a, File>, + }, } impl<'a> Iterator for EkoSlicesIter<'a> { type Item = Result<(OperatorSliceInfo, Array4)>; fn next(&mut self) -> Option { - let mut fun = || { - while let Some(entry) = self.entries.next() { - let entry = entry?; - let path = entry.path()?; - - // here we're only interested in the operators themselves - if path.starts_with("./operators") - && (path.extension() == Some(OsStr::new("lz4"))) - && (path.with_extension("").extension() == Some(OsStr::new("npz"))) - { - // TODO: use let-else when available in MSRV - let file_stem = if let Some(file_stem) = path.with_extension("").file_stem() - { - file_stem.to_os_string() - } else { - continue; - }; - - let mut reader = BufReader::new(FrameDecoder::new(BufReader::new(entry))); - let mut buffer = Vec::new(); - io::copy(&mut reader, &mut buffer)?; - let mut npz = NpzReader::new(Cursor::new(buffer))?; - let operator: Array4 = npz.by_name("operator.npy")?; - - let mut info = self.info.clone(); - info.fac1 = self.fac1.get(&file_stem).copied().ok_or_else(|| anyhow!("file '{}.yaml' not found, could not determine the operator's factorization scale", file_stem.to_string_lossy()))?; + match self { + Self::V2 { + fac1, + info, + entries, + } => { + let mut fun = || { + while let Some(entry) = entries.next() { + let entry = entry?; + let path = entry.path()?; + + // here we're only interested in the operators themselves + if path.starts_with("./operators") + && (path.extension() == Some(OsStr::new("lz4"))) + && (path.with_extension("").extension() == Some(OsStr::new("npz"))) + { + // TODO: use let-else when available in MSRV + let file_stem = + if let Some(file_stem) = path.with_extension("").file_stem() { + file_stem.to_os_string() + } else { + continue; + }; + + let mut reader = + BufReader::new(FrameDecoder::new(BufReader::new(entry))); + let mut buffer = Vec::new(); + io::copy(&mut reader, &mut buffer)?; + let mut npz = NpzReader::new(Cursor::new(buffer))?; + let operator: Array4 = npz.by_name("operator.npy")?; + + let mut info = info.clone(); + info.fac1 = fac1.get(&file_stem).copied().ok_or_else(|| anyhow!("file '{}.yaml' not found, could not determine the operator's factorization scale", file_stem.to_string_lossy()))?; + + return Ok(Some((info, operator))); + } + } + + Ok(None) + }; - return Ok(Some((info, operator))); - } + fun().transpose() } - - Ok(None) - }; - - fun().transpose() + } } } From 7935cea73e304a4248ed3638c2146956c69768ca Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 12 Oct 2023 15:15:14 +0200 Subject: [PATCH 09/40] Add support for v0 EKO file format --- pineappl_cli/src/evolve.rs | 77 +++++++++++++++++++++++++++++++++++++- 1 file changed, 75 insertions(+), 2 deletions(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index b4a9a656..554da5f9 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -17,7 +17,8 @@ mod eko { use base64::Engine; use either::Either; use lz4_flex::frame::FrameDecoder; - use ndarray::{Array4, Array5, Axis}; + use ndarray::iter::AxisIter; + use ndarray::{Array4, Array5, Axis, Ix4}; use ndarray_npy::{NpzReader, ReadNpyExt}; use pineappl::evolution::{OperatorInfo, OperatorSliceInfo}; use pineappl::pids; @@ -26,7 +27,9 @@ mod eko { use std::ffi::{OsStr, OsString}; use std::fs::File; use std::io::{self, BufReader, Cursor}; + use std::iter::Zip; use std::path::Path; + use std::slice::Iter; use tar::{Archive, Entries}; #[derive(Deserialize)] @@ -95,6 +98,11 @@ mod eko { } pub enum EkoSlices { + V0 { + fac1: Vec, + info: OperatorSliceInfo, + operator: Array5, + }, V2 { fac1: HashMap, info: OperatorSliceInfo, @@ -128,12 +136,53 @@ mod eko { let metadata = Self::read_metadata(eko_path)?; match metadata { - Metadata::V0(v0) => todo!(), + Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas, xir, xif), Metadata::V1(v1) => todo!(), Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas, xir, xif), } } + fn with_v0( + metadata: MetadataV0, + eko_path: &Path, + ren1: &[f64], + alphas: &[f64], + xir: f64, + xif: f64, + ) -> Result { + let mut operator = None; + + for entry in Archive::new(File::open(eko_path)?).entries_with_seek()? { + let entry = entry?; + let path = entry.path()?; + + if path.ends_with("operators.npy.lz4") { + operator = Some(Array5::read_npy(FrameDecoder::new(BufReader::new(entry)))?); + } + } + + let operator = + operator.ok_or_else(|| anyhow!("no file 'operator.yaml' in EKO archive found"))?; + + Ok(Self::V0 { + fac1: metadata.q2_grid, + info: OperatorSliceInfo { + lumi_id_types: pids::determine_lumi_id_types(&metadata.inputpids), + fac0: metadata.q2_ref, + pids0: metadata.inputpids, + x0: metadata.inputgrid, + fac1: 0.0, + pids1: metadata.targetpids, + x1: metadata.targetgrid, + ren1: ren1.to_vec(), + alphas: alphas.to_vec(), + xir, + xif, + }, + operator, + }) + } + fn with_v2( metadata: MetadataV2, eko_path: &Path, @@ -219,12 +268,21 @@ mod eko { // TODO: we could make this a Cow<'_, [f64]> pub fn fac1(&self) -> Vec { match self { + Self::V0 { fac1, .. } => fac1.clone(), Self::V2 { fac1, .. } => fac1.values().copied().collect(), } } pub fn iter_mut(&mut self) -> EkoSlicesIter { match self { + Self::V0 { + fac1, + info, + operator, + } => EkoSlicesIter::V0 { + info: info.clone(), + iter: fac1.iter().zip(operator.axis_iter(Axis(0))), + }, Self::V2 { fac1, info, @@ -243,6 +301,10 @@ mod eko { } pub enum EkoSlicesIter<'a> { + V0 { + info: OperatorSliceInfo, + iter: Zip, AxisIter<'a, f64, Ix4>>, + }, V2 { fac1: HashMap, info: OperatorSliceInfo, @@ -255,6 +317,17 @@ mod eko { fn next(&mut self) -> Option { match self { + Self::V0 { info, iter } => { + if let Some((fac1, operator)) = iter.next() { + let mut info = info.clone(); + info.fac1 = *fac1; + + // TODO: see if we can replace this some kind of Cow structure + Some(Ok((info, operator.to_owned()))) + } else { + None + } + } Self::V2 { fac1, info, From 5d1a4e9e2b25b2fbc69c503b2ebe34daf75b7fd7 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 12 Oct 2023 16:16:19 +0200 Subject: [PATCH 10/40] Add support for v1 EKO file format --- pineappl_cli/src/evolve.rs | 89 +++++++++++++++++++++++++++++++++++++- 1 file changed, 88 insertions(+), 1 deletion(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 554da5f9..3c9684b8 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -103,6 +103,7 @@ mod eko { info: OperatorSliceInfo, operator: Array5, }, + // V1 is a special case of V2 V2 { fac1: HashMap, info: OperatorSliceInfo, @@ -137,7 +138,7 @@ mod eko { match metadata { Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas, xir, xif), - Metadata::V1(v1) => todo!(), + Metadata::V1(v1) => Self::with_v1(v1, eko_path, ren1, alphas, xir, xif), Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas, xir, xif), } } @@ -183,6 +184,92 @@ mod eko { }) } + fn with_v1( + metadata: MetadataV1, + eko_path: &Path, + ren1: &[f64], + alphas: &[f64], + xir: f64, + xif: f64, + ) -> Result { + let mut fac1 = HashMap::new(); + let base64 = GeneralPurpose::new(&URL_SAFE, PAD); + + for entry in Archive::new(File::open(eko_path)?).entries_with_seek()? { + let entry = entry?; + let path = entry.path()?; + + if path.starts_with("./operators") + && (path.extension() == Some(OsStr::new("lz4"))) + && (path.with_extension("").extension() == Some(OsStr::new("npz"))) + { + // TODO: use let-else when available in MSRV + let file_stem = if let Some(file_stem) = path.with_extension("").file_stem() { + file_stem.to_os_string() + } else { + continue; + }; + + let bytes = base64.decode(file_stem.to_string_lossy().as_bytes())?; + let array: [u8; 8] = bytes.as_slice().try_into().unwrap(); + let scale = f64::from_le_bytes(array); + + fac1.insert(file_stem, scale); + } + } + + let pids0 = metadata.rotations.inputpids.map_or_else( + || metadata.rotations.pids.clone(), + |either| { + either.right_or_else(|basis| { + basis + .into_iter() + .map(|factors| { + let tuples: Vec<_> = metadata + .rotations + .pids + .iter() + .copied() + .zip(factors.into_iter()) + .collect(); + + // UNWRAP: we assume that an evolution basis is specified, if + // that's not the case we must make the algorithm more generic + pids::pdg_mc_ids_to_evol(&tuples).unwrap() + }) + .collect() + }) + }, + ); + + Ok(Self::V2 { + fac1, + info: OperatorSliceInfo { + lumi_id_types: pids::determine_lumi_id_types(&pids0), + fac0: metadata.mu20, + pids0, + x0: metadata + .rotations + .inputgrid + .unwrap_or_else(|| metadata.rotations.xgrid.clone()), + fac1: 0.0, + pids1: metadata + .rotations + .targetpids + .unwrap_or(metadata.rotations.pids), + x1: metadata + .rotations + .targetgrid + .unwrap_or(metadata.rotations.xgrid), + ren1: ren1.to_vec(), + alphas: alphas.to_vec(), + xir, + xif, + }, + archive: Archive::new(File::open(eko_path)?), + }) + } + fn with_v2( metadata: MetadataV2, eko_path: &Path, From 695e0d5674fdf6ccfe9c41ac0e32e594249c9d2b Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 12 Oct 2023 16:17:01 +0200 Subject: [PATCH 11/40] Add evolution for single-hadron initial-state grids --- pineappl/src/evolution.rs | 148 ++++++++++++++++++++++++++++++++++++++ pineappl/src/grid.rs | 3 +- 2 files changed, 149 insertions(+), 2 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 695681f0..deeb57f1 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -687,6 +687,154 @@ pub(crate) fn evolve_with_one( )) } +pub(crate) fn evolve_slice_with_one( + grid: &Grid, + operator: &ArrayView4, + info: &OperatorSliceInfo, + order_mask: &[bool], +) -> Result<(Array3, Vec), GridError> { + let op5 = operator.insert_axis(Axis(0)); + let inf = OperatorInfo { + fac0: info.fac0, + pids0: info.pids0.clone(), + x0: info.x0.clone(), + fac1: vec![info.fac1], + pids1: info.pids1.clone(), + x1: info.x1.clone(), + ren1: info.ren1.clone(), + alphas: info.alphas.clone(), + xir: info.xir, + xif: info.xif, + lumi_id_types: info.lumi_id_types.clone(), + }; + let gluon_has_pid_zero = gluon_has_pid_zero(grid); + let has_pdf1 = grid.has_pdf1(); + + let (pid_indices, pids) = pids(&op5, &inf, gluon_has_pid_zero, &|pid| { + grid.lumi() + .iter() + .flat_map(LumiEntry::entry) + .any(|&(a, b, _)| if has_pdf1 { a } else { b } == pid) + })?; + + let lumi0 = lumi0_with_one(&pids); + let mut sub_fk_tables = Vec::with_capacity(grid.bin_info().bins() * lumi0.len()); + let new_axis = if has_pdf1 { 2 } else { 1 }; + + let mut last_x1 = Vec::new(); + let mut ops = Vec::new(); + + for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { + let mut tables = vec![Array1::zeros(info.x0.len()); lumi0.len()]; + + for (lumi1, subgrids_o) in subgrids_ol.axis_iter(Axis(1)).enumerate() { + let (x1_a, x1_b, array) = + ndarray_from_subgrid_orders_slice(info, &subgrids_o, grid.orders(), order_mask)?; + + let x1 = if has_pdf1 { x1_a } else { x1_b }; + + if x1.is_empty() { + continue; + } + + if (last_x1.len() != x1.len()) + || last_x1 + .iter() + .zip(x1.iter()) + .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) + { + ops = operators(&op5, &inf, &inf.fac1, &pid_indices, &x1)?; + last_x1 = x1; + } + + // TODO: get rid of array-index access + for (&pid1, &factor) in + grid.lumi()[lumi1].entry().iter().map( + |(a, b, f)| { + if has_pdf1 { + (a, f) + } else { + (b, f) + } + }, + ) + { + for (fk_table, op) in + lumi0 + .iter() + .zip(tables.iter_mut()) + .filter_map(|(&pid0, fk_table)| { + pids.iter() + .zip(ops.iter()) + .find_map(|(&(p0, p1), op)| { + (p0 == pid0 && p1 == pid1).then_some(op) + }) + .map(|op| (fk_table, op)) + }) + { + let mut result = Array1::zeros(info.x0.len()); + + for imu2 in 0..array.dim().0 { + let op = op.index_axis(Axis(0), imu2); + + result += &op.dot( + &array + .index_axis(Axis(0), imu2) + .index_axis(Axis(new_axis - 1), 0), + ); + } + + fk_table.scaled_add(factor, &result); + } + } + } + + sub_fk_tables.extend(tables.into_iter().map(|table| { + ImportOnlySubgridV2::new( + SparseArray3::from_ndarray( + table + .insert_axis(Axis(0)) + .insert_axis(Axis(new_axis)) + .view(), + 0, + 1, + ), + vec![Mu2 { + // TODO: FK tables don't depend on the renormalization scale + //ren: -1.0, + ren: info.fac0, + fac: info.fac0, + }], + if has_pdf1 { info.x0.clone() } else { vec![1.0] }, + if has_pdf1 { vec![1.0] } else { info.x0.clone() }, + ) + .into() + })); + } + + let pid = if has_pdf1 { + grid.initial_state_2() + } else { + grid.initial_state_1() + }; + + Ok(( + Array1::from_iter(sub_fk_tables) + .into_shape((1, grid.bin_info().bins(), lumi0.len())) + .unwrap(), + lumi0 + .iter() + .map(|&a| { + lumi_entry![ + if has_pdf1 { a } else { pid }, + if has_pdf1 { pid } else { a }, + 1.0 + ] + }) + .collect(), + )) +} + pub(crate) fn evolve_with_two( grid: &Grid, operator: &ArrayView5, diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 023761cb..757c3bea 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -1979,8 +1979,7 @@ impl Grid { let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { evolution::evolve_slice_with_two(self, &operator, info, order_mask) } else { - todo!() - //evolution::evolve_slice_with_one(self, &operator, info, order_mask) + evolution::evolve_slice_with_one(self, &operator, info, order_mask) }?; let mut grid = Self { From b94ceeddeb1b1088873cd574587db23e9a416195 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 12 Oct 2023 16:20:25 +0200 Subject: [PATCH 12/40] Replace `evolve` with `evolve_slice` in `pineappl evolve` --- pineappl_cli/src/evolve.rs | 218 ++++++------------------------------- 1 file changed, 34 insertions(+), 184 deletions(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 3c9684b8..a2c4c1e8 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -10,7 +10,7 @@ use std::process::ExitCode; #[cfg(feature = "evolve")] mod eko { - use anyhow::{anyhow, bail, Result}; + use anyhow::{anyhow, Result}; use base64::alphabet::URL_SAFE; use base64::engine::general_purpose::PAD; use base64::engine::GeneralPurpose; @@ -20,7 +20,7 @@ mod eko { use ndarray::iter::AxisIter; use ndarray::{Array4, Array5, Axis, Ix4}; use ndarray_npy::{NpzReader, ReadNpyExt}; - use pineappl::evolution::{OperatorInfo, OperatorSliceInfo}; + use pineappl::evolution::OperatorSliceInfo; use pineappl::pids; use serde::Deserialize; use std::collections::HashMap; @@ -460,178 +460,6 @@ mod eko { } } } - - // --- - - pub fn read(eko: &Path) -> Result<(OperatorInfo, Array5)> { - let mut archive = Archive::new(File::open(eko)?); - - let mut metadata = None; - let mut operator = None; - let mut operators = Vec::new(); - let mut fac1 = Vec::new(); - let mut fac0 = None; - let mut fac1_order = Vec::new(); - let mut fac1_translator = HashMap::new(); - - let base64 = GeneralPurpose::new(&URL_SAFE, PAD); - - for entry in archive.entries()? { - let file = entry?; - let path = file.header().path()?; - - if let Some(file_name) = path.file_name() { - // TODO: get rid of the unwrap - match file_name.to_str().unwrap() { - "metadata.yaml" => metadata = Some(serde_yaml::from_reader(file)?), - "operator.yaml" => { - let operator: OperatorV1 = serde_yaml::from_reader(file)?; - fac0 = Some(operator.mu0 * operator.mu0); - } - "operators.npy.lz4" => { - operator = Some(Array5::::read_npy(FrameDecoder::new( - BufReader::new(file), - ))?); - } - x if x.ends_with(".npz.lz4") => { - fac1_order.push(x.strip_suffix(".npz.lz4").unwrap().to_string()); - - let name = x.strip_suffix(".npz.lz4").unwrap(); - let bytes = base64.decode(name.as_bytes())?; - let array: [u8; 8] = bytes.as_slice().try_into().unwrap(); - let muf2 = f64::from_le_bytes(array); - - let mut reader = BufReader::new(FrameDecoder::new(BufReader::new(file))); - let mut buffer = Vec::new(); - std::io::copy(&mut reader, &mut buffer)?; - let mut npz = NpzReader::new(std::io::Cursor::new(buffer))?; - let operator: Array4 = npz.by_name("operator.npy")?; - - fac1.push(muf2); - operators.push(operator); - } - x if x.ends_with("=.yaml") => { - let name = x.strip_suffix(".yaml").unwrap().to_string(); - let op_info: OperatorInfoV1 = serde_yaml::from_reader(file)?; - fac1_translator.insert(name, op_info.scale); - } - _ => {} - } - } - } - - if !operators.is_empty() { - let ops: Vec<_> = operators.iter().map(Array4::view).collect(); - operator = Some(ndarray::stack(Axis(0), &ops).unwrap()); - } - - if !fac1_translator.is_empty() { - for (fac1, fac1_name) in fac1.iter_mut().zip(&fac1_order) { - *fac1 = fac1_translator[fac1_name]; - } - } - - let mut info = match metadata { - Some(Metadata::V0(metadata)) => OperatorInfo { - fac1: metadata.q2_grid.clone(), - pids0: metadata.inputpids, - x0: metadata.inputgrid, - pids1: metadata.targetpids, - x1: metadata.targetgrid, - fac0: metadata.q2_ref, - ren1: vec![], - alphas: vec![], - xir: 1.0, - xif: 1.0, - lumi_id_types: String::new(), - }, - Some(Metadata::V1(metadata)) => OperatorInfo { - fac1, - pids0: metadata.rotations.inputpids.map_or_else( - || metadata.rotations.pids.clone(), - |either| { - either.right_or_else(|basis| { - basis - .into_iter() - .map(|factors| { - let tuples: Vec<_> = metadata - .rotations - .pids - .iter() - .copied() - .zip(factors.into_iter()) - .collect(); - - pids::pdg_mc_ids_to_evol(&tuples).unwrap() - }) - .collect() - }) - }, - ), - x0: metadata - .rotations - .inputgrid - .unwrap_or_else(|| metadata.rotations.xgrid.clone()), - pids1: metadata - .rotations - .targetpids - .unwrap_or(metadata.rotations.pids), - x1: metadata - .rotations - .targetgrid - .unwrap_or(metadata.rotations.xgrid), - fac0: metadata.mu20, - ren1: vec![], - alphas: vec![], - xir: 1.0, - xif: 1.0, - lumi_id_types: String::new(), - }, - Some(Metadata::V2(metadata)) => OperatorInfo { - fac1, - pids0: metadata.bases.inputpids.map_or_else( - || BASES_V1_DEFAULT_PIDS.to_vec(), - |basis| { - basis - .into_iter() - .map(|factors| { - let tuples: Vec<_> = BASES_V1_DEFAULT_PIDS - .iter() - .copied() - .zip(factors.into_iter()) - .collect(); - - pids::pdg_mc_ids_to_evol(&tuples).unwrap() - }) - .collect() - }, - ), - x0: metadata - .bases - .inputgrid - .unwrap_or_else(|| metadata.bases.xgrid.clone()), - pids1: metadata - .bases - .targetpids - .unwrap_or_else(|| BASES_V1_DEFAULT_PIDS.to_vec()), - x1: metadata - .bases - .targetgrid - .unwrap_or_else(|| metadata.bases.xgrid.clone()), - fac0: fac0.unwrap(), - ren1: vec![], - alphas: vec![], - xir: 1.0, - xif: 1.0, - lumi_id_types: String::new(), - }, - None => bail!("no `metadata.yaml` file found"), - }; - - info.lumi_id_types = pids::determine_lumi_id_types(&info.pids0); - - Ok((info, operator.unwrap())) - } } #[cfg(feature = "evolve")] @@ -643,11 +471,10 @@ fn evolve_grid( xir: f64, xif: f64, ) -> Result { + use eko::EkoSlices; + use float_cmp::approx_eq; use pineappl::subgrid::{Mu2, Subgrid}; - let (mut info, operator) = eko::read(eko)?; - - // TODO: the following should probably be a method of `Grid` let mut ren1: Vec<_> = grid .subgrids() .iter() @@ -664,12 +491,7 @@ fn evolve_grid( let ren1 = ren1; let alphas: Vec<_> = ren1.iter().map(|&mur2| pdf.alphas_q2(mur2)).collect(); - info.ren1 = ren1; - info.alphas = alphas; - info.xir = xir; - info.xif = xif; - - let orders: Vec<_> = grid + let order_mask: Vec<_> = grid .orders() .iter() .map(|order| { @@ -680,7 +502,35 @@ fn evolve_grid( }) .collect(); - Ok(grid.evolve(operator.view(), &info, &orders)?) + let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas, xir, xif)?; + let grid_muf2 = grid.evolve_info(&order_mask).fac1; + + // check that every factorization-scale value in the grid has an operator + assert!(!grid_muf2 + .iter() + .any(|&muf2| !eko_slices.fac1().iter().any(|&op_muf2| approx_eq!( + f64, + muf2, + op_muf2, + ulps = 64 + )))); + + let mut lhs: Option = None; + + for item in eko_slices.iter_mut() { + let (info, operator) = item?; + let rhs = grid + .evolve_slice(operator.view(), &info, &order_mask)? + .into_grid(); + + if let Some(lhs) = &mut lhs { + lhs.merge(rhs)?; + } else { + lhs = Some(rhs); + } + } + + Ok(FkTable::try_from(lhs.unwrap()).unwrap()) } #[cfg(not(feature = "evolve"))] From 3bc3beb03cd8a5358e5042640f8f398469cdc8b7 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Fri, 13 Oct 2023 12:58:51 +0200 Subject: [PATCH 13/40] Add missing UNWRAP comment --- pineappl_cli/src/evolve.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index a2c4c1e8..1aee301e 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -211,6 +211,7 @@ mod eko { }; let bytes = base64.decode(file_stem.to_string_lossy().as_bytes())?; + // UNWRAP: we assume that the filenames represent exactly 8 bytes let array: [u8; 8] = bytes.as_slice().try_into().unwrap(); let scale = f64::from_le_bytes(array); From b7872a77aa6530a086a28f8e29cf1a02c0957b8f Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Fri, 13 Oct 2023 12:59:43 +0200 Subject: [PATCH 14/40] Fix E906 evolution integration test --- pineappl_cli/tests/evolve.rs | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/pineappl_cli/tests/evolve.rs b/pineappl_cli/tests/evolve.rs index e69d2769..42a4b717 100644 --- a/pineappl_cli/tests/evolve.rs +++ b/pineappl_cli/tests/evolve.rs @@ -26,11 +26,8 @@ Options: const E906NLO_BIN_00_STR: &str = "b Grid FkTable rel. diff -+------------+------------+------------- 0 1.0659807e-1 1.0657904e-1 -1.7851986e-4 -1 1.0659807e-1 1.0657904e-1 -1.7851986e-4 -2 1.0659807e-1 1.0657904e-1 -1.7851986e-4 -3 1.0659807e-1 1.0657904e-1 -1.7851986e-4 -4 3.2698655e0 3.2711890e0 4.0477586e-4 -5 1.6039253e0 1.6047566e0 5.1825508e-4 +1 3.2698655e0 3.2711890e0 4.0477586e-4 +2 1.6039253e0 1.6047566e0 5.1825508e-4 "; const LHCB_DY_8TEV_STR: &str = "b Grid FkTable rel. diff @@ -242,14 +239,28 @@ fn lhcb_wp_7tev_v2() { #[test] fn e906nlo_bin_00() { + let input = NamedTempFile::new("E906nlo_bin_00_unique_bin_limits.pineappl.lz4").unwrap(); let output = NamedTempFile::new("fktable2.lz4").unwrap(); + // we need to delete bins with the same bin limits for `Grid::merge` to work properly + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "write", + "--delete-bins=1-3", + "../test-data/E906nlo_bin_00.pineappl.lz4", + input.path().to_str().unwrap(), + ]) + .assert() + .success() + .stdout(""); + Command::cargo_bin("pineappl") .unwrap() .args([ "--silence-lhapdf", "evolve", - "../test-data/E906nlo_bin_00.pineappl.lz4", + input.path().to_str().unwrap(), "../test-data/E906nlo_bin_00.tar", output.path().to_str().unwrap(), "NNPDF40_nlo_as_01180", From 6ab3a6b3964a6cf116e1acf899926537194d7ce1 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Fri, 13 Oct 2023 14:35:45 +0200 Subject: [PATCH 15/40] Fix some clippy warnings --- pineappl_cli/src/evolve.rs | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 1aee301e..a993f93d 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -231,7 +231,7 @@ mod eko { .pids .iter() .copied() - .zip(factors.into_iter()) + .zip(factors) .collect(); // UNWRAP: we assume that an evolution basis is specified, if @@ -311,11 +311,8 @@ mod eko { basis .into_iter() .map(|factors| { - let tuples: Vec<_> = BASES_V1_DEFAULT_PIDS - .iter() - .copied() - .zip(factors.into_iter()) - .collect(); + let tuples: Vec<_> = + BASES_V1_DEFAULT_PIDS.iter().copied().zip(factors).collect(); // UNWRAP: we assume that an evolution basis is specified, if that's // not the case we must make the algorithm more generic @@ -388,6 +385,15 @@ mod eko { } } + impl<'a> IntoIterator for &'a mut EkoSlices { + type Item = Result<(OperatorSliceInfo, Array4)>; + type IntoIter = EkoSlicesIter<'a>; + + fn into_iter(self) -> Self::IntoIter { + self.iter_mut() + } + } + pub enum EkoSlicesIter<'a> { V0 { info: OperatorSliceInfo, @@ -421,8 +427,8 @@ mod eko { info, entries, } => { - let mut fun = || { - while let Some(entry) = entries.next() { + let fun = || { + for entry in entries { let entry = entry?; let path = entry.path()?; From 4e637e4fdfae8229179fd27c1c644f93f8e58cfd Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Fri, 13 Oct 2023 19:29:35 +0200 Subject: [PATCH 16/40] Simplify slice-based evolution code --- pineappl/src/evolution.rs | 158 +++++++++++++++++++++++--------------- 1 file changed, 96 insertions(+), 62 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index deeb57f1..0ff1ed11 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -209,6 +209,54 @@ pub(crate) fn pids( Ok((pid_indices, pids)) } +pub(crate) fn pid_slices( + operator: &ArrayView4, + info: &OperatorSliceInfo, + gluon_has_pid_zero: bool, + pid1_nonzero: &dyn Fn(i32) -> bool, +) -> Result<(Pid01IndexTuples, Pid01Tuples), GridError> { + // list of all non-zero PID indices + let pid_indices: Vec<_> = (0..operator.dim().2) + .cartesian_product(0..operator.dim().0) + .filter(|&(pid0_idx, pid1_idx)| { + // 1) at least one element of the operator must be non-zero, and 2) the pid must be + // contained in the lumi somewhere + operator + .slice(s![pid1_idx, .., pid0_idx, ..]) + .iter() + .any(|&value| value != 0.0) + && pid1_nonzero(if gluon_has_pid_zero && info.pids1[pid1_idx] == 21 { + 0 + } else { + info.pids1[pid1_idx] + }) + }) + .collect(); + + if pid_indices.is_empty() { + return Err(GridError::EvolutionFailure( + "no non-zero operator found; result would be an empty FkTable".to_string(), + )); + } + + // list of all non-zero (pid0, pid1) combinations + let pids = pid_indices + .iter() + .map(|&(pid0_idx, pid1_idx)| { + ( + info.pids0[pid0_idx], + if gluon_has_pid_zero && info.pids1[pid1_idx] == 21 { + 0 + } else { + info.pids1[pid1_idx] + }, + ) + }) + .collect(); + + Ok((pid_indices, pids)) +} + pub(crate) fn lumi0_with_one(pids: &[(i32, i32)]) -> Vec { let mut pids0: Vec<_> = pids.iter().map(|&(pid0, _)| pid0).collect(); pids0.sort_unstable(); @@ -284,6 +332,42 @@ pub(crate) fn operators( Ok(operators) } +pub(crate) fn operator_slices( + operator: &ArrayView4, + info: &OperatorSliceInfo, + pid_indices: &[(usize, usize)], + x1: &[f64], +) -> Result>, GridError> { + // permutation between the grid x values and the operator x1 values + let x1_indices: Vec<_> = x1 + .iter() + .map(|&x1p| { + info.x1 + .iter() + .position(|&x1| approx_eq!(f64, x1p, x1, ulps = EVOLUTION_TOL_ULPS)) + .ok_or_else(|| { + GridError::EvolutionFailure(format!("no operator for x = {x1p} found")) + }) + }) + // TODO: use `try_collect` once stabilized + .collect::>()?; + + // create the corresponding operators accessible in the form [muf2, x0, x1] + let operators: Vec<_> = pid_indices + .iter() + .map(|&(pid0_idx, pid1_idx)| { + operator + .slice(s![pid1_idx, .., pid0_idx, ..]) + .select(Axis(0), &x1_indices) + .permuted_axes([1, 0]) + .as_standard_layout() + .into_owned() + }) + .collect(); + + Ok(operators) +} + type Fac1X1aX1bOp3Tuple = (Vec, Vec, Vec, Array3); pub(crate) fn ndarray_from_subgrid_orders( @@ -425,14 +509,14 @@ pub(crate) fn ndarray_from_subgrid_orders( Ok((fac1, x1_a, x1_b, array)) } -type X1aX1bOp3Tuple = (Vec, Vec, Array3); +type X1aX1bOp2Tuple = (Vec, Vec, Array2); fn ndarray_from_subgrid_orders_slice( info: &OperatorSliceInfo, subgrids: &ArrayView1, orders: &[Order], order_mask: &[bool], -) -> Result { +) -> Result { // TODO: skip empty subgrids let fac1 = info.xif * info.xif * info.fac1; @@ -454,7 +538,7 @@ fn ndarray_from_subgrid_orders_slice( x1_b.sort_by(f64::total_cmp); x1_b.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLUTION_TOL_ULPS)); - let mut array = Array3::::zeros((1, x1_a.len(), x1_b.len())); + let mut array = Array2::::zeros((x1_a.len(), x1_b.len())); // add subgrids for different orders, but the same bin and lumi, using the right // couplings @@ -539,7 +623,7 @@ fn ndarray_from_subgrid_orders_slice( ))); }; - array[[0, xa_indices[ix1], xb_indices[ix2]]] += als * logs * value; + array[[xa_indices[ix1], xb_indices[ix2]]] += als * logs * value; } } @@ -693,24 +777,10 @@ pub(crate) fn evolve_slice_with_one( info: &OperatorSliceInfo, order_mask: &[bool], ) -> Result<(Array3, Vec), GridError> { - let op5 = operator.insert_axis(Axis(0)); - let inf = OperatorInfo { - fac0: info.fac0, - pids0: info.pids0.clone(), - x0: info.x0.clone(), - fac1: vec![info.fac1], - pids1: info.pids1.clone(), - x1: info.x1.clone(), - ren1: info.ren1.clone(), - alphas: info.alphas.clone(), - xir: info.xir, - xif: info.xif, - lumi_id_types: info.lumi_id_types.clone(), - }; let gluon_has_pid_zero = gluon_has_pid_zero(grid); let has_pdf1 = grid.has_pdf1(); - let (pid_indices, pids) = pids(&op5, &inf, gluon_has_pid_zero, &|pid| { + let (pid_indices, pids) = pid_slices(operator, info, gluon_has_pid_zero, &|pid| { grid.lumi() .iter() .flat_map(LumiEntry::entry) @@ -743,7 +813,7 @@ pub(crate) fn evolve_slice_with_one( .zip(x1.iter()) .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) { - ops = operators(&op5, &inf, &inf.fac1, &pid_indices, &x1)?; + ops = operator_slices(operator, info, &pid_indices, &x1)?; last_x1 = x1; } @@ -772,19 +842,7 @@ pub(crate) fn evolve_slice_with_one( .map(|op| (fk_table, op)) }) { - let mut result = Array1::zeros(info.x0.len()); - - for imu2 in 0..array.dim().0 { - let op = op.index_axis(Axis(0), imu2); - - result += &op.dot( - &array - .index_axis(Axis(0), imu2) - .index_axis(Axis(new_axis - 1), 0), - ); - } - - fk_table.scaled_add(factor, &result); + fk_table.scaled_add(factor, &op.dot(&array.index_axis(Axis(new_axis - 1), 0))); } } } @@ -967,29 +1025,15 @@ pub(crate) fn evolve_slice_with_two( info: &OperatorSliceInfo, order_mask: &[bool], ) -> Result<(Array3, Vec), GridError> { - let op5 = operator.insert_axis(Axis(0)); - let inf = OperatorInfo { - fac0: info.fac0, - pids0: info.pids0.clone(), - x0: info.x0.clone(), - fac1: vec![info.fac1], - pids1: info.pids1.clone(), - x1: info.x1.clone(), - ren1: info.ren1.clone(), - alphas: info.alphas.clone(), - xir: info.xir, - xif: info.xif, - lumi_id_types: info.lumi_id_types.clone(), - }; let gluon_has_pid_zero = gluon_has_pid_zero(grid); - let (pid_indices_a, pids_a) = pids(&op5, &inf, gluon_has_pid_zero, &|pid1| { + let (pid_indices_a, pids_a) = pid_slices(operator, info, gluon_has_pid_zero, &|pid1| { grid.lumi() .iter() .flat_map(LumiEntry::entry) .any(|&(a, _, _)| a == pid1) })?; - let (pid_indices_b, pids_b) = pids(&op5, &inf, gluon_has_pid_zero, &|pid1| { + let (pid_indices_b, pids_b) = pid_slices(operator, info, gluon_has_pid_zero, &|pid1| { grid.lumi() .iter() .flat_map(LumiEntry::entry) @@ -1017,7 +1061,7 @@ pub(crate) fn evolve_slice_with_two( .zip(x1_a.iter()) .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) { - operators_a = operators(&op5, &inf, &inf.fac1, &pid_indices_a, &x1_a)?; + operators_a = operator_slices(operator, info, &pid_indices_a, &x1_a)?; last_x1a = x1_a; } @@ -1027,7 +1071,7 @@ pub(crate) fn evolve_slice_with_two( .zip(x1_b.iter()) .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) { - operators_b = operators(&op5, &inf, &inf.fac1, &pid_indices_b, &x1_b)?; + operators_b = operator_slices(operator, info, &pid_indices_b, &x1_b)?; last_x1b = x1_b; } @@ -1049,17 +1093,7 @@ pub(crate) fn evolve_slice_with_two( .map(|(opa, opb)| (fk_table, opa, opb)) }) { - let mut result = Array2::zeros((info.x0.len(), info.x0.len())); - - for imu2 in 0..array.dim().0 { - let opa = opa.index_axis(Axis(0), imu2); - let opb = opb.index_axis(Axis(0), imu2); - let arr = array.index_axis(Axis(0), imu2); - - result += &opa.dot(&arr.dot(&opb.t())); - } - - fk_table.scaled_add(factor, &result); + fk_table.scaled_add(factor, &opa.dot(&array.dot(&opb.t()))); } } } From 8dc1eb4a1800cffcefb7ed6c544a21a627850468 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 14 Oct 2023 12:32:33 +0200 Subject: [PATCH 17/40] Replace index-access with iterators --- pineappl/src/evolution.rs | 44 +++++++++++++-------------------------- 1 file changed, 14 insertions(+), 30 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 0ff1ed11..a9225678 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -657,7 +657,7 @@ pub(crate) fn evolve_with_one( for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { let mut tables = vec![Array1::zeros(info.x0.len()); lumi0.len()]; - for (lumi1, subgrids_o) in subgrids_ol.axis_iter(Axis(1)).enumerate() { + for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { let (fac1, x1_a, x1_b, array) = ndarray_from_subgrid_orders(info, &subgrids_o, grid.orders(), order_mask)?; @@ -683,17 +683,10 @@ pub(crate) fn evolve_with_one( last_x1 = x1; } - // TODO: get rid of array-index access - for (&pid1, &factor) in - grid.lumi()[lumi1].entry().iter().map( - |(a, b, f)| { - if has_pdf1 { - (a, f) - } else { - (b, f) - } - }, - ) + for (&pid1, &factor) in lumi1 + .entry() + .iter() + .map(|(a, b, f)| if has_pdf1 { (a, f) } else { (b, f) }) { for (fk_table, op) in lumi0 @@ -797,7 +790,7 @@ pub(crate) fn evolve_slice_with_one( for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { let mut tables = vec![Array1::zeros(info.x0.len()); lumi0.len()]; - for (lumi1, subgrids_o) in subgrids_ol.axis_iter(Axis(1)).enumerate() { + for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { let (x1_a, x1_b, array) = ndarray_from_subgrid_orders_slice(info, &subgrids_o, grid.orders(), order_mask)?; @@ -817,17 +810,10 @@ pub(crate) fn evolve_slice_with_one( last_x1 = x1; } - // TODO: get rid of array-index access - for (&pid1, &factor) in - grid.lumi()[lumi1].entry().iter().map( - |(a, b, f)| { - if has_pdf1 { - (a, f) - } else { - (b, f) - } - }, - ) + for (&pid1, &factor) in lumi1 + .entry() + .iter() + .map(|(a, b, f)| if has_pdf1 { (a, f) } else { (b, f) }) { for (fk_table, op) in lumi0 @@ -926,7 +912,7 @@ pub(crate) fn evolve_with_two( for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { let mut tables = vec![Array2::zeros((info.x0.len(), info.x0.len())); lumi0.len()]; - for (lumi1, subgrids_o) in subgrids_ol.axis_iter(Axis(1)).enumerate() { + for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { let (fac1, x1_a, x1_b, array) = ndarray_from_subgrid_orders(info, &subgrids_o, grid.orders(), order_mask)?; @@ -962,8 +948,7 @@ pub(crate) fn evolve_with_two( last_fac1 = fac1; }; - // TODO: get rid of array-index access - for &(pida1, pidb1, factor) in grid.lumi()[lumi1].entry() { + for &(pida1, pidb1, factor) in lumi1.entry() { for (fk_table, opa, opb) in lumi0 .iter() @@ -1051,7 +1036,7 @@ pub(crate) fn evolve_slice_with_two( for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { let mut tables = vec![Array2::zeros((info.x0.len(), info.x0.len())); lumi0.len()]; - for (lumi1, subgrids_o) in subgrids_ol.axis_iter(Axis(1)).enumerate() { + for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { let (x1_a, x1_b, array) = ndarray_from_subgrid_orders_slice(info, &subgrids_o, grid.orders(), order_mask)?; @@ -1075,8 +1060,7 @@ pub(crate) fn evolve_slice_with_two( last_x1b = x1_b; } - // TODO: get rid of array-index access - for &(pida1, pidb1, factor) in grid.lumi()[lumi1].entry() { + for &(pida1, pidb1, factor) in lumi1.entry() { for (fk_table, opa, opb) in lumi0 .iter() From ef962fa3896f6b8bb0c73b40794ff2e4969308e9 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 14 Oct 2023 13:21:20 +0200 Subject: [PATCH 18/40] Improve evolution performance for double-hadronic grids Instead of using a cartesian product of two PIDs which is subsequently filtered use two iterators of PIDs first filtered and subsequently zipped. This is a change of O(n^2) to O(n) which makes a big difference if n is large enough. This is the case for dijet grids, for example. --- pineappl/src/evolution.rs | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index a9225678..a3dc0217 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -957,11 +957,14 @@ pub(crate) fn evolve_with_two( pids_a .iter() .zip(operators_a.iter()) - .cartesian_product(pids_b.iter().zip(operators_b.iter())) - .find_map(|((&(pa0, pa1), opa), (&(pb0, pb1), opb))| { - (pa0 == pida0 && pa1 == pida1 && pb0 == pidb0 && pb1 == pidb1) - .then_some((opa, opb)) + .find_map(|(&(pa0, pa1), opa)| { + (pa0 == pida0 && pa1 == pida1).then_some(opa) }) + .zip(pids_b.iter().zip(operators_b.iter()).find_map( + |(&(pb0, pb1), opb)| { + (pb0 == pidb0 && pb1 == pidb1).then_some(opb) + }, + )) .map(|(opa, opb)| (fk_table, opa, opb)) }) { @@ -1069,11 +1072,14 @@ pub(crate) fn evolve_slice_with_two( pids_a .iter() .zip(operators_a.iter()) - .cartesian_product(pids_b.iter().zip(operators_b.iter())) - .find_map(|((&(pa0, pa1), opa), (&(pb0, pb1), opb))| { - (pa0 == pida0 && pa1 == pida1 && pb0 == pidb0 && pb1 == pidb1) - .then_some((opa, opb)) + .find_map(|(&(pa0, pa1), opa)| { + (pa0 == pida0 && pa1 == pida1).then_some(opa) }) + .zip(pids_b.iter().zip(operators_b.iter()).find_map( + |(&(pb0, pb1), opb)| { + (pb0 == pidb0 && pb1 == pidb1).then_some(opb) + }, + )) .map(|(opa, opb)| (fk_table, opa, opb)) }) { From d9cb34f3493919245067e031efd0c3b89fb5822a Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 17 Oct 2023 15:15:29 +0200 Subject: [PATCH 19/40] Make evolution functions private --- pineappl/src/evolution.rs | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index a3dc0217..f7ddc23e 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -161,7 +161,7 @@ fn gluon_has_pid_zero(grid: &Grid) -> bool { type Pid01IndexTuples = Vec<(usize, usize)>; type Pid01Tuples = Vec<(i32, i32)>; -pub(crate) fn pids( +fn pids( operator: &ArrayView5, info: &OperatorInfo, gluon_has_pid_zero: bool, @@ -209,7 +209,7 @@ pub(crate) fn pids( Ok((pid_indices, pids)) } -pub(crate) fn pid_slices( +fn pid_slices( operator: &ArrayView4, info: &OperatorSliceInfo, gluon_has_pid_zero: bool, @@ -257,7 +257,7 @@ pub(crate) fn pid_slices( Ok((pid_indices, pids)) } -pub(crate) fn lumi0_with_one(pids: &[(i32, i32)]) -> Vec { +fn lumi0_with_one(pids: &[(i32, i32)]) -> Vec { let mut pids0: Vec<_> = pids.iter().map(|&(pid0, _)| pid0).collect(); pids0.sort_unstable(); pids0.dedup(); @@ -265,7 +265,7 @@ pub(crate) fn lumi0_with_one(pids: &[(i32, i32)]) -> Vec { pids0 } -pub(crate) fn lumi0_with_two(pids_a: &[(i32, i32)], pids_b: &[(i32, i32)]) -> Vec<(i32, i32)> { +fn lumi0_with_two(pids_a: &[(i32, i32)], pids_b: &[(i32, i32)]) -> Vec<(i32, i32)> { let mut pids0_a: Vec<_> = pids_a.iter().map(|&(pid0, _)| pid0).collect(); pids0_a.sort_unstable(); pids0_a.dedup(); @@ -280,7 +280,7 @@ pub(crate) fn lumi0_with_two(pids_a: &[(i32, i32)], pids_b: &[(i32, i32)]) -> Ve .collect() } -pub(crate) fn operators( +fn operators( operator: &ArrayView5, info: &OperatorInfo, fac1: &[f64], @@ -332,7 +332,7 @@ pub(crate) fn operators( Ok(operators) } -pub(crate) fn operator_slices( +fn operator_slices( operator: &ArrayView4, info: &OperatorSliceInfo, pid_indices: &[(usize, usize)], @@ -370,7 +370,7 @@ pub(crate) fn operator_slices( type Fac1X1aX1bOp3Tuple = (Vec, Vec, Vec, Array3); -pub(crate) fn ndarray_from_subgrid_orders( +fn ndarray_from_subgrid_orders( info: &OperatorInfo, subgrids: &ArrayView1, orders: &[Order], From 11723aa0fe688727818d1917dded3ecb90f1d6f8 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 18 Oct 2023 14:39:08 +0200 Subject: [PATCH 20/40] Use `reversed_axes` instead of `permuted_axes` --- pineappl/src/evolution.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index f7ddc23e..99deaf5b 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -359,7 +359,7 @@ fn operator_slices( operator .slice(s![pid1_idx, .., pid0_idx, ..]) .select(Axis(0), &x1_indices) - .permuted_axes([1, 0]) + .reversed_axes() .as_standard_layout() .into_owned() }) From aa8bb9131320b296725e4f8fb0b676eacf54d670 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 18 Oct 2023 14:40:05 +0200 Subject: [PATCH 21/40] Use `general_mat_mul` instead of `scaled_add`/`dot`/`t` --- pineappl/src/evolution.rs | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 99deaf5b..91bcaef2 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -8,6 +8,7 @@ use super::sparse_array3::SparseArray3; use super::subgrid::{Mu2, Subgrid, SubgridEnum}; use float_cmp::approx_eq; use itertools::Itertools; +use ndarray::linalg; use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView4, ArrayView5, Axis}; use std::iter; @@ -1063,6 +1064,8 @@ pub(crate) fn evolve_slice_with_two( last_x1b = x1_b; } + let mut tmp = Array2::zeros((last_x1a.len(), info.x0.len())); + for &(pida1, pidb1, factor) in lumi1.entry() { for (fk_table, opa, opb) in lumi0 @@ -1083,7 +1086,8 @@ pub(crate) fn evolve_slice_with_two( .map(|(opa, opb)| (fk_table, opa, opb)) }) { - fk_table.scaled_add(factor, &opa.dot(&array.dot(&opb.t()))); + linalg::general_mat_mul(1.0, &array, &opb.t(), 0.0, &mut tmp); + linalg::general_mat_mul(factor, opa, &tmp, 1.0, fk_table); } } } From 890c2fb85ae7b7198736af74c5945f5763a5666d Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 21 Oct 2023 13:05:27 +0200 Subject: [PATCH 22/40] Replace `evolve_slice` with `evolve_with_slice_iter` --- pineappl/Cargo.toml | 1 + pineappl/src/grid.rs | 119 +++++++++++++++++++++++++------------ pineappl_cli/src/evolve.rs | 29 +-------- 3 files changed, 84 insertions(+), 65 deletions(-) diff --git a/pineappl/Cargo.toml b/pineappl/Cargo.toml index 280fe4c4..8c6df8a8 100644 --- a/pineappl/Cargo.toml +++ b/pineappl/Cargo.toml @@ -13,6 +13,7 @@ rust-version.workspace = true version.workspace = true [dependencies] +anyhow = "1.0.48" arrayvec = "0.7.2" bincode = "1.3.3" bitflags = "1.3.2" diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 757c3bea..ee283a2a 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -18,7 +18,7 @@ use git_version::git_version; use indicatif::{ProgressBar, ProgressStyle}; use itertools::Itertools; use lz4_flex::frame::{FrameDecoder, FrameEncoder}; -use ndarray::{s, Array3, Array5, ArrayView4, ArrayView5, Axis, Dimension}; +use ndarray::{s, Array3, Array4, Array5, ArrayView5, Axis, Dimension}; use serde::{Deserialize, Serialize}; use std::borrow::Cow; use std::cmp::Ordering; @@ -280,6 +280,9 @@ pub enum GridError { /// Returned from [`Grid::evolve`] if the evolution failed. #[error("failed to evolve grid: {0}")] EvolutionFailure(String), + /// Errors that do no originate from this crate itself. + #[error(transparent)] + Other(#[from] anyhow::Error), } #[derive(Clone, Deserialize, Serialize)] @@ -1945,55 +1948,97 @@ impl Grid { Ok(FkTable::try_from(grid).unwrap_or_else(|_| unreachable!())) } - /// Converts this `Grid` into an [`FkTable`] using an evolution kernel operator (EKO) slice - /// given as `operator`. The dimensions and properties of this operator must be described using - /// `info`. 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. + // TODO: + // - try to find a better solution than to require that E must be convertible into + // anyhow::Error + // - change the iterator to a `CowArray`? + + /// 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`. - pub fn evolve_slice( + /// incompatible with this `Grid`. Returns a [`GridError::Other`] if the iterator from `slices` + /// return an error. + pub fn evolve_with_slice_iter>( &self, - operator: ArrayView4, - info: &OperatorSliceInfo, + slices: impl IntoIterator), E>>, order_mask: &[bool], ) -> Result { - let op_info_dim = ( - info.pids1.len(), - info.x1.len(), - info.pids0.len(), - info.x0.len(), - ); + use super::evolution::EVOLVE_INFO_TOL_ULPS; - if operator.dim() != op_info_dim { - return Err(GridError::EvolutionFailure(format!( - "operator information {:?} does not match the operator's dimensions: {:?}", - op_info_dim, - operator.dim(), - ))); + let mut lhs: Option = None; + let mut fac1 = Vec::new(); + + for result in slices { + let (info, operator) = result.map_err(|err| GridError::Other(err.into()))?; + + let op_info_dim = ( + info.pids1.len(), + info.x1.len(), + info.pids0.len(), + info.x0.len(), + ); + + if operator.dim() != op_info_dim { + return Err(GridError::EvolutionFailure(format!( + "operator information {:?} does not match the operator's dimensions: {:?}", + op_info_dim, + operator.dim(), + ))); + } + + let view = operator.view(); + + let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { + evolution::evolve_slice_with_two(self, &view, &info, order_mask) + } else { + evolution::evolve_slice_with_one(self, &view, &info, order_mask) + }?; + + let mut rhs = Self { + subgrids, + lumi, + bin_limits: self.bin_limits.clone(), + orders: vec![Order::new(0, 0, 0, 0)], + subgrid_params: SubgridParams::default(), + more_members: self.more_members.clone(), + }; + + // write additional metadata + rhs.set_key_value("lumi_id_types", &info.lumi_id_types); + + if let Some(lhs) = &mut lhs { + lhs.merge(rhs)?; + } else { + lhs = Some(rhs); + } + + fac1.push(info.fac1); } - let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { - evolution::evolve_slice_with_two(self, &operator, info, order_mask) - } else { - evolution::evolve_slice_with_one(self, &operator, info, order_mask) - }?; + // UNWRAP: if we can't compare two numbers there's a bug + fac1.sort_by(|a, b| a.partial_cmp(b).unwrap()); - let mut grid = Self { - subgrids, - lumi, - bin_limits: self.bin_limits.clone(), - orders: vec![Order::new(0, 0, 0, 0)], - subgrid_params: SubgridParams::default(), - more_members: self.more_members.clone(), - }; + // make sure we've evolved all slices + if let Some(muf2) = self.evolve_info(&order_mask).fac1.into_iter().find(|&x| { + !fac1 + .iter() + .any(|&y| approx_eq!(f64, x, y, ulps = EVOLVE_INFO_TOL_ULPS)) + }) { + return Err(GridError::EvolutionFailure(format!( + "no operator for muf2 = {muf2} found in {fac1:#?}" + ))); + } - // write additional metadata - grid.set_key_value("lumi_id_types", &info.lumi_id_types); + // UNWRAP: should panick when all subgrids are empty + let grid = lhs.unwrap(); + // UNWRAP: merging evolved slices should be a proper FkTable again Ok(FkTable::try_from(grid).unwrap_or_else(|_| unreachable!())) } diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index a993f93d..9cbbd0e6 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -479,7 +479,6 @@ fn evolve_grid( xif: f64, ) -> Result { use eko::EkoSlices; - use float_cmp::approx_eq; use pineappl::subgrid::{Mu2, Subgrid}; let mut ren1: Vec<_> = grid @@ -510,34 +509,8 @@ fn evolve_grid( .collect(); let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas, xir, xif)?; - let grid_muf2 = grid.evolve_info(&order_mask).fac1; - // check that every factorization-scale value in the grid has an operator - assert!(!grid_muf2 - .iter() - .any(|&muf2| !eko_slices.fac1().iter().any(|&op_muf2| approx_eq!( - f64, - muf2, - op_muf2, - ulps = 64 - )))); - - let mut lhs: Option = None; - - for item in eko_slices.iter_mut() { - let (info, operator) = item?; - let rhs = grid - .evolve_slice(operator.view(), &info, &order_mask)? - .into_grid(); - - if let Some(lhs) = &mut lhs { - lhs.merge(rhs)?; - } else { - lhs = Some(rhs); - } - } - - Ok(FkTable::try_from(lhs.unwrap()).unwrap()) + Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask)?) } #[cfg(not(feature = "evolve"))] From 082a764c8889c5fecf3ad47356c67849aebc0c8c Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 24 Oct 2023 15:49:43 +0200 Subject: [PATCH 23/40] Add comment to check for possible bug --- pineappl/src/grid.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index ee283a2a..1c7343e2 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -2024,6 +2024,8 @@ impl Grid { // UNWRAP: if we can't compare two numbers there's a bug fac1.sort_by(|a, b| a.partial_cmp(b).unwrap()); + // TODO: here there's possibly is a bug if xif isn't one + // make sure we've evolved all slices if let Some(muf2) = self.evolve_info(&order_mask).fac1.into_iter().find(|&x| { !fac1 From 36147118fc1822b9931d3d101bcd8f196c52f140 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 24 Oct 2023 15:50:11 +0200 Subject: [PATCH 24/40] Remove old functions and migrate `Grid::evolve` to use new method --- pineappl/src/evolution.rs | 505 +------------------------------------- pineappl/src/grid.rs | 59 ++--- 2 files changed, 25 insertions(+), 539 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 91bcaef2..da2d456d 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -9,7 +9,7 @@ use super::subgrid::{Mu2, Subgrid, SubgridEnum}; use float_cmp::approx_eq; use itertools::Itertools; use ndarray::linalg; -use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView4, ArrayView5, Axis}; +use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView4, Axis}; use std::iter; /// Number of ULPS used to de-duplicate grid values in [`Grid::evolve_info`]. @@ -162,54 +162,6 @@ fn gluon_has_pid_zero(grid: &Grid) -> bool { type Pid01IndexTuples = Vec<(usize, usize)>; type Pid01Tuples = Vec<(i32, i32)>; -fn pids( - operator: &ArrayView5, - info: &OperatorInfo, - gluon_has_pid_zero: bool, - pid1_nonzero: &dyn Fn(i32) -> bool, -) -> Result<(Pid01IndexTuples, Pid01Tuples), GridError> { - // list of all non-zero PID indices - let pid_indices: Vec<_> = (0..operator.dim().3) - .cartesian_product(0..operator.dim().1) - .filter(|&(pid0_idx, pid1_idx)| { - // 1) at least one element of the operator must be non-zero, and 2) the pid must be - // contained in the lumi somewhere - operator - .slice(s![.., pid1_idx, .., pid0_idx, ..]) - .iter() - .any(|&value| value != 0.0) - && pid1_nonzero(if gluon_has_pid_zero && info.pids1[pid1_idx] == 21 { - 0 - } else { - info.pids1[pid1_idx] - }) - }) - .collect(); - - if pid_indices.is_empty() { - return Err(GridError::EvolutionFailure( - "no non-zero operator found; result would be an empty FkTable".to_string(), - )); - } - - // list of all non-zero (pid0, pid1) combinations - let pids = pid_indices - .iter() - .map(|&(pid0_idx, pid1_idx)| { - ( - info.pids0[pid0_idx], - if gluon_has_pid_zero && info.pids1[pid1_idx] == 21 { - 0 - } else { - info.pids1[pid1_idx] - }, - ) - }) - .collect(); - - Ok((pid_indices, pids)) -} - fn pid_slices( operator: &ArrayView4, info: &OperatorSliceInfo, @@ -281,58 +233,6 @@ fn lumi0_with_two(pids_a: &[(i32, i32)], pids_b: &[(i32, i32)]) -> Vec<(i32, i32 .collect() } -fn operators( - operator: &ArrayView5, - info: &OperatorInfo, - fac1: &[f64], - pid_indices: &[(usize, usize)], - x1: &[f64], -) -> Result>, GridError> { - // permutation between the grid fac1 values and the operator fac1 values - let fac1_indices: Vec<_> = fac1 - .iter() - .map(|&fac1p| { - info.fac1 - .iter() - .position(|&fac1| approx_eq!(f64, fac1p, fac1, ulps = EVOLUTION_TOL_ULPS)) - .ok_or_else(|| { - GridError::EvolutionFailure(format!("no operator for muf2 = {fac1p} found")) - }) - }) - // TODO: use `try_collect` once stabilized - .collect::>()?; - - // permutation between the grid x values and the operator x1 values - let x1_indices: Vec<_> = x1 - .iter() - .map(|&x1p| { - info.x1 - .iter() - .position(|&x1| approx_eq!(f64, x1p, x1, ulps = EVOLUTION_TOL_ULPS)) - .ok_or_else(|| { - GridError::EvolutionFailure(format!("no operator for x = {x1p} found")) - }) - }) - // TODO: use `try_collect` once stabilized - .collect::>()?; - - // create the corresponding operators accessible in the form [muf2, x0, x1] - let operators: Vec<_> = pid_indices - .iter() - .map(|&(pid0_idx, pid1_idx)| { - operator - .slice(s![.., pid1_idx, .., pid0_idx, ..]) - .select(Axis(0), &fac1_indices) - .select(Axis(1), &x1_indices) - .permuted_axes([0, 2, 1]) - .as_standard_layout() - .into_owned() - }) - .collect(); - - Ok(operators) -} - fn operator_slices( operator: &ArrayView4, info: &OperatorSliceInfo, @@ -369,147 +269,6 @@ fn operator_slices( Ok(operators) } -type Fac1X1aX1bOp3Tuple = (Vec, Vec, Vec, Array3); - -fn ndarray_from_subgrid_orders( - info: &OperatorInfo, - subgrids: &ArrayView1, - orders: &[Order], - order_mask: &[bool], -) -> Result { - // TODO: skip empty subgrids - - let mut fac1: Vec<_> = subgrids - .iter() - .enumerate() - .filter(|(index, _)| order_mask.get(*index).copied().unwrap_or(true)) - .flat_map(|(_, subgrid)| { - subgrid - .mu2_grid() - .iter() - .map(|mu2| info.xif * info.xif * mu2.fac) - .collect::>() - }) - .collect(); - let mut x1_a: Vec<_> = subgrids - .iter() - .enumerate() - .filter(|(index, _)| order_mask.get(*index).copied().unwrap_or(true)) - .flat_map(|(_, subgrid)| subgrid.x1_grid().into_owned()) - .collect(); - let mut x1_b: Vec<_> = subgrids - .iter() - .enumerate() - .filter(|(index, _)| order_mask.get(*index).copied().unwrap_or(true)) - .flat_map(|(_, subgrid)| subgrid.x2_grid().into_owned()) - .collect(); - - fac1.sort_by(f64::total_cmp); - fac1.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLUTION_TOL_ULPS)); - x1_a.sort_by(f64::total_cmp); - x1_a.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLUTION_TOL_ULPS)); - x1_b.sort_by(f64::total_cmp); - x1_b.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLUTION_TOL_ULPS)); - - let mut array = Array3::::zeros((fac1.len(), x1_a.len(), x1_b.len())); - - // add subgrids for different orders, but the same bin and lumi, using the right - // couplings - for (subgrid, order) in subgrids - .iter() - .zip(orders.iter()) - .zip(order_mask.iter().chain(iter::repeat(&true))) - .filter_map(|((subgrid, order), &enabled)| { - (enabled && !subgrid.is_empty()).then_some((subgrid, order)) - }) - { - let mut logs = 1.0; - - if order.logxir > 0 { - if approx_eq!(f64, info.xir, 1.0, ulps = 4) { - continue; - } - - logs *= (info.xir * info.xir).ln(); - } - - if order.logxif > 0 { - if approx_eq!(f64, info.xif, 1.0, ulps = 4) { - continue; - } - - logs *= (info.xif * info.xif).ln(); - } - - // TODO: use `try_collect` once stabilized - let fac1_indices: Vec<_> = subgrid - .mu2_grid() - .iter() - .map(|&Mu2 { fac, .. }| { - fac1.iter() - .position(|&scale| { - approx_eq!( - f64, - info.xif * info.xif * fac, - scale, - ulps = EVOLUTION_TOL_ULPS - ) - }) - .ok_or_else(|| { - GridError::EvolutionFailure(format!("no operator for muf2 = {fac} found")) - }) - }) - .collect::>()?; - let xa_indices: Vec<_> = subgrid - .x1_grid() - .iter() - .map(|&xa| { - x1_a.iter() - .position(|&x1a| approx_eq!(f64, x1a, xa, ulps = EVOLUTION_TOL_ULPS)) - .ok_or_else(|| { - GridError::EvolutionFailure(format!("no operator for x1 = {xa} found")) - }) - }) - .collect::>()?; - let xb_indices: Vec<_> = subgrid - .x2_grid() - .iter() - .map(|&xb| { - x1_b.iter() - .position(|&x1b| approx_eq!(f64, x1b, xb, ulps = EVOLUTION_TOL_ULPS)) - .ok_or_else(|| { - GridError::EvolutionFailure(format!("no operator for x1 = {xb} found")) - }) - }) - .collect::>()?; - - for ((ifac1, ix1, ix2), value) in subgrid.indexed_iter() { - let mur2 = info.xir * info.xir * subgrid.mu2_grid()[ifac1].ren; - - let als = if order.alphas == 0 { - 1.0 - } else if let Some(alphas) = - info.ren1 - .iter() - .zip(info.alphas.iter()) - .find_map(|(&ren1, &alphas)| { - approx_eq!(f64, ren1, mur2, ulps = EVOLUTION_TOL_ULPS).then(|| alphas) - }) - { - alphas.powi(order.alphas.try_into().unwrap()) - } else { - return Err(GridError::EvolutionFailure(format!( - "no alphas for mur2 = {mur2} found" - ))); - }; - - array[[fac1_indices[ifac1], xa_indices[ix1], xb_indices[ix2]]] += als * logs * value; - } - } - - Ok((fac1, x1_a, x1_b, array)) -} - type X1aX1bOp2Tuple = (Vec, Vec, Array2); fn ndarray_from_subgrid_orders_slice( @@ -631,140 +390,6 @@ fn ndarray_from_subgrid_orders_slice( Ok((x1_a, x1_b, array)) } -pub(crate) fn evolve_with_one( - grid: &Grid, - operator: &ArrayView5, - info: &OperatorInfo, - order_mask: &[bool], -) -> Result<(Array3, Vec), GridError> { - let gluon_has_pid_zero = gluon_has_pid_zero(grid); - let has_pdf1 = grid.has_pdf1(); - - let (pid_indices, pids) = pids(operator, info, gluon_has_pid_zero, &|pid| { - grid.lumi() - .iter() - .flat_map(LumiEntry::entry) - .any(|&(a, b, _)| if has_pdf1 { a } else { b } == pid) - })?; - - let lumi0 = lumi0_with_one(&pids); - let mut sub_fk_tables = Vec::with_capacity(grid.bin_info().bins() * lumi0.len()); - let new_axis = if has_pdf1 { 2 } else { 1 }; - - let mut last_x1 = Vec::new(); - let mut last_fac1 = Vec::new(); - let mut ops = Vec::new(); - - for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { - let mut tables = vec![Array1::zeros(info.x0.len()); lumi0.len()]; - - for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { - let (fac1, x1_a, x1_b, array) = - ndarray_from_subgrid_orders(info, &subgrids_o, grid.orders(), order_mask)?; - - let x1 = if has_pdf1 { x1_a } else { x1_b }; - - if x1.is_empty() { - continue; - } - - if (last_fac1.len() != fac1.len()) - || last_fac1 - .iter() - .zip(fac1.iter()) - .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) - || (last_x1.len() != x1.len()) - || last_x1 - .iter() - .zip(x1.iter()) - .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) - { - ops = operators(operator, info, &fac1, &pid_indices, &x1)?; - last_fac1 = fac1; - last_x1 = x1; - } - - for (&pid1, &factor) in lumi1 - .entry() - .iter() - .map(|(a, b, f)| if has_pdf1 { (a, f) } else { (b, f) }) - { - for (fk_table, op) in - lumi0 - .iter() - .zip(tables.iter_mut()) - .filter_map(|(&pid0, fk_table)| { - pids.iter() - .zip(ops.iter()) - .find_map(|(&(p0, p1), op)| { - (p0 == pid0 && p1 == pid1).then_some(op) - }) - .map(|op| (fk_table, op)) - }) - { - let mut result = Array1::zeros(info.x0.len()); - - for imu2 in 0..array.dim().0 { - let op = op.index_axis(Axis(0), imu2); - - result += &op.dot( - &array - .index_axis(Axis(0), imu2) - .index_axis(Axis(new_axis - 1), 0), - ); - } - - fk_table.scaled_add(factor, &result); - } - } - } - - sub_fk_tables.extend(tables.into_iter().map(|table| { - ImportOnlySubgridV2::new( - SparseArray3::from_ndarray( - table - .insert_axis(Axis(0)) - .insert_axis(Axis(new_axis)) - .view(), - 0, - 1, - ), - vec![Mu2 { - // TODO: FK tables don't depend on the renormalization scale - //ren: -1.0, - ren: info.fac0, - fac: info.fac0, - }], - if has_pdf1 { info.x0.clone() } else { vec![1.0] }, - if has_pdf1 { vec![1.0] } else { info.x0.clone() }, - ) - .into() - })); - } - - let pid = if has_pdf1 { - grid.initial_state_2() - } else { - grid.initial_state_1() - }; - - Ok(( - Array1::from_iter(sub_fk_tables) - .into_shape((1, grid.bin_info().bins(), lumi0.len())) - .unwrap(), - lumi0 - .iter() - .map(|&a| { - lumi_entry![ - if has_pdf1 { a } else { pid }, - if has_pdf1 { pid } else { a }, - 1.0 - ] - }) - .collect(), - )) -} - pub(crate) fn evolve_slice_with_one( grid: &Grid, operator: &ArrayView4, @@ -880,134 +505,6 @@ pub(crate) fn evolve_slice_with_one( )) } -pub(crate) fn evolve_with_two( - grid: &Grid, - operator: &ArrayView5, - info: &OperatorInfo, - order_mask: &[bool], -) -> Result<(Array3, Vec), GridError> { - let gluon_has_pid_zero = gluon_has_pid_zero(grid); - - let (pid_indices_a, pids_a) = pids(operator, info, gluon_has_pid_zero, &|pid1| { - grid.lumi() - .iter() - .flat_map(LumiEntry::entry) - .any(|&(a, _, _)| a == pid1) - })?; - let (pid_indices_b, pids_b) = pids(operator, info, gluon_has_pid_zero, &|pid1| { - grid.lumi() - .iter() - .flat_map(LumiEntry::entry) - .any(|&(_, b, _)| b == pid1) - })?; - - let lumi0 = lumi0_with_two(&pids_a, &pids_b); - let mut sub_fk_tables = Vec::with_capacity(grid.bin_info().bins() * lumi0.len()); - - let mut last_fac1 = Vec::new(); - let mut last_x1a = Vec::new(); - let mut last_x1b = Vec::new(); - let mut operators_a = Vec::new(); - let mut operators_b = Vec::new(); - - for subgrids_ol in grid.subgrids().axis_iter(Axis(1)) { - let mut tables = vec![Array2::zeros((info.x0.len(), info.x0.len())); lumi0.len()]; - - for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { - let (fac1, x1_a, x1_b, array) = - ndarray_from_subgrid_orders(info, &subgrids_o, grid.orders(), order_mask)?; - - let fac1_diff = (last_fac1.len() != fac1.len()) - || last_fac1 - .iter() - .zip(fac1.iter()) - .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)); - - if fac1_diff - || (last_x1a.len() != x1_a.len()) - || last_x1a - .iter() - .zip(x1_a.iter()) - .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) - { - operators_a = operators(operator, info, &fac1, &pid_indices_a, &x1_a)?; - last_x1a = x1_a; - } - - if fac1_diff - || (last_x1b.len() != x1_b.len()) - || last_x1b - .iter() - .zip(x1_b.iter()) - .any(|(&lhs, &rhs)| !approx_eq!(f64, lhs, rhs, ulps = EVOLUTION_TOL_ULPS)) - { - operators_b = operators(operator, info, &fac1, &pid_indices_b, &x1_b)?; - last_x1b = x1_b; - } - - if fac1_diff { - last_fac1 = fac1; - }; - - for &(pida1, pidb1, factor) in lumi1.entry() { - for (fk_table, opa, opb) in - lumi0 - .iter() - .zip(tables.iter_mut()) - .filter_map(|(&(pida0, pidb0), fk_table)| { - pids_a - .iter() - .zip(operators_a.iter()) - .find_map(|(&(pa0, pa1), opa)| { - (pa0 == pida0 && pa1 == pida1).then_some(opa) - }) - .zip(pids_b.iter().zip(operators_b.iter()).find_map( - |(&(pb0, pb1), opb)| { - (pb0 == pidb0 && pb1 == pidb1).then_some(opb) - }, - )) - .map(|(opa, opb)| (fk_table, opa, opb)) - }) - { - let mut result = Array2::zeros((info.x0.len(), info.x0.len())); - - for imu2 in 0..array.dim().0 { - let opa = opa.index_axis(Axis(0), imu2); - let opb = opb.index_axis(Axis(0), imu2); - let arr = array.index_axis(Axis(0), imu2); - - result += &opa.dot(&arr.dot(&opb.t())); - } - - fk_table.scaled_add(factor, &result); - } - } - } - - 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: info.fac0, - fac: info.fac0, - }], - info.x0.clone(), - info.x0.clone(), - ) - .into() - })); - } - - Ok(( - Array1::from_iter(sub_fk_tables) - .into_shape((1, grid.bin_info().bins(), lumi0.len())) - .unwrap(), - lumi0.iter().map(|&(a, b)| lumi_entry![a, b, 1.0]).collect(), - )) -} - pub(crate) fn evolve_slice_with_two( grid: &Grid, operator: &ArrayView4, diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 1c7343e2..c5e41134 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -1911,41 +1911,30 @@ impl Grid { info: &OperatorInfo, order_mask: &[bool], ) -> Result { - let op_info_dim = ( - info.fac1.len(), - info.pids1.len(), - info.x1.len(), - info.pids0.len(), - info.x0.len(), - ); - - if operator.dim() != op_info_dim { - return Err(GridError::EvolutionFailure(format!( - "operator information {:?} does not match the operator's dimensions: {:?}", - op_info_dim, - operator.dim(), - ))); - } - - let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { - evolution::evolve_with_two(self, &operator, info, order_mask) - } else { - evolution::evolve_with_one(self, &operator, info, order_mask) - }?; - - let mut grid = Self { - subgrids, - lumi, - bin_limits: self.bin_limits.clone(), - orders: vec![Order::new(0, 0, 0, 0)], - subgrid_params: SubgridParams::default(), - more_members: self.more_members.clone(), - }; - - // write additional metadata - grid.set_key_value("lumi_id_types", &info.lumi_id_types); - - Ok(FkTable::try_from(grid).unwrap_or_else(|_| unreachable!())) + self.evolve_with_slice_iter( + info.fac1 + .iter() + .zip(operator.axis_iter(Axis(0))) + .map(|(&fac1, op)| { + Ok::<_, GridError>(( + OperatorSliceInfo { + fac0: info.fac0.clone(), + pids0: info.pids0.clone(), + x0: info.x0.clone(), + fac1, + pids1: info.pids1.clone(), + x1: info.x1.clone(), + ren1: info.ren1.clone(), + alphas: info.alphas.clone(), + xir: info.xir, + xif: info.xif, + lumi_id_types: info.lumi_id_types.clone(), + }, + op.into_owned(), + )) + }), + order_mask, + ) } // TODO: From 36dce35588e97cbe839f5e49756a333157bd9b8a Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 12 Dec 2023 13:45:38 +0100 Subject: [PATCH 25/40] Fix PID = 0 gluon detection --- pineappl/src/evolution.rs | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index da2d456d..cf964dca 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -147,16 +147,16 @@ pub struct OperatorSliceInfo { } fn gluon_has_pid_zero(grid: &Grid) -> bool { - grid.key_values() - .and_then(|key_values| key_values.get("lumi_id_types")) - .and_then(|value| { - (value == "pdg_mc_ids").then(|| { - grid.lumi() - .iter() - .any(|entry| entry.entry().iter().any(|&(a, b, _)| (a == 0) || (b == 0))) - }) - }) - .unwrap_or(false) + // if there are any PID zero particles ... + grid.lumi() + .iter() + .any(|entry| entry.entry().iter().any(|&(a, b, _)| (a == 0) || (b == 0))) + // and if lumi_id_types = pdg_mc_ids or if the key-value pair doesn't exist + && grid + .key_values() + .and_then(|key_values| key_values.get("lumi_id_types")) + .map(|value| value == "pdg_mc_ids") + .unwrap_or(true) } type Pid01IndexTuples = Vec<(usize, usize)>; From ded265b632698708933b8aebb3be42b7fa1d59fa Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 1 Feb 2024 10:00:05 +0100 Subject: [PATCH 26/40] Use `CowArray` instead of `Array` to possibly avoid copying large arrays --- pineappl/src/grid.rs | 9 ++++----- pineappl_cli/src/evolve.rs | 11 +++++------ 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index c5e41134..dfb61d9b 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -18,7 +18,7 @@ use git_version::git_version; use indicatif::{ProgressBar, ProgressStyle}; use itertools::Itertools; use lz4_flex::frame::{FrameDecoder, FrameEncoder}; -use ndarray::{s, Array3, Array4, Array5, ArrayView5, Axis, Dimension}; +use ndarray::{s, Array3, Array5, ArrayView5, Axis, CowArray, Dimension, Ix4}; use serde::{Deserialize, Serialize}; use std::borrow::Cow; use std::cmp::Ordering; @@ -1930,7 +1930,7 @@ impl Grid { xif: info.xif, lumi_id_types: info.lumi_id_types.clone(), }, - op.into_owned(), + CowArray::from(op), )) }), order_mask, @@ -1940,7 +1940,6 @@ impl Grid { // TODO: // - try to find a better solution than to require that E must be convertible into // anyhow::Error - // - change the iterator to a `CowArray`? /// 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 @@ -1953,9 +1952,9 @@ impl Grid { /// 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_iter>( + pub fn evolve_with_slice_iter<'a, E: Into>( &self, - slices: impl IntoIterator), E>>, + slices: impl IntoIterator), E>>, order_mask: &[bool], ) -> Result { use super::evolution::EVOLVE_INFO_TOL_ULPS; diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 9cbbd0e6..67bd3036 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -18,7 +18,7 @@ mod eko { use either::Either; use lz4_flex::frame::FrameDecoder; use ndarray::iter::AxisIter; - use ndarray::{Array4, Array5, Axis, Ix4}; + use ndarray::{Array4, Array5, Axis, CowArray, Ix4}; use ndarray_npy::{NpzReader, ReadNpyExt}; use pineappl::evolution::OperatorSliceInfo; use pineappl::pids; @@ -386,7 +386,7 @@ mod eko { } impl<'a> IntoIterator for &'a mut EkoSlices { - type Item = Result<(OperatorSliceInfo, Array4)>; + type Item = Result<(OperatorSliceInfo, CowArray<'a, f64, Ix4>)>; type IntoIter = EkoSlicesIter<'a>; fn into_iter(self) -> Self::IntoIter { @@ -407,7 +407,7 @@ mod eko { } impl<'a> Iterator for EkoSlicesIter<'a> { - type Item = Result<(OperatorSliceInfo, Array4)>; + type Item = Result<(OperatorSliceInfo, CowArray<'a, f64, Ix4>)>; fn next(&mut self) -> Option { match self { @@ -416,8 +416,7 @@ mod eko { let mut info = info.clone(); info.fac1 = *fac1; - // TODO: see if we can replace this some kind of Cow structure - Some(Ok((info, operator.to_owned()))) + Some(Ok((info, CowArray::from(operator)))) } else { None } @@ -455,7 +454,7 @@ mod eko { let mut info = info.clone(); info.fac1 = fac1.get(&file_stem).copied().ok_or_else(|| anyhow!("file '{}.yaml' not found, could not determine the operator's factorization scale", file_stem.to_string_lossy()))?; - return Ok(Some((info, operator))); + return Ok(Some((info, CowArray::from(operator)))); } } From d2a4f3d72e1edf898dbb653db848fa7ab9d10d33 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 13:18:21 +0100 Subject: [PATCH 27/40] Make `xif` parameter of `Grid::evolve_with_slice_iter` --- pineappl/src/evolution.rs | 36 +++++++++++++++++++++--------------- pineappl/src/grid.rs | 7 ++++--- pineappl_cli/src/evolve.rs | 24 ++++++------------------ 3 files changed, 31 insertions(+), 36 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index cf964dca..1c6c2481 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -140,8 +140,6 @@ pub struct OperatorSliceInfo { pub alphas: Vec, /// Multiplicative factor for the central renormalization scale. pub xir: f64, - /// Multiplicative factor for the central factorization scale. - pub xif: f64, /// Identifier of the particle basis for the `FkTable`. pub lumi_id_types: String, } @@ -276,10 +274,11 @@ fn ndarray_from_subgrid_orders_slice( subgrids: &ArrayView1, orders: &[Order], order_mask: &[bool], + xif: f64, ) -> Result { // TODO: skip empty subgrids - let fac1 = info.xif * info.xif * info.fac1; + let fac1 = xif * xif * info.fac1; let mut x1_a: Vec<_> = subgrids .iter() .enumerate() @@ -321,11 +320,11 @@ fn ndarray_from_subgrid_orders_slice( } if order.logxif > 0 { - if approx_eq!(f64, info.xif, 1.0, ulps = 4) { + if approx_eq!(f64, xif, 1.0, ulps = 4) { continue; } - logs *= (info.xif * info.xif).ln(); + logs *= (xif * xif).ln(); } // TODO: use `try_collect` once stabilized @@ -355,12 +354,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, - info.xif * info.xif * fac, - fac1, - ulps = EVOLUTION_TOL_ULPS - ) { + if !approx_eq!(f64, xif * xif * fac, fac1, ulps = EVOLUTION_TOL_ULPS) { continue; } @@ -395,6 +389,7 @@ pub(crate) fn evolve_slice_with_one( operator: &ArrayView4, info: &OperatorSliceInfo, order_mask: &[bool], + xif: f64, ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); let has_pdf1 = grid.has_pdf1(); @@ -417,8 +412,13 @@ pub(crate) fn evolve_slice_with_one( let mut tables = vec![Array1::zeros(info.x0.len()); lumi0.len()]; for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { - let (x1_a, x1_b, array) = - ndarray_from_subgrid_orders_slice(info, &subgrids_o, grid.orders(), order_mask)?; + let (x1_a, x1_b, array) = ndarray_from_subgrid_orders_slice( + info, + &subgrids_o, + grid.orders(), + order_mask, + xif, + )?; let x1 = if has_pdf1 { x1_a } else { x1_b }; @@ -510,6 +510,7 @@ pub(crate) fn evolve_slice_with_two( operator: &ArrayView4, info: &OperatorSliceInfo, order_mask: &[bool], + xif: f64, ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); @@ -538,8 +539,13 @@ pub(crate) fn evolve_slice_with_two( let mut tables = vec![Array2::zeros((info.x0.len(), info.x0.len())); lumi0.len()]; for (subgrids_o, lumi1) in subgrids_ol.axis_iter(Axis(1)).zip(grid.lumi()) { - let (x1_a, x1_b, array) = - ndarray_from_subgrid_orders_slice(info, &subgrids_o, grid.orders(), order_mask)?; + let (x1_a, x1_b, array) = ndarray_from_subgrid_orders_slice( + info, + &subgrids_o, + grid.orders(), + order_mask, + xif, + )?; if (last_x1a.len() != x1_a.len()) || last_x1a diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index dfb61d9b..82c8d391 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -1927,13 +1927,13 @@ impl Grid { ren1: info.ren1.clone(), alphas: info.alphas.clone(), xir: info.xir, - xif: info.xif, lumi_id_types: info.lumi_id_types.clone(), }, CowArray::from(op), )) }), order_mask, + info.xif, ) } @@ -1956,6 +1956,7 @@ impl Grid { &self, slices: impl IntoIterator), E>>, order_mask: &[bool], + xif: f64, ) -> Result { use super::evolution::EVOLVE_INFO_TOL_ULPS; @@ -1983,9 +1984,9 @@ impl Grid { let view = operator.view(); let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { - evolution::evolve_slice_with_two(self, &view, &info, order_mask) + evolution::evolve_slice_with_two(self, &view, &info, order_mask, xif) } else { - evolution::evolve_slice_with_one(self, &view, &info, order_mask) + evolution::evolve_slice_with_one(self, &view, &info, order_mask, xif) }?; let mut rhs = Self { diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 67bd3036..cb3f5d97 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -127,19 +127,13 @@ mod eko { Err(anyhow!("no file 'metadata.yaml' in EKO archive found")) } - pub fn new( - eko_path: &Path, - ren1: &[f64], - alphas: &[f64], - xir: f64, - xif: f64, - ) -> Result { + pub fn new(eko_path: &Path, ren1: &[f64], alphas: &[f64], xir: f64) -> Result { let metadata = Self::read_metadata(eko_path)?; match metadata { - Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas, xir, xif), - Metadata::V1(v1) => Self::with_v1(v1, eko_path, ren1, alphas, xir, xif), - Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas, xir, xif), + Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas, xir), + Metadata::V1(v1) => Self::with_v1(v1, eko_path, ren1, alphas, xir), + Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas, xir), } } @@ -149,7 +143,6 @@ mod eko { ren1: &[f64], alphas: &[f64], xir: f64, - xif: f64, ) -> Result { let mut operator = None; @@ -178,7 +171,6 @@ mod eko { ren1: ren1.to_vec(), alphas: alphas.to_vec(), xir, - xif, }, operator, }) @@ -190,7 +182,6 @@ mod eko { ren1: &[f64], alphas: &[f64], xir: f64, - xif: f64, ) -> Result { let mut fac1 = HashMap::new(); let base64 = GeneralPurpose::new(&URL_SAFE, PAD); @@ -265,7 +256,6 @@ mod eko { ren1: ren1.to_vec(), alphas: alphas.to_vec(), xir, - xif, }, archive: Archive::new(File::open(eko_path)?), }) @@ -277,7 +267,6 @@ mod eko { ren1: &[f64], alphas: &[f64], xir: f64, - xif: f64, ) -> Result { let mut fac1 = HashMap::new(); let mut operator: Option = None; @@ -344,7 +333,6 @@ mod eko { ren1: ren1.to_vec(), alphas: alphas.to_vec(), xir, - xif, }, archive: Archive::new(File::open(eko_path)?), }) @@ -507,9 +495,9 @@ fn evolve_grid( }) .collect(); - let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas, xir, xif)?; + let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas, xir)?; - Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask)?) + Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, xif)?) } #[cfg(not(feature = "evolve"))] From b3450a5c27e65eeb088daff20e78b143bd735e0a Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 13:36:11 +0100 Subject: [PATCH 28/40] Make `xir` parameter of `Grid::evolve_with_slice_iter` --- pineappl/src/evolution.rs | 15 ++++++++------- pineappl/src/grid.rs | 9 ++++----- pineappl_cli/src/evolve.rs | 18 ++++++------------ 3 files changed, 18 insertions(+), 24 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 1c6c2481..cf2ee2fe 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -138,8 +138,6 @@ pub struct OperatorSliceInfo { pub ren1: Vec, /// Strong couplings corresponding to the order given in [`ren1`](Self::ren1). pub alphas: Vec, - /// Multiplicative factor for the central renormalization scale. - pub xir: f64, /// Identifier of the particle basis for the `FkTable`. pub lumi_id_types: String, } @@ -274,6 +272,7 @@ fn ndarray_from_subgrid_orders_slice( subgrids: &ArrayView1, orders: &[Order], order_mask: &[bool], + xir: f64, xif: f64, ) -> Result { // TODO: skip empty subgrids @@ -312,11 +311,11 @@ fn ndarray_from_subgrid_orders_slice( let mut logs = 1.0; if order.logxir > 0 { - if approx_eq!(f64, info.xir, 1.0, ulps = 4) { + if approx_eq!(f64, xir, 1.0, ulps = 4) { continue; } - logs *= (info.xir * info.xir).ln(); + logs *= (xir * xir).ln(); } if order.logxif > 0 { @@ -358,7 +357,7 @@ fn ndarray_from_subgrid_orders_slice( continue; } - let mur2 = info.xir * info.xir * ren; + let mur2 = xir * xir * ren; let als = if order.alphas == 0 { 1.0 @@ -389,7 +388,7 @@ pub(crate) fn evolve_slice_with_one( operator: &ArrayView4, info: &OperatorSliceInfo, order_mask: &[bool], - xif: f64, + (xir, xif): (f64, f64), ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); let has_pdf1 = grid.has_pdf1(); @@ -417,6 +416,7 @@ pub(crate) fn evolve_slice_with_one( &subgrids_o, grid.orders(), order_mask, + xir, xif, )?; @@ -510,7 +510,7 @@ pub(crate) fn evolve_slice_with_two( operator: &ArrayView4, info: &OperatorSliceInfo, order_mask: &[bool], - xif: f64, + (xir, xif): (f64, f64), ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); @@ -544,6 +544,7 @@ pub(crate) fn evolve_slice_with_two( &subgrids_o, grid.orders(), order_mask, + xir, xif, )?; diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 82c8d391..bc823f7f 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -1926,14 +1926,13 @@ impl Grid { x1: info.x1.clone(), ren1: info.ren1.clone(), alphas: info.alphas.clone(), - xir: info.xir, lumi_id_types: info.lumi_id_types.clone(), }, CowArray::from(op), )) }), order_mask, - info.xif, + (info.xir, info.xif), ) } @@ -1956,7 +1955,7 @@ impl Grid { &self, slices: impl IntoIterator), E>>, order_mask: &[bool], - xif: f64, + xi: (f64, f64), ) -> Result { use super::evolution::EVOLVE_INFO_TOL_ULPS; @@ -1984,9 +1983,9 @@ impl Grid { let view = operator.view(); let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { - evolution::evolve_slice_with_two(self, &view, &info, order_mask, xif) + evolution::evolve_slice_with_two(self, &view, &info, order_mask, xi) } else { - evolution::evolve_slice_with_one(self, &view, &info, order_mask, xif) + evolution::evolve_slice_with_one(self, &view, &info, order_mask, xi) }?; let mut rhs = Self { diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index cb3f5d97..e4dde369 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -127,13 +127,13 @@ mod eko { Err(anyhow!("no file 'metadata.yaml' in EKO archive found")) } - pub fn new(eko_path: &Path, ren1: &[f64], alphas: &[f64], xir: f64) -> Result { + pub fn new(eko_path: &Path, ren1: &[f64], alphas: &[f64]) -> Result { let metadata = Self::read_metadata(eko_path)?; match metadata { - Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas, xir), - Metadata::V1(v1) => Self::with_v1(v1, eko_path, ren1, alphas, xir), - Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas, xir), + Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas), + Metadata::V1(v1) => Self::with_v1(v1, eko_path, ren1, alphas), + Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas), } } @@ -142,7 +142,6 @@ mod eko { eko_path: &Path, ren1: &[f64], alphas: &[f64], - xir: f64, ) -> Result { let mut operator = None; @@ -170,7 +169,6 @@ mod eko { x1: metadata.targetgrid, ren1: ren1.to_vec(), alphas: alphas.to_vec(), - xir, }, operator, }) @@ -181,7 +179,6 @@ mod eko { eko_path: &Path, ren1: &[f64], alphas: &[f64], - xir: f64, ) -> Result { let mut fac1 = HashMap::new(); let base64 = GeneralPurpose::new(&URL_SAFE, PAD); @@ -255,7 +252,6 @@ mod eko { .unwrap_or(metadata.rotations.xgrid), ren1: ren1.to_vec(), alphas: alphas.to_vec(), - xir, }, archive: Archive::new(File::open(eko_path)?), }) @@ -266,7 +262,6 @@ mod eko { eko_path: &Path, ren1: &[f64], alphas: &[f64], - xir: f64, ) -> Result { let mut fac1 = HashMap::new(); let mut operator: Option = None; @@ -332,7 +327,6 @@ mod eko { .unwrap_or_else(|| metadata.bases.xgrid.clone()), ren1: ren1.to_vec(), alphas: alphas.to_vec(), - xir, }, archive: Archive::new(File::open(eko_path)?), }) @@ -495,9 +489,9 @@ fn evolve_grid( }) .collect(); - let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas, xir)?; + let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas)?; - Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, xif)?) + Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, (xir, xif))?) } #[cfg(not(feature = "evolve"))] From 3cfe6dcf987d5332805a16861375e17ed12b0093 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 13:51:42 +0100 Subject: [PATCH 29/40] Pack `xir` and `xif` into a tuple --- pineappl/src/evolution.rs | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index cf2ee2fe..54cde746 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -272,8 +272,7 @@ fn ndarray_from_subgrid_orders_slice( subgrids: &ArrayView1, orders: &[Order], order_mask: &[bool], - xir: f64, - xif: f64, + (xir, xif): (f64, f64), ) -> Result { // TODO: skip empty subgrids @@ -388,7 +387,7 @@ pub(crate) fn evolve_slice_with_one( operator: &ArrayView4, info: &OperatorSliceInfo, order_mask: &[bool], - (xir, xif): (f64, f64), + xi: (f64, f64), ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); let has_pdf1 = grid.has_pdf1(); @@ -416,8 +415,7 @@ pub(crate) fn evolve_slice_with_one( &subgrids_o, grid.orders(), order_mask, - xir, - xif, + xi, )?; let x1 = if has_pdf1 { x1_a } else { x1_b }; @@ -510,7 +508,7 @@ pub(crate) fn evolve_slice_with_two( operator: &ArrayView4, info: &OperatorSliceInfo, order_mask: &[bool], - (xir, xif): (f64, f64), + xi: (f64, f64), ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); @@ -544,8 +542,7 @@ pub(crate) fn evolve_slice_with_two( &subgrids_o, grid.orders(), order_mask, - xir, - xif, + xi, )?; if (last_x1a.len() != x1_a.len()) From 6a86a88b578d953f3e8b1d352bf412e375b21283 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 14:10:31 +0100 Subject: [PATCH 30/40] Introduce new type `AlphasTable` --- pineappl/src/evolution.rs | 54 ++++++++++++++++++++++++++++++------ pineappl/src/grid.rs | 13 +++++---- pineappl_cli/src/evolve.rs | 57 ++++++++------------------------------ 3 files changed, 64 insertions(+), 60 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 54cde746..c708aaa0 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -134,12 +134,43 @@ pub struct OperatorSliceInfo { /// `x`-grid coordinates of the `Grid`. pub x1: Vec, + /// Identifier of the particle basis for the `FkTable`. + pub lumi_id_types: String, +} + +/// A mapping of squared renormalization scales in `ren1` to strong couplings in `alphas`. The +/// ordering of both members defines the mapping. +pub struct AlphasTable { /// Renormalization scales of the `Grid`. pub ren1: Vec, /// Strong couplings corresponding to the order given in [`ren1`](Self::ren1). pub alphas: Vec, - /// Identifier of the particle basis for the `FkTable`. - pub lumi_id_types: String, +} + +impl AlphasTable { + /// Create an `AlphasTable` for `grid`, varying the renormalization scale by `xir` for the + /// strong couplings given by `alphas`. The only argument of `alphas` must be the squared + /// renormalization scale. + pub fn from_grid(grid: &Grid, xir: f64, alphas: &dyn Fn(f64) -> f64) -> Self { + let mut ren1: Vec<_> = grid + .subgrids() + .iter() + .flat_map(|subgrid| { + subgrid + .mu2_grid() + .iter() + .map(|Mu2 { ren, .. }| xir * xir * ren) + .collect::>() + }) + .collect(); + // UNWRAP: if we can't sort numbers the grid is fishy + ren1.sort_by(|a, b| a.partial_cmp(b).unwrap()); + ren1.dedup(); + let ren1 = ren1; + let alphas: Vec<_> = ren1.iter().map(|&mur2| alphas(mur2)).collect(); + + Self { ren1, alphas } + } } fn gluon_has_pid_zero(grid: &Grid) -> bool { @@ -273,6 +304,7 @@ fn ndarray_from_subgrid_orders_slice( orders: &[Order], order_mask: &[bool], (xir, xif): (f64, f64), + alphas_table: &AlphasTable, ) -> Result { // TODO: skip empty subgrids @@ -360,13 +392,13 @@ fn ndarray_from_subgrid_orders_slice( let als = if order.alphas == 0 { 1.0 - } else if let Some(alphas) = - info.ren1 - .iter() - .zip(info.alphas.iter()) - .find_map(|(&ren1, &alphas)| { - approx_eq!(f64, ren1, mur2, ulps = EVOLUTION_TOL_ULPS).then(|| alphas) - }) + } else if let Some(alphas) = alphas_table + .ren1 + .iter() + .zip(alphas_table.alphas.iter()) + .find_map(|(&ren1, &alphas)| { + approx_eq!(f64, ren1, mur2, ulps = EVOLUTION_TOL_ULPS).then(|| alphas) + }) { alphas.powi(order.alphas.try_into().unwrap()) } else { @@ -388,6 +420,7 @@ pub(crate) fn evolve_slice_with_one( info: &OperatorSliceInfo, order_mask: &[bool], xi: (f64, f64), + alphas_table: &AlphasTable, ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); let has_pdf1 = grid.has_pdf1(); @@ -416,6 +449,7 @@ pub(crate) fn evolve_slice_with_one( grid.orders(), order_mask, xi, + alphas_table, )?; let x1 = if has_pdf1 { x1_a } else { x1_b }; @@ -509,6 +543,7 @@ pub(crate) fn evolve_slice_with_two( info: &OperatorSliceInfo, order_mask: &[bool], xi: (f64, f64), + alphas_table: &AlphasTable, ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); @@ -543,6 +578,7 @@ pub(crate) fn evolve_slice_with_two( grid.orders(), order_mask, xi, + alphas_table, )?; if (last_x1a.len() != x1_a.len()) diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index bc823f7f..5563965a 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -2,7 +2,7 @@ use super::bin::{BinInfo, BinLimits, BinRemapper}; use super::empty_subgrid::EmptySubgridV1; -use super::evolution::{self, EvolveInfo, OperatorInfo, OperatorSliceInfo}; +use super::evolution::{self, AlphasTable, EvolveInfo, OperatorInfo, OperatorSliceInfo}; use super::fk_table::FkTable; use super::import_only_subgrid::ImportOnlySubgridV2; use super::lagrange_subgrid::{LagrangeSparseSubgridV1, LagrangeSubgridV1, LagrangeSubgridV2}; @@ -1924,8 +1924,6 @@ impl Grid { fac1, pids1: info.pids1.clone(), x1: info.x1.clone(), - ren1: info.ren1.clone(), - alphas: info.alphas.clone(), lumi_id_types: info.lumi_id_types.clone(), }, CowArray::from(op), @@ -1933,6 +1931,10 @@ impl Grid { }), order_mask, (info.xir, info.xif), + &AlphasTable { + ren1: info.ren1.clone(), + alphas: info.alphas.clone(), + }, ) } @@ -1956,6 +1958,7 @@ impl Grid { slices: impl IntoIterator), E>>, order_mask: &[bool], xi: (f64, f64), + alphas_table: &AlphasTable, ) -> Result { use super::evolution::EVOLVE_INFO_TOL_ULPS; @@ -1983,9 +1986,9 @@ impl Grid { let view = operator.view(); let (subgrids, lumi) = if self.has_pdf1() && self.has_pdf2() { - evolution::evolve_slice_with_two(self, &view, &info, order_mask, xi) + evolution::evolve_slice_with_two(self, &view, &info, order_mask, xi, alphas_table) } else { - evolution::evolve_slice_with_one(self, &view, &info, order_mask, xi) + evolution::evolve_slice_with_one(self, &view, &info, order_mask, xi, alphas_table) }?; let mut rhs = Self { diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index e4dde369..c781bbca 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -127,22 +127,17 @@ mod eko { Err(anyhow!("no file 'metadata.yaml' in EKO archive found")) } - pub fn new(eko_path: &Path, ren1: &[f64], alphas: &[f64]) -> Result { + pub fn new(eko_path: &Path) -> Result { let metadata = Self::read_metadata(eko_path)?; match metadata { - Metadata::V0(v0) => Self::with_v0(v0, eko_path, ren1, alphas), - Metadata::V1(v1) => Self::with_v1(v1, eko_path, ren1, alphas), - Metadata::V2(v2) => Self::with_v2(v2, eko_path, ren1, alphas), + Metadata::V0(v0) => Self::with_v0(v0, eko_path), + Metadata::V1(v1) => Self::with_v1(v1, eko_path), + Metadata::V2(v2) => Self::with_v2(v2, eko_path), } } - fn with_v0( - metadata: MetadataV0, - eko_path: &Path, - ren1: &[f64], - alphas: &[f64], - ) -> Result { + fn with_v0(metadata: MetadataV0, eko_path: &Path) -> Result { let mut operator = None; for entry in Archive::new(File::open(eko_path)?).entries_with_seek()? { @@ -167,19 +162,12 @@ mod eko { fac1: 0.0, pids1: metadata.targetpids, x1: metadata.targetgrid, - ren1: ren1.to_vec(), - alphas: alphas.to_vec(), }, operator, }) } - fn with_v1( - metadata: MetadataV1, - eko_path: &Path, - ren1: &[f64], - alphas: &[f64], - ) -> Result { + fn with_v1(metadata: MetadataV1, eko_path: &Path) -> Result { let mut fac1 = HashMap::new(); let base64 = GeneralPurpose::new(&URL_SAFE, PAD); @@ -250,19 +238,12 @@ mod eko { .rotations .targetgrid .unwrap_or(metadata.rotations.xgrid), - ren1: ren1.to_vec(), - alphas: alphas.to_vec(), }, archive: Archive::new(File::open(eko_path)?), }) } - fn with_v2( - metadata: MetadataV2, - eko_path: &Path, - ren1: &[f64], - alphas: &[f64], - ) -> Result { + fn with_v2(metadata: MetadataV2, eko_path: &Path) -> Result { let mut fac1 = HashMap::new(); let mut operator: Option = None; @@ -325,8 +306,6 @@ mod eko { .bases .targetgrid .unwrap_or_else(|| metadata.bases.xgrid.clone()), - ren1: ren1.to_vec(), - alphas: alphas.to_vec(), }, archive: Archive::new(File::open(eko_path)?), }) @@ -460,23 +439,9 @@ fn evolve_grid( xif: f64, ) -> Result { use eko::EkoSlices; - use pineappl::subgrid::{Mu2, Subgrid}; + use pineappl::evolution::AlphasTable; - let mut ren1: Vec<_> = grid - .subgrids() - .iter() - .flat_map(|subgrid| { - subgrid - .mu2_grid() - .iter() - .map(|Mu2 { ren, .. }| xir * xir * ren) - .collect::>() - }) - .collect(); - ren1.sort_by(|a, b| a.partial_cmp(b).unwrap()); - ren1.dedup(); - let ren1 = ren1; - let alphas: Vec<_> = ren1.iter().map(|&mur2| pdf.alphas_q2(mur2)).collect(); + let alphas_table = AlphasTable::from_grid(grid, xir, &|q2| pdf.alphas_q2(q2)); let order_mask: Vec<_> = grid .orders() @@ -489,9 +454,9 @@ fn evolve_grid( }) .collect(); - let mut eko_slices = EkoSlices::new(eko, &ren1, &alphas)?; + let mut eko_slices = EkoSlices::new(eko)?; - Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, (xir, xif))?) + Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, (xir, xif), &alphas_table)?) } #[cfg(not(feature = "evolve"))] From e556a3b5e74cd7f894e4534b30fd8d8dc1eab5e7 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 14:22:51 +0100 Subject: [PATCH 31/40] Fix documentation of `OperatorSliceInfo` --- pineappl/src/evolution.rs | 42 +++++++++++++-------------------------- 1 file changed, 14 insertions(+), 28 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index c708aaa0..3fa837ed 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -89,35 +89,21 @@ pub struct OperatorInfo { pub lumi_id_types: String, } -/// Information about the evolution kernel operator slice (EKO) passed to [`Grid::evolve_slice`] as -/// `operator`, which is used to convert a [`Grid`] into an [`FkTable`]. The dimensions of the EKO -/// must correspond to the values given in [`fac1`], [`pids0`], [`x0`], [`pids1`] and [`x1`], -/// exactly in this order. Members with a `1` are defined at the squared factorization scale given -/// as [`fac1`] (often called process scale) and are found in the [`Grid`] that -/// [`Grid::evolve_slice`] is called with. Members with a `0` are defined at the squared -/// factorization scale [`fac0`] (often called fitting scale or starting scale) and are found in -/// the [`FkTable`] resulting from [`Grid::evolve`]. +/// Information about the evolution kernel operator slice (EKO) passed to +/// [`Grid::evolve_with_slice_iter`](super::grid::Grid::evolve_with_slice_iter) as `operator`, +/// which is used to convert a [`Grid`] into an [`FkTable`](super::fk_table::FkTable). The +/// dimensions of the EKO must correspond to the values given in [`fac1`](Self::fac1), +/// [`pids0`](Self::pids0), [`x0`](Self::x0), [`pids1`](Self::pids1) and [`x1`](Self::x1), exactly +/// in this order. Members with a `1` are defined at the squared factorization scale given as +/// `fac1` (often called process scale) and are found in the [`Grid`] that +/// `Grid::evolve_with_slice_iter` is called with. Members with a `0` are defined at the squared +/// factorization scale [`fac0`](Self::fac0) (often called fitting scale or starting scale) and are +/// found in the `FkTable` resulting from [`Grid::evolve`]. /// -/// The EKO slice may convert a `Grid` from a basis given by the particle identifiers [`pids1`] to -/// a possibly different basis given by [`pids0`]. This basis must also be identified using -/// [`lumi_id_types`], which tells [`FkTable::convolute`] how to perform a convolution. The members -/// [`ren1`] and [`alphas`] must be the strong couplings given at the respective renormalization -/// scales. Finally, [`xir`] and [`xif`] can be used to vary the renormalization and factorization -/// scales, respectively, around their central values. -/// -/// [`FkTable::convolute`]: super::fk_table::FkTable::convolute -/// [`FkTable`]: super::fk_table::FkTable -/// [`alphas`]: Self::alphas -/// [`fac0`]: Self::fac0 -/// [`fac1`]: Self::fac1 -/// [`lumi_id_types`]: Self::lumi_id_types -/// [`pids0`]: Self::pids0 -/// [`pids1`]: Self::pids1 -/// [`ren1`]: Self::ren1 -/// [`x0`]: Self::x0 -/// [`x1`]: Self::x1 -/// [`xif`]: Self::xif -/// [`xir`]: Self::xir +/// The EKO slice may convert a `Grid` from a basis given by the particle identifiers `pids1` to a +/// possibly different basis given by `pids0`. This basis must also be identified using +/// [`lumi_id_types`](Self::lumi_id_types), which tells +/// [`FkTable::convolute`](super::fk_table::FkTable::convolute) how to perform a convolution. #[derive(Clone)] pub struct OperatorSliceInfo { /// Squared factorization scale of the `FkTable`. From 2a56cce3e278ca248297cf569c24e740dcdbae59 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 15:20:54 +0100 Subject: [PATCH 32/40] Fix two clippy warnings --- pineappl/src/evolution.rs | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 3fa837ed..61104eca 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -150,7 +150,7 @@ impl AlphasTable { }) .collect(); // UNWRAP: if we can't sort numbers the grid is fishy - ren1.sort_by(|a, b| a.partial_cmp(b).unwrap()); + ren1.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!())); ren1.dedup(); let ren1 = ren1; let alphas: Vec<_> = ren1.iter().map(|&mur2| alphas(mur2)).collect(); @@ -168,8 +168,7 @@ fn gluon_has_pid_zero(grid: &Grid) -> bool { && grid .key_values() .and_then(|key_values| key_values.get("lumi_id_types")) - .map(|value| value == "pdg_mc_ids") - .unwrap_or(true) + .map_or(true, |value| value == "pdg_mc_ids") } type Pid01IndexTuples = Vec<(usize, usize)>; From 5ae6f9c99a5743f7505e6843841f86241bbc40b2 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 16:35:07 +0100 Subject: [PATCH 33/40] Fix and test evolution with scale variations --- pineappl/src/evolution.rs | 3 +- pineappl/src/grid.rs | 20 +++-- pineappl_cli/tests/evolve.rs | 143 +++++++++++++++++++++++++++++++++++ 3 files changed, 156 insertions(+), 10 deletions(-) diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 61104eca..070e96b4 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -293,7 +293,6 @@ fn ndarray_from_subgrid_orders_slice( ) -> Result { // TODO: skip empty subgrids - let fac1 = xif * xif * info.fac1; let mut x1_a: Vec<_> = subgrids .iter() .enumerate() @@ -369,7 +368,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, fac1, ulps = EVOLUTION_TOL_ULPS) { + if !approx_eq!(f64, xif * xif * fac, info.fac1, ulps = EVOLUTION_TOL_ULPS) { continue; } diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 5563965a..39afad04 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -2015,16 +2015,20 @@ impl Grid { // UNWRAP: if we can't compare two numbers there's a bug fac1.sort_by(|a, b| a.partial_cmp(b).unwrap()); - // TODO: here there's possibly is a bug if xif isn't one - // make sure we've evolved all slices - if let Some(muf2) = self.evolve_info(&order_mask).fac1.into_iter().find(|&x| { - !fac1 - .iter() - .any(|&y| approx_eq!(f64, x, y, ulps = EVOLVE_INFO_TOL_ULPS)) - }) { + 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:#?}" + "no operator for muf2 = {muf2} found in {fac1:?}" ))); } diff --git a/pineappl_cli/tests/evolve.rs b/pineappl_cli/tests/evolve.rs index 42a4b717..405759a6 100644 --- a/pineappl_cli/tests/evolve.rs +++ b/pineappl_cli/tests/evolve.rs @@ -76,6 +76,34 @@ const LHCB_WP_7TEV_V2_STR: &str = "b Grid FkTable 7 2.9272102213930090e1 2.9268443366651141e1 -1.2499434622803562e-4 "; +const LHCB_WP_7TEV_V2_XIR2_STR: &str = "b Grid FkTable rel. diff +-+--------------------+--------------------+---------------------- +0 7.7634833292737017e2 7.7614037816519419e2 -2.6786270203205120e-4 +1 7.0866199875124983e2 7.0847444839781781e2 -2.6465417048249229e-4 +2 6.1427556024981789e2 6.1411417374531106e2 -2.6272655946324441e-4 +3 4.9482819982783724e2 4.9469964081143053e2 -2.5980535557890150e-4 +4 3.6756257449354945e2 3.6746967569489709e2 -2.5274281196974169e-4 +5 2.4912642701834142e2 2.4906651029915440e2 -2.4050727939273209e-4 +6 1.1776254040032327e2 1.1773772039493417e2 -2.1076316207790935e-4 +7 2.8749891297668260e1 2.8746299479656258e1 -1.2493327278395583e-4 +"; + +const LHCB_WP_7TEV_V2_XIF_2_STR: &str = + "b Grid FkTable rel. diff +-+--------------------+--------------------+---------------------- +0 8.0902449713533758e2 8.0880109089579207e2 -2.7614273774967391e-4 +1 7.3869242569893402e2 7.3849113100483919e2 -2.7250136469769703e-4 +2 6.4102495904778243e2 6.4085178025871448e2 -2.7015919836448354e-4 +3 5.1668563837653949e2 5.1654786167667771e2 -2.6665478896348294e-4 +4 3.8405066991124284e2 3.8395127677619655e2 -2.5880213949180941e-4 +5 2.6047697125229388e2 2.6041295913273854e2 -2.4574963094659008e-4 +6 1.2324364745022301e2 1.2321715784184289e2 -2.1493690691698486e-4 +7 3.0134629982656573e1 3.0130872371345841e1 -1.2469412476256991e-4 +"; + +const LHCB_WP_7TEV_V2_XIF_2_ERROR_STR: &str = "Error: failed to evolve grid: no operator for muf2 = 25825.775616000003 found in [6456.443904000001] +"; + const LHCB_WP_7TEV_OPTIMIZED_STR: &str = "b etal dsig/detal [] [pb] -+----+----+----------- @@ -237,6 +265,121 @@ fn lhcb_wp_7tev_v2() { .stdout(LHCB_WP_7TEV_V2_STR); } +#[test] +fn lhcb_wp_7tev_v2_xir_2() { + let input = NamedTempFile::new("optimized.pineappl.lz4").unwrap(); + + // we first need to optimize the grid, to strip empty x-grid values not contained in the EKO + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "write", + "--optimize", + "../test-data/LHCB_WP_7TEV.pineappl.lz4", + input.path().to_str().unwrap(), + ]) + .assert() + .success() + .stdout(""); + + let output = NamedTempFile::new("fktable2b.lz4").unwrap(); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--silence-lhapdf", + "evolve", + "--digits-abs=16", + "--digits-rel=16", + input.path().to_str().unwrap(), + "../test-data/LHCB_WP_7TEV_v2.tar", + output.path().to_str().unwrap(), + "NNPDF40_nlo_as_01180", + "--orders=a2,as1a2", + "--xir=2", + ]) + .assert() + .success() + .stdout(LHCB_WP_7TEV_V2_XIR2_STR); +} + +#[test] +fn lhcb_wp_7tev_v2_xif_2() { + let input = NamedTempFile::new("optimized.pineappl.lz4").unwrap(); + + // we first need to optimize the grid, to strip empty x-grid values not contained in the EKO + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "write", + "--optimize", + "../test-data/LHCB_WP_7TEV.pineappl.lz4", + input.path().to_str().unwrap(), + ]) + .assert() + .success() + .stdout(""); + + let output = NamedTempFile::new("fktable2c.lz4").unwrap(); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--silence-lhapdf", + "evolve", + "--digits-abs=16", + "--digits-rel=16", + input.path().to_str().unwrap(), + "../test-data/myeko.tar", + output.path().to_str().unwrap(), + "NNPDF40_nlo_as_01180", + "--orders=a2,as1a2", + "--xif=2", + ]) + .assert() + .success() + .stdout(LHCB_WP_7TEV_V2_XIF_2_STR); +} + +#[test] +fn lhcb_wp_7tev_v2_xif_2_error() { + let input = NamedTempFile::new("optimized.pineappl.lz4").unwrap(); + + // we first need to optimize the grid, to strip empty x-grid values not contained in the EKO + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "write", + "--optimize", + "../test-data/LHCB_WP_7TEV.pineappl.lz4", + input.path().to_str().unwrap(), + ]) + .assert() + .success() + .stdout(""); + + let output = NamedTempFile::new("fktable2c.lz4").unwrap(); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--silence-lhapdf", + "evolve", + "--digits-abs=16", + "--digits-rel=16", + input.path().to_str().unwrap(), + "../test-data/LHCB_WP_7TEV_v2.tar", + output.path().to_str().unwrap(), + "NNPDF40_nlo_as_01180", + "--orders=a2,as1a2", + "--xif=2", + ]) + .assert() + .failure() + .stderr(LHCB_WP_7TEV_V2_XIF_2_ERROR_STR) + .stdout(""); +} + #[test] fn e906nlo_bin_00() { let input = NamedTempFile::new("E906nlo_bin_00_unique_bin_limits.pineappl.lz4").unwrap(); From c6d9e76c8ace36c837fdebe41ceed190b4f26cb0 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 16:39:55 +0100 Subject: [PATCH 34/40] Add missing test data --- .github/workflows/rust.yml | 3 ++- pineappl_cli/tests/evolve.rs | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index a61db8ab..2f90f0b4 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -23,7 +23,7 @@ jobs: uses: actions/cache@v3 with: path: test-data - key: test-data-v10 + key: test-data-v11 - name: Download test data if: steps.cache-test-data.outputs.cache-hit != 'true' run: | @@ -43,6 +43,7 @@ jobs: curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.tar' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_v2.tar' + curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_v2_xif_2.tar' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/NJetEvents_0-0-2.tab.gz' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/NUTEV_CC_NU_FE_SIGMARED.pineappl.lz4' curl -s -C - -O 'https://data.nnpdf.science/pineappl/test-data/NUTEV_CC_NU_FE_SIGMARED.tar' diff --git a/pineappl_cli/tests/evolve.rs b/pineappl_cli/tests/evolve.rs index 405759a6..2a03636b 100644 --- a/pineappl_cli/tests/evolve.rs +++ b/pineappl_cli/tests/evolve.rs @@ -330,7 +330,7 @@ fn lhcb_wp_7tev_v2_xif_2() { "--digits-abs=16", "--digits-rel=16", input.path().to_str().unwrap(), - "../test-data/myeko.tar", + "../test-data/LHCB_WP_7TEV_v2_xif_2.tar", output.path().to_str().unwrap(), "NNPDF40_nlo_as_01180", "--orders=a2,as1a2", From 0292ac89f8c7ced44ef4f4a30e9a599808aab254 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 3 Mar 2024 17:00:28 +0100 Subject: [PATCH 35/40] Test old method `Grid::evolve` --- pineappl_cli/src/evolve.rs | 57 ++++++++++++++++++++++++++++++++---- pineappl_cli/tests/evolve.rs | 21 +++++++++++++ 2 files changed, 72 insertions(+), 6 deletions(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index c781bbca..555e10a7 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -437,11 +437,10 @@ fn evolve_grid( orders: &[(u32, u32)], xir: f64, xif: f64, + use_old_evolve: bool, ) -> Result { use eko::EkoSlices; - use pineappl::evolution::AlphasTable; - - let alphas_table = AlphasTable::from_grid(grid, xir, &|q2| pdf.alphas_q2(q2)); + use pineappl::evolution::{AlphasTable, OperatorInfo}; let order_mask: Vec<_> = grid .orders() @@ -455,12 +454,48 @@ fn evolve_grid( .collect(); let mut eko_slices = EkoSlices::new(eko)?; + let alphas_table = AlphasTable::from_grid(grid, xir, &|q2| pdf.alphas_q2(q2)); - Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, (xir, xif), &alphas_table)?) + if use_old_evolve { + if let EkoSlices::V0 { + fac1, + info, + operator, + } = eko_slices + { + let op_info = OperatorInfo { + fac0: info.fac0, + pids0: info.pids0.clone(), + x0: info.x0.clone(), + fac1: fac1.clone(), + pids1: info.pids1.clone(), + x1: info.x1.clone(), + ren1: alphas_table.ren1, + alphas: alphas_table.alphas, + xir, + xif, + lumi_id_types: info.lumi_id_types, + }; + + Ok(grid.evolve(operator.view(), &op_info, &order_mask)?) + } else { + unimplemented!(); + } + } else { + Ok(grid.evolve_with_slice_iter(&mut eko_slices, &order_mask, (xir, xif), &alphas_table)?) + } } #[cfg(not(feature = "evolve"))] -fn evolve_grid(_: &Grid, _: &Path, _: &Pdf, _: &[(u32, u32)], _: f64, _: f64) -> Result { +fn evolve_grid( + _: &Grid, + _: &Path, + _: &Pdf, + _: &[(u32, u32)], + _: f64, + _: f64, + _: bool, +) -> Result { Err(anyhow!( "you need to install `pineappl` with feature `evolve`" )) @@ -505,6 +540,8 @@ pub struct Opts { /// Rescale the factorization scale with this factor. #[arg(default_value_t = 1.0, long)] xif: f64, + #[arg(hide = true, long)] + use_old_evolve: bool, } impl Subcommand for Opts { @@ -524,7 +561,15 @@ impl Subcommand for Opts { cfg.force_positive, ); - let fk_table = evolve_grid(&grid, &self.eko, &pdf, &self.orders, self.xir, self.xif)?; + let fk_table = evolve_grid( + &grid, + &self.eko, + &pdf, + &self.orders, + self.xir, + self.xif, + self.use_old_evolve, + )?; let evolved_results = helpers::convolute_scales( fk_table.grid(), &mut pdf, diff --git a/pineappl_cli/tests/evolve.rs b/pineappl_cli/tests/evolve.rs index 2a03636b..65ca27bd 100644 --- a/pineappl_cli/tests/evolve.rs +++ b/pineappl_cli/tests/evolve.rs @@ -228,6 +228,27 @@ fn lhcb_wp_7tev() { .stdout(LHCB_WP_7TEV_OPTIMIZED_STR); } +#[test] +fn lhcb_wp_7tev_use_old_evolve() { + let output = NamedTempFile::new("fktable1c.lz4").unwrap(); + + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--silence-lhapdf", + "evolve", + "../test-data/LHCB_WP_7TEV.pineappl.lz4", + "../test-data/LHCB_WP_7TEV.tar", + output.path().to_str().unwrap(), + "NNPDF40_nlo_as_01180", + "--orders=a2,as1a2", + "--use-old-evolve", + ]) + .assert() + .success() + .stdout(LHCB_WP_7TEV_STR); +} + #[test] fn lhcb_wp_7tev_v2() { let input = NamedTempFile::new("optimized.pineappl.lz4").unwrap(); From 91246246b406d18923b66d00c67c0ac4f7332092 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 4 Mar 2024 10:18:10 +0100 Subject: [PATCH 36/40] Mark `Grid::evolve` deprecated --- CHANGELOG.md | 4 ++++ pineappl/src/grid.rs | 1 + pineappl_cli/src/evolve.rs | 1 + 3 files changed, 6 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6b222190..1be1470b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - added support for new EKO format introduced with EKO v0.13 +### Changed + +- `Grid::evolve` has now been marked deprecated + ## [0.6.2] - 09/10/2023 ### Added diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 39afad04..4b5fbdde 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -1905,6 +1905,7 @@ impl Grid { /// /// Returns a [`GridError::EvolutionFailure`] if either the `operator` or its `info` is /// incompatible with this `Grid`. + #[deprecated(since = "0.7.4", note = "use evolve_with_slice_iter instead")] pub fn evolve( &self, operator: ArrayView5, diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index 555e10a7..c39ff006 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -477,6 +477,7 @@ fn evolve_grid( lumi_id_types: info.lumi_id_types, }; + #[allow(deprecated)] Ok(grid.evolve(operator.view(), &op_info, &order_mask)?) } else { unimplemented!(); From 58972f20655eea49d38e05870ac1986cb4b72cf0 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 4 Mar 2024 10:21:37 +0100 Subject: [PATCH 37/40] Improve changelog entry --- CHANGELOG.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1be1470b..eb8620f8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- added support for new EKO format introduced with EKO v0.13 +- added `Grid::evolve_with_slice_iter`, `AlphasTable` and `OperatorSliceInfo`, + which define a new interface supporting very large evolution kernels that + have been introduced in EKO v0.13. This interface will replace `Grid::evolve` ### Changed From 90d9efdda34092786a59f1cb5c4091575acc9a72 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 4 Mar 2024 10:28:15 +0100 Subject: [PATCH 38/40] Add tentative Python interface for `Grid::evolve_with_slice_iter` --- pineappl_py/src/evolution.rs | 10 +++++++- pineappl_py/src/grid.rs | 46 ++++++++++++++++++++++++++++++++++-- 2 files changed, 53 insertions(+), 3 deletions(-) diff --git a/pineappl_py/src/evolution.rs b/pineappl_py/src/evolution.rs index b554aa1d..c1c83992 100644 --- a/pineappl_py/src/evolution.rs +++ b/pineappl_py/src/evolution.rs @@ -1,5 +1,5 @@ use numpy::{IntoPyArray, PyArray1}; -use pineappl::evolution::EvolveInfo; +use pineappl::evolution::{EvolveInfo, OperatorSliceInfo}; use pyo3::prelude::*; @@ -36,3 +36,11 @@ impl PyEvolveInfo { self.evolve_info.ren1.clone().into_pyarray(py) } } + +/// TODO +#[pyclass] +#[repr(transparent)] +#[derive(Clone)] +pub struct PyOperatorSliceInfo { + pub(crate) slice_info: OperatorSliceInfo, +} diff --git a/pineappl_py/src/grid.rs b/pineappl_py/src/grid.rs index 102b65f4..fc7dd54e 100644 --- a/pineappl_py/src/grid.rs +++ b/pineappl_py/src/grid.rs @@ -7,21 +7,25 @@ use pineappl::grid::{EkoInfo, GridAxes}; use pineappl::lumi::LumiCache; use super::bin::PyBinRemapper; -use super::evolution::PyEvolveInfo; +use super::evolution::{PyEvolveInfo, PyOperatorSliceInfo}; use super::fk_table::PyFkTable; use super::lumi::PyLumiEntry; use super::subgrid::{PySubgridEnum, PySubgridParams}; use itertools::izip; -use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1, PyReadonlyArray5}; +use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1, PyReadonlyArray4, PyReadonlyArray5}; use std::collections::HashMap; +use std::error::Error; use std::fs::File; use std::io::BufReader; use std::path::PathBuf; use pyo3::exceptions::PyValueError; use pyo3::prelude::*; +use pyo3::types::{PyIterator, PyTuple}; + +use ndarray::CowArray; /// PyO3 wrapper to :rustdoc:`pineappl::grid::Order ` /// @@ -562,6 +566,44 @@ impl PyGrid { } } + /// TODO + /// + /// Parameters + /// ---------- + /// slices : TODO + /// order_mask : TODO + /// + /// Returns + /// ------- + /// TODO + pub fn evolve_with_slice_iter( + &self, + slices: &PyIterator, + order_mask: PyReadonlyArray1, + ) -> PyResult { + Ok(self + .grid + .evolve_with_slice_iter( + slices.map(|result| { + // TODO: check whether we can avoid the `.unwrap` calls + let any = result.unwrap(); + let tuple = any.downcast::().unwrap(); + let item0 = tuple.get_item(0).unwrap(); + let item1 = tuple.get_item(1).unwrap(); + let slice_info = item0.extract::().unwrap(); + let operator = item1.extract::>().unwrap(); + + // TODO: change `PyErr` into something appropriate + Ok::<_, PyErr>((slice_info.slice_info, CowArray::from(operator.as_array()))) + }), + // TODO: what if it's non-contiguous? + order_mask.as_slice().unwrap(), + ) + .map(|fk_table| PyFkTable { fk_table }) + // TODO: get rid of this `.unwrap` call + .unwrap()) + } + /// Load grid from file. /// /// **Usage:** `pineko`, FKTable generation From 1e417a0f6aace1a69256a61034b1f5d6c36c51eb Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 4 Mar 2024 10:39:54 +0100 Subject: [PATCH 39/40] Fix compilation problems in the previous commit --- pineappl_py/src/grid.rs | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/pineappl_py/src/grid.rs b/pineappl_py/src/grid.rs index fc7dd54e..d50312de 100644 --- a/pineappl_py/src/grid.rs +++ b/pineappl_py/src/grid.rs @@ -1,4 +1,4 @@ -use pineappl::evolution::OperatorInfo; +use pineappl::evolution::{AlphasTable, OperatorInfo}; use pineappl::grid::{Grid, Ntuple, Order}; #[allow(deprecated)] @@ -16,7 +16,6 @@ use itertools::izip; use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1, PyReadonlyArray4, PyReadonlyArray5}; use std::collections::HashMap; -use std::error::Error; use std::fs::File; use std::io::BufReader; use std::path::PathBuf; @@ -580,6 +579,9 @@ impl PyGrid { &self, slices: &PyIterator, order_mask: PyReadonlyArray1, + xi: (f64, f64), + ren1: Vec, + alphas: Vec, ) -> PyResult { Ok(self .grid @@ -592,12 +594,16 @@ impl PyGrid { let item1 = tuple.get_item(1).unwrap(); let slice_info = item0.extract::().unwrap(); let operator = item1.extract::>().unwrap(); + // TODO: can we get rid of the `into_owned` call? + let array = CowArray::from(operator.as_array().into_owned()); // TODO: change `PyErr` into something appropriate - Ok::<_, PyErr>((slice_info.slice_info, CowArray::from(operator.as_array()))) + Ok::<_, PyErr>((slice_info.slice_info, array)) }), // TODO: what if it's non-contiguous? order_mask.as_slice().unwrap(), + xi, + &AlphasTable { ren1, alphas }, ) .map(|fk_table| PyFkTable { fk_table }) // TODO: get rid of this `.unwrap` call From ad42065ad5061724fa870b5d153db7005d0b06b4 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 4 Mar 2024 16:38:58 +0100 Subject: [PATCH 40/40] Remove unused method --- pineappl_cli/src/evolve.rs | 8 -------- 1 file changed, 8 deletions(-) diff --git a/pineappl_cli/src/evolve.rs b/pineappl_cli/src/evolve.rs index c39ff006..8d5a7a73 100644 --- a/pineappl_cli/src/evolve.rs +++ b/pineappl_cli/src/evolve.rs @@ -311,14 +311,6 @@ mod eko { }) } - // TODO: we could make this a Cow<'_, [f64]> - pub fn fac1(&self) -> Vec { - match self { - Self::V0 { fac1, .. } => fac1.clone(), - Self::V2 { fac1, .. } => fac1.values().copied().collect(), - } - } - pub fn iter_mut(&mut self) -> EkoSlicesIter { match self { Self::V0 {