diff --git a/pineappl/src/boc.rs b/pineappl/src/boc.rs index 9965459d..b71f9a95 100644 --- a/pineappl/src/boc.rs +++ b/pineappl/src/boc.rs @@ -339,17 +339,17 @@ impl Channel { /// use pineappl::boc::Channel; /// use pineappl::channel; /// - /// let entry = Channel::translate(&channel![103, 11, 1.0], &|evol_id| match evol_id { + /// let entry = Channel::translate(&channel![103, 11, 10], &|evol_id| match evol_id { /// 103 => vec![(2, 1.0), (-2, -1.0), (1, -1.0), (-1, 1.0)], /// _ => vec![(evol_id, 1.0)], /// }); /// - /// assert_eq!(entry, channel![2, 11, 1.0; -2, 11, -1.0; 1, 11, -1.0; -1, 11, 1.0]); + /// assert_eq!(entry, channel![2, 11, 10.0; -2, 11, -10.0; 1, 11, -10.0; -1, 11, 10.0]); /// ``` pub fn translate(entry: &Self, translator: &dyn Fn(i32) -> Vec<(i32, f64)>) -> Self { let mut result = Vec::new(); - for &(pids, factor) in &entry.entry { + for (pids, factor) in &entry.entry { for tuples in pids .iter() .map(|&pid| translator(pid)) @@ -357,7 +357,7 @@ impl Channel { { result.push(( tuples.iter().map(|&(pid, _)| pid).collect(), - tuples.iter().map(|(_, f)| f).product::(), + tuples.iter().map(|(_, f)| factor * f).product::(), )); } } @@ -382,11 +382,20 @@ impl Channel { &self.entry } - // /// Creates a new object with the initial states transposed. - // #[must_use] - // pub fn transpose(&self) -> Self { - // Self::new(self.entry.iter().map(|(a, b, c)| (*b, *a, *c)).collect()) - // } + /// Create a new object with the PIDs at index `i` and `j` transposed. + #[must_use] + pub fn transpose(&self, i: usize, j: usize) -> Self { + Self::new( + self.entry + .iter() + .map(|(pids, c)| { + let mut transposed = pids.clone(); + transposed.swap(i, j); + (transposed, *c) + }) + .collect(), + ) + } /// If `other` is the same channel when only comparing PIDs and neglecting the factors, return /// the number `f1 / f2`, where `f1` is the factor from `self` and `f2` is the factor from @@ -498,7 +507,7 @@ impl FromStr for Channel { #[macro_export] macro_rules! channel { ($a:expr, $b:expr, $factor:expr $(; $c:expr, $d:expr, $fac:expr)*) => { - $crate::boc::Channel::new(vec![($a, $b, $factor), $(($c, $d, $fac)),*]) + $crate::boc::Channel::new(vec![(vec![$a, $b], $factor), $((vec![$c, $d], $fac)),*]) }; } diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 5167bbb7..1bb7433e 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -386,6 +386,7 @@ fn ndarray_from_subgrid_orders_slice( Ok((vec![x1_a, x1_b], (!zero).then_some(array))) } +// TODO: merge this method into evolve_slice_with_two pub(crate) fn evolve_slice_with_one( grid: &Grid, operator: &ArrayView4, @@ -395,18 +396,22 @@ pub(crate) fn evolve_slice_with_one( alphas_table: &AlphasTable, ) -> Result<(Array3, Vec), GridError> { let gluon_has_pid_zero = gluon_has_pid_zero(grid); - let has_pdf1 = grid.convolutions()[0] != Convolution::None; - + // UNWRAP: there must be exactly one convolution that's not None + let index = grid + .convolutions() + .iter() + .position(|c| *c != Convolution::None) + .unwrap(); let (pid_indices, pids) = pid_slices(operator, info, gluon_has_pid_zero, &|pid| { grid.channels() .iter() .flat_map(Channel::entry) - .any(|&(a, b, _)| if has_pdf1 { a } else { b } == pid) + .any(|(pids, _)| pids[index] == pid) })?; let channels0 = channels0_with_one(&pids); let mut sub_fk_tables = Vec::with_capacity(grid.bin_info().bins() * channels0.len()); - let new_axis = if has_pdf1 { 2 } else { 1 }; + let new_axis = 2 - index; let mut last_x1 = Vec::new(); let mut ops = Vec::new(); @@ -429,7 +434,7 @@ pub(crate) fn evolve_slice_with_one( continue; }; - let x1 = if has_pdf1 { x1.remove(0) } else { x1.remove(1) }; + let x1 = x1.remove(index); if x1.is_empty() { continue; @@ -445,12 +450,7 @@ pub(crate) fn evolve_slice_with_one( last_x1 = x1; } - for (&pid1, &factor) in - channel1 - .entry() - .iter() - .map(|(a, b, f)| if has_pdf1 { (a, f) } else { (b, f) }) - { + for (pid1, &factor) in channel1.entry().iter().map(|(pids, f)| (pids[index], f)) { for (fk_table, op) in channels0 .iter() @@ -485,18 +485,22 @@ pub(crate) fn evolve_slice_with_one( 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() }, + if index == 0 { + info.x0.clone() + } else { + vec![1.0] + }, + if index == 0 { + vec![1.0] + } else { + info.x0.clone() + }, ) .into() })); } - let pid = if grid.convolutions()[0] == Convolution::None { - grid.channels()[0].entry()[0].0 - } else { - grid.channels()[0].entry()[0].1 - }; + let pid = grid.channels()[0].entry()[0].0[index]; Ok(( Array1::from_iter(sub_fk_tables) @@ -506,8 +510,8 @@ pub(crate) fn evolve_slice_with_one( .iter() .map(|&a| { channel![ - if has_pdf1 { a } else { pid }, - if has_pdf1 { pid } else { a }, + if index == 0 { a } else { pid }, + if index == 0 { pid } else { a }, 1.0 ] }) @@ -532,12 +536,7 @@ pub(crate) fn evolve_slice_with_two( grid.channels() .iter() .flat_map(Channel::entry) - .any(|tuple| match d { - // TODO: `Channel::entry` should return a tuple of a `Vec` and an `f64` - 0 => tuple.0 == pid1, - 1 => tuple.1 == pid1, - _ => unreachable!(), - }) + .any(|(pids, _)| pids[d] == pid1) }) }) .collect::, _>>()? @@ -596,7 +595,7 @@ pub(crate) fn evolve_slice_with_two( for (pids1, factor) in channel1 .entry() .iter() - .map(|&(pida1, pidb1, factor)| ([pida1, pidb1], factor)) + .map(|(pids1, factor)| ([pids1[0], pids1[1]], factor)) { for (fk_table, ops) in channels0 @@ -615,7 +614,7 @@ pub(crate) fn evolve_slice_with_two( }) { linalg::general_mat_mul(1.0, &array, &ops[1].t(), 0.0, &mut tmp); - linalg::general_mat_mul(factor, ops[0], &tmp, 1.0, fk_table); + linalg::general_mat_mul(*factor, ops[0], &tmp, 1.0, fk_table); } } } @@ -675,12 +674,7 @@ pub(crate) fn evolve_slice_with_two2( grid.channels() .iter() .flat_map(Channel::entry) - .any(|tuple| match d { - // TODO: `Channel::entry` should return a tuple of a `Vec` and an `f64` - 0 => tuple.0 == pid1, - 1 => tuple.1 == pid1, - _ => unreachable!(), - }) + .any(|(pids, _)| pids[d] == pid1) }) }) .collect::, _>>()? @@ -747,7 +741,7 @@ pub(crate) fn evolve_slice_with_two2( for (pids1, factor) in channel1 .entry() .iter() - .map(|&(pida1, pidb1, factor)| ([pida1, pidb1], factor)) + .map(|(pids1, factor)| ([pids1[0], pids1[1]], factor)) { for (fk_table, ops) in channels0 @@ -768,7 +762,7 @@ pub(crate) fn evolve_slice_with_two2( // tmp = array * ops[1]^T linalg::general_mat_mul(1.0, &array, &ops[1].t(), 0.0, &mut tmp); // fk_table += factor * ops[0] * tmp - linalg::general_mat_mul(factor, ops[0], &tmp, 1.0, fk_table); + linalg::general_mat_mul(*factor, ops[0], &tmp, 1.0, fk_table); } } } diff --git a/pineappl/src/fk_table.rs b/pineappl/src/fk_table.rs index b7d04519..6b3bcb12 100644 --- a/pineappl/src/fk_table.rs +++ b/pineappl/src/fk_table.rs @@ -226,11 +226,11 @@ impl FkTable { /// Return the channel definition for this `FkTable`. All factors are `1.0`. #[must_use] - pub fn channels(&self) -> Vec<(i32, i32)> { + pub fn channels(&self) -> Vec> { self.grid .channels() .iter() - .map(|entry| (entry.entry()[0].0, entry.entry()[0].1)) + .map(|entry| entry.entry()[0].0.clone()) .collect() } @@ -390,7 +390,7 @@ impl TryFrom for FkTable { for channel in grid.channels() { let entry = channel.entry(); - if entry.len() != 1 || entry[0].2 != 1.0 { + if entry.len() != 1 || entry[0].1 != 1.0 { return Err(TryFromGridError::InvalidChannel); } } diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 86e744b9..18c36847 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -388,9 +388,10 @@ impl Grid { let mut lumi = 0.0; for entry in channel.entry() { - let xfx1 = lumi_cache.xfx1(entry.0, ix1, imu2); - let xfx2 = lumi_cache.xfx2(entry.1, ix2, imu2); - lumi += xfx1 * xfx2 * entry.2 / (x1 * x2); + debug_assert_eq!(entry.0.len(), 2); + let xfx1 = lumi_cache.xfx1(entry.0[0], ix1, imu2); + let xfx2 = lumi_cache.xfx2(entry.0[1], ix2, imu2); + lumi += xfx1 * xfx2 * entry.1 / (x1 * x2); } let alphas = lumi_cache.alphas(imu2); @@ -454,9 +455,10 @@ impl Grid { let mut lumi = 0.0; for entry in channel.entry() { - let xfx1 = lumi_cache.xfx1(entry.0, ix1, imu2); - let xfx2 = lumi_cache.xfx2(entry.1, ix2, imu2); - lumi += xfx1 * xfx2 * entry.2 / (x1 * x2); + debug_assert_eq!(entry.0.len(), 2); + let xfx1 = lumi_cache.xfx1(entry.0[0], ix1, imu2); + let xfx2 = lumi_cache.xfx2(entry.0[1], ix2, imu2); + lumi += xfx1 * xfx2 * entry.1 / (x1 * x2); } let alphas = lumi_cache.alphas(imu2); @@ -804,11 +806,7 @@ impl Grid { { Some(Ok(pid)) => { let condition = !self.channels().iter().all(|entry| { - entry.entry().iter().all(|&channels| match index { - 1 => channels.0 == pid, - 2 => channels.1 == pid, - _ => unreachable!(), - }) + entry.entry().iter().all(|(pids, _)| pids[index] == pid) }); if condition { @@ -1147,7 +1145,7 @@ impl Grid { // only keep channels that have non-zero factors and for which at least one subgrid is // non-empty for (channel, entry) in self.channels.iter().enumerate() { - if !entry.entry().iter().all(|&(_, _, factor)| factor == 0.0) + if !entry.entry().iter().all(|&(_, factor)| factor == 0.0) && !self .subgrids .slice(s![.., .., channel]) @@ -1196,6 +1194,8 @@ impl Grid { fn symmetrize_channels(&mut self) { let convolutions = self.convolutions(); + // TODO: generalize this method to n convolutions + assert_eq!(convolutions.len(), 2); if convolutions[0] != convolutions[1] { return; } @@ -1205,7 +1205,7 @@ impl Grid { while let Some(index) = indices.pop() { let channel_entry = &self.channels[index]; - if *channel_entry == channel_entry.transpose() { + if *channel_entry == channel_entry.transpose(0, 1) { // check if in all cases the limits are compatible with merging self.subgrids .slice_mut(s![.., .., index]) @@ -1218,7 +1218,7 @@ impl Grid { } else if let Some((j, &other_index)) = indices .iter() .enumerate() - .find(|(_, i)| self.channels[**i] == channel_entry.transpose()) + .find(|(_, i)| self.channels[**i] == channel_entry.transpose(0, 1)) { indices.remove(j); @@ -1289,6 +1289,9 @@ impl Grid { pub fn evolve_info(&self, order_mask: &[bool]) -> EvolveInfo { use super::evolution::EVOLVE_INFO_TOL_ULPS; + // TODO: generalize this method to n convolutions and different EKOs + assert_eq!(self.convolutions().len(), 2); + let has_pdf1 = self.convolutions()[0] != Convolution::None; let has_pdf2 = self.convolutions()[1] != Convolution::None; @@ -1324,10 +1327,20 @@ impl Grid { x1.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = EVOLVE_INFO_TOL_ULPS)); if has_pdf1 { - pids1.extend(self.channels()[channel].entry().iter().map(|(a, _, _)| a)); + pids1.extend( + self.channels()[channel] + .entry() + .iter() + .map(|(pids, _)| pids[0]), + ); } if has_pdf2 { - pids1.extend(self.channels()[channel].entry().iter().map(|(_, b, _)| b)); + pids1.extend( + self.channels()[channel] + .entry() + .iter() + .map(|(pids, _)| pids[1]), + ); } pids1.sort_unstable(); @@ -1753,6 +1766,9 @@ impl Grid { } pub(crate) fn rewrite_channels(&mut self, add: &[(i32, i32)], del: &[i32]) { + // TODO: generalize this method to n convolutions + assert_eq!(self.convolutions().len(), 2); + self.channels = self .channels() .iter() @@ -1761,21 +1777,29 @@ impl Grid { entry .entry() .iter() - .map(|(a, b, f)| { + .map(|(pids, f)| { ( - // if `a` is to be added to another pid replace it with this pid - add.iter().fold( - *a, - |id, &(source, target)| if id == source { target } else { id }, - ), - // if `b` is to be added to another pid replace it with this pid - add.iter().fold( - *b, - |id, &(source, target)| if id == source { target } else { id }, - ), + vec![ + // if `a` is to be added to another pid replace it with this pid + add.iter().fold(pids[0], |id, &(source, target)| { + if id == source { + target + } else { + id + } + }), + // if `b` is to be added to another pid replace it with this pid + add.iter().fold(pids[1], |id, &(source, target)| { + if id == source { + target + } else { + id + } + }), + ], // if any of the pids `a` or `b` are to b deleted set the factor to // zero - if del.iter().any(|id| id == a || id == b) { + if del.iter().any(|&id| id == pids[0] || id == pids[1]) { 0.0 } else { *f @@ -1805,7 +1829,7 @@ impl Grid { entry .entry() .iter() - .copied() + .cloned() .map(move |entry| Channel::new(vec![entry])) }) .collect(); diff --git a/pineappl/src/pids.rs b/pineappl/src/pids.rs index 8e23eaa4..fca5a5b3 100644 --- a/pineappl/src/pids.rs +++ b/pineappl/src/pids.rs @@ -894,9 +894,9 @@ mod tests { ); assert_eq!(result.entry().len(), 1); - assert_eq!(result.entry()[0].0, pid); - assert_eq!(result.entry()[0].1, pid); - assert_approx_eq!(f64, result.entry()[0].2, 1.0, ulps = 8); + assert_eq!(result.entry()[0].0[0], pid); + assert_eq!(result.entry()[0].0[1], pid); + assert_approx_eq!(f64, result.entry()[0].1, 1.0, ulps = 8); } } } diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs index 38c45d2d..cab56b19 100644 --- a/pineappl_capi/src/lib.rs +++ b/pineappl_capi/src/lib.rs @@ -1017,7 +1017,7 @@ pub unsafe extern "C" fn pineappl_lumi_add( pdg_id_pairs .chunks(2) .zip(factors) - .map(|x| (x.0[0], x.0[1], x.1)) + .map(|x| (vec![x.0[0], x.0[1]], x.1)) .collect(), )); } @@ -1076,12 +1076,12 @@ pub unsafe extern "C" fn pineappl_lumi_entry( entry .iter() - .flat_map(|(id1, id2, _)| vec![id1, id2]) + .flat_map(|(pids, _)| pids) .zip(pdg_ids.iter_mut()) .for_each(|(from, to)| *to = *from); entry .iter() - .map(|(_, _, factor)| factor) + .map(|(_, factor)| factor) .zip(factors.iter_mut()) .for_each(|(from, to)| *to = *from); } diff --git a/pineappl_cli/src/plot.rs b/pineappl_cli/src/plot.rs index 87acbe1b..ea948508 100644 --- a/pineappl_cli/src/plot.rs +++ b/pineappl_cli/src/plot.rs @@ -104,15 +104,16 @@ fn map_format_parton(parton: i32) -> &'static str { } } +// TODO: generalize this function to n convolutions fn map_format_channel(channel: &Channel, has_pdf1: bool, has_pdf2: bool) -> String { channel .entry() .iter() - .map(|&(a, b, _)| { + .map(|(pids, _)| { format!( "{}{}", - if has_pdf1 { map_format_parton(a) } else { "" }, - if has_pdf2 { map_format_parton(b) } else { "" } + if has_pdf1 { map_format_parton(pids[0]) } else { "" }, + if has_pdf2 { map_format_parton(pids[1]) } else { "" } ) }) .join(" + ") diff --git a/pineappl_cli/src/read.rs b/pineappl_cli/src/read.rs index b08e4e15..33230c7d 100644 --- a/pineappl_cli/src/read.rs +++ b/pineappl_cli/src/read.rs @@ -124,8 +124,8 @@ impl Subcommand for Opts { row.add_cell(cell!(format!("{index}"))); - for (id1, id2, factor) in channel.entry() { - row.add_cell(cell!(format!("{factor} \u{d7} ({id1:2}, {id2:2})"))); + for (pids, factor) in channel.entry() { + row.add_cell(cell!(format!("{factor} \u{d7} ({})", pids.iter().map(ToString::to_string).join(", ")))); } } } else if self.group.ew || self.group.qcd { diff --git a/pineappl_cli/src/write.rs b/pineappl_cli/src/write.rs index 3de51c3b..522cf53f 100644 --- a/pineappl_cli/src/write.rs +++ b/pineappl_cli/src/write.rs @@ -494,9 +494,13 @@ impl Subcommand for Opts { for arg in &self.more_args.args { match arg { + // TODO: generalize to arbitrary convolutions OpsArg::Cc1(true) | OpsArg::Cc2(true) => { - let cc1 = matches!(arg, OpsArg::Cc1(true)); - let cc2 = matches!(arg, OpsArg::Cc2(true)); + let index = match arg { + OpsArg::Cc1(true) => 0, + OpsArg::Cc2(true) => 1, + _ => unreachable!(), + }; let pid_basis = grid.pid_basis(); @@ -505,29 +509,17 @@ impl Subcommand for Opts { channel .entry() .iter() - .map(|&(a, b, f)| { - let (ap, f1) = if cc1 { - pid_basis.charge_conjugate(a) - } else { - (a, 1.0) - }; - let (bp, f2) = if cc2 { - pid_basis.charge_conjugate(b) - } else { - (b, 1.0) - }; - (ap, bp, f * f1 * f2) + .map(|(pids, f)| { + let mut cc_pids = pids.clone(); + let (cc_pid, f1) = pid_basis.charge_conjugate(pids[index]); + cc_pids[index] = cc_pid; + (cc_pids, f * f1) }) .collect(), ); } - if cc1 { - grid.set_convolution(0, grid.convolutions()[0].charge_conjugate()); - } - if cc2 { - grid.set_convolution(1, grid.convolutions()[1].charge_conjugate()); - } + grid.set_convolution(index, grid.convolutions()[index].charge_conjugate()); } OpsArg::DedupChannels(ulps) => { grid.dedup_channels(*ulps);