Skip to content

Commit

Permalink
Remove Subgrid::convolve
Browse files Browse the repository at this point in the history
  • Loading branch information
t7phy committed Aug 16, 2024
1 parent df2cb57 commit 5d71c36
Show file tree
Hide file tree
Showing 5 changed files with 18 additions and 216 deletions.
11 changes: 0 additions & 11 deletions pineappl/src/empty_subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -10,16 +10,6 @@ use std::iter;
pub struct EmptySubgridV1;

impl Subgrid for EmptySubgridV1 {
fn convolve(
&self,
_: &[f64],
_: &[f64],
_: &[Mu2],
_: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
0.0
}

fn fill(&mut self, _: &[f64], _: f64) {
panic!("EmptySubgridV1 doesn't support the fill operation");
}
Expand Down Expand Up @@ -81,7 +71,6 @@ mod tests {
#[test]
fn create_empty() {
let mut subgrid = EmptySubgridV1;
assert_eq!(subgrid.convolve(&[], &[], &[], &mut |_, _, _| 0.0), 0.0,);
assert!(subgrid.is_empty());
subgrid.merge(&mut EmptySubgridV1.into(), false);
subgrid.scale(2.0);
Expand Down
34 changes: 18 additions & 16 deletions pineappl/src/grid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -309,24 +309,26 @@ impl Grid {

lumi_cache.set_grids(&mu2_grid, &x1_grid, &x2_grid, xir, xif);

let mut value =
subgrid.convolve(&x1_grid, &x2_grid, &mu2_grid, &mut |ix1, ix2, imu2| {
let x1 = x1_grid[ix1];
let x2 = x2_grid[ix2];
let mut lumi = 0.0;

for entry in channel.entry() {
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 mut value = 0.0;

for ((imu2, ix1, ix2), v) in subgrid.indexed_iter() {
let x1 = x1_grid[ix1];
let x2 = x2_grid[ix2];
let mut lumi = 0.0;

for entry in channel.entry() {
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);
let alphas = lumi_cache.alphas(imu2);

lumi *= alphas.powi(order.alphas.try_into().unwrap());
lumi
});
lumi *= alphas.powi(order.alphas.try_into().unwrap());

value += lumi * v;
}

if order.logxir > 0 {
value *= (xir * xir).ln().powi(order.logxir.try_into().unwrap());
Expand Down
155 changes: 0 additions & 155 deletions pineappl/src/lagrange_subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -162,33 +162,6 @@ impl LagrangeSubgridV2 {
}

impl Subgrid for LagrangeSubgridV2 {
fn convolve(
&self,
x1: &[f64],
x2: &[f64],
_: &[Mu2],
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.grid.as_ref().map_or(0.0, |grid| {
grid.indexed_iter()
.map(|((imu2, ix1, ix2), &sigma)| {
if sigma == 0.0 {
0.0
} else {
let mut value = sigma * lumi(ix1, ix2, imu2 + self.itaumin);
if self.reweight1 {
value *= weightfun(x1[ix1]);
}
if self.reweight2 {
value *= weightfun(x2[ix2]);
}
value
}
})
.sum()
})
}

fn fill(&mut self, ntuple: &[f64], weight: f64) {
if weight == 0.0 {
return;
Expand Down Expand Up @@ -458,117 +431,6 @@ impl Subgrid for LagrangeSubgridV2 {
#[cfg(test)]
mod tests {
use super::*;
use float_cmp::assert_approx_eq;

fn test_q2_slice_methods<G: Subgrid>(mut grid: G) -> G {
grid.fill(&[0.1, 0.2, 90.0_f64.powi(2)], 1.0);
grid.fill(&[0.9, 0.1, 90.0_f64.powi(2)], 1.0);
grid.fill(&[0.009, 0.01, 90.0_f64.powi(2)], 1.0);
grid.fill(&[0.009, 0.5, 90.0_f64.powi(2)], 1.0);

// the grid must not be empty
assert!(!grid.is_empty());

let x1 = grid.x1_grid();
let x2 = grid.x2_grid();
let mu2 = grid.mu2_grid();

let reference = grid.convolve(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));

let mut test = 0.0;

// check `reference` against manually calculated result from q2 slices
for ((_, ix1, ix2), value) in grid.indexed_iter() {
test += value / (x1[ix1] * x2[ix2]);
}

assert_approx_eq!(f64, test, reference, ulps = 8);

grid
}

fn test_merge_method<G: Subgrid>(mut grid1: G, mut grid2: G, mut grid3: G)
where
SubgridEnum: From<G>,
{
grid1.fill(&[0.1, 0.2, 90.0_f64.powi(2)], 1.0);
grid1.fill(&[0.9, 0.1, 90.0_f64.powi(2)], 1.0);
grid1.fill(&[0.009, 0.01, 90.0_f64.powi(2)], 1.0);
grid1.fill(&[0.009, 0.5, 90.0_f64.powi(2)], 1.0);

assert!(!grid1.is_empty());
assert!(grid2.is_empty());

let x1 = grid1.x1_grid().into_owned();
let x2 = grid1.x2_grid().into_owned();
let mu2 = grid1.mu2_grid().into_owned();

let reference =
grid1.convolve(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));

// merge filled grid into empty one
grid2.merge(&mut grid1.into(), false);
assert!(!grid2.is_empty());

let merged = grid2.convolve(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));

assert_approx_eq!(f64, reference, merged, ulps = 8);

grid3.fill(&[0.1, 0.2, 90.0_f64.powi(2)], 1.0);
grid3.fill(&[0.9, 0.1, 90.0_f64.powi(2)], 1.0);
grid3.fill(&[0.009, 0.01, 90.0_f64.powi(2)], 1.0);
grid3.fill(&[0.009, 0.5, 90.0_f64.powi(2)], 1.0);

grid2.merge(&mut grid3.into(), false);

let merged = grid2.convolve(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));

assert_approx_eq!(f64, 2.0 * reference, merged, ulps = 8);
}

fn test_empty_subgrid<G: Subgrid>(mut grid: G) {
// this following events should be skipped

// q2 is too large
grid.fill(&[0.5, 0.5, 2e+8], 1.0);
// q2 is too small
grid.fill(&[0.5, 0.5, 5e+1], 1.0);
// x1 is too large
grid.fill(&[1.1, 0.5, 1e+3], 1.0);
// x1 is too small
grid.fill(&[0.5, 1e-7, 1e+3], 1.0);
// x1 is too large
grid.fill(&[0.5, 1.1, 1e+3], 1.0);
// x1 is too small
grid.fill(&[1e-7, 0.5, 1e+3], 1.0);

let x1 = grid.x1_grid();
let x2 = grid.x2_grid();
let mu2 = grid.mu2_grid();

let result = grid.convolve(&x1, &x2, &mu2, &mut |_, _, _| 1.0);

assert_eq!(result, 0.0);
}

#[test]
fn q2_slice_v2() {
let subgrid = test_q2_slice_methods(LagrangeSubgridV2::new(
&SubgridParams::default(),
&ExtraSubgridParams::default(),
));

assert_eq!(
subgrid.stats(),
Stats {
total: 10000,
allocated: 10000,
zeros: 256,
overhead: 0,
bytes_per_value: 8
}
);
}

#[test]
fn fill_zero_v2() {
Expand All @@ -580,21 +442,4 @@ mod tests {
assert!(subgrid.is_empty());
assert_eq!(subgrid.indexed_iter().count(), 0);
}

#[test]
fn merge_dense_v2() {
test_merge_method(
LagrangeSubgridV2::new(&SubgridParams::default(), &ExtraSubgridParams::default()),
LagrangeSubgridV2::new(&SubgridParams::default(), &ExtraSubgridParams::default()),
LagrangeSubgridV2::new(&SubgridParams::default(), &ExtraSubgridParams::default()),
);
}

#[test]
fn empty_v2() {
test_empty_subgrid(LagrangeSubgridV2::new(
&SubgridParams::default(),
&ExtraSubgridParams::default(),
));
}
}
24 changes: 0 additions & 24 deletions pineappl/src/packed_subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -39,19 +39,6 @@ impl PackedQ1X2SubgridV1 {
}

impl Subgrid for PackedQ1X2SubgridV1 {
fn convolve(
&self,
_: &[f64],
_: &[f64],
_: &[Mu2],
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.array
.indexed_iter()
.map(|([imu2, ix1, ix2], sigma)| sigma * lumi(ix1, ix2, imu2))
.sum()
}

fn fill(&mut self, _: &[f64], _: f64) {
panic!("PackedQ1X2SubgridV1 doesn't support the fill operation");
}
Expand Down Expand Up @@ -312,12 +299,6 @@ mod tests {
assert_eq!(grid1.indexed_iter().nth(2), Some(((0, 4, 3), 4.0)));
assert_eq!(grid1.indexed_iter().nth(3), Some(((0, 7, 1), 8.0)));

// symmetric luminosity function
let lumi =
&mut (|ix1, ix2, _| x[ix1] * x[ix2]) as &mut dyn FnMut(usize, usize, usize) -> f64;

assert_eq!(grid1.convolve(&x, &x, &mu2, lumi), 0.228515625);

// create grid with transposed entries, but different q2
let mut grid2: SubgridEnum = PackedQ1X2SubgridV1::new(
PackedArray::new([1, 10, 10]),
Expand All @@ -338,7 +319,6 @@ mod tests {
} else {
unreachable!();
}
assert_eq!(grid2.convolve(&x, &x, &mu2, lumi), 0.228515625);

assert_eq!(grid2.indexed_iter().next(), Some(((0, 1, 7), 8.0)));
assert_eq!(grid2.indexed_iter().nth(1), Some(((0, 2, 1), 1.0)));
Expand All @@ -347,8 +327,6 @@ mod tests {

grid1.merge(&mut grid2, false);

assert_eq!(grid1.convolve(&x, &x, &mu2, lumi), 2.0 * 0.228515625);

let mut grid1 = {
let mut g = grid1.clone_empty();
g.merge(&mut grid1, false);
Expand All @@ -358,10 +336,8 @@ mod tests {
// the luminosity function is symmetric, so after symmetrization the result must be
// unchanged
grid1.symmetrize();
assert_eq!(grid1.convolve(&x, &x, &mu2, lumi), 2.0 * 0.228515625);

grid1.scale(2.0);
assert_eq!(grid1.convolve(&x, &x, &mu2, lumi), 4.0 * 0.228515625);

assert_eq!(
grid1.stats(),
Expand Down
10 changes: 0 additions & 10 deletions pineappl/src/subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -68,16 +68,6 @@ pub trait Subgrid {
/// return an empty slice.
fn x2_grid(&self) -> Cow<[f64]>;

/// Convolute the subgrid with a luminosity function, which takes indices as arguments that
/// correspond to the entries given in the slices `x1`, `x2` and `mu2`.
fn convolve(
&self,
x1: &[f64],
x2: &[f64],
mu2: &[Mu2],
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64;

/// Fills the subgrid with `weight` for the parton momentum fractions `x1` and `x2`, and the
/// scale `q2`. Filling is currently only support where both renormalization and factorization
/// scale have the same value.
Expand Down

0 comments on commit 5d71c36

Please sign in to comment.