Skip to content

Commit

Permalink
Start getting rid of Convolution::None
Browse files Browse the repository at this point in the history
  • Loading branch information
cschwan committed Oct 10, 2024
1 parent 0632b03 commit ae4b9a1
Show file tree
Hide file tree
Showing 8 changed files with 157 additions and 131 deletions.
5 changes: 0 additions & 5 deletions pineappl/src/convolutions.rs
Original file line number Diff line number Diff line change
Expand Up @@ -52,11 +52,6 @@ impl<'a> ConvolutionCache<'a> {
}

pub(crate) fn setup(&mut self, grid: &Grid, xi: &[(f64, f64, f64)]) -> Result<(), ()> {
let convolutions = grid.convolutions();

// TODO: the following code only works with exactly two convolutions
assert_eq!(convolutions.len(), 2);

self.perm = grid
.convolutions()
.iter()
Expand Down
2 changes: 1 addition & 1 deletion pineappl/src/empty_subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ mod tests {
#[should_panic(expected = "EmptySubgridV1 doesn't support the fill operation")]
fn fill() {
let mut subgrid = EmptySubgridV1;
subgrid.fill(&v0::default_interps(), &[0.0; 3], 0.0);
subgrid.fill(&v0::default_interps(2), &[0.0; 3], 0.0);
}

#[test]
Expand Down
63 changes: 33 additions & 30 deletions pineappl/src/evolution.rs
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ use float_cmp::approx_eq;
use itertools::izip;
use itertools::Itertools;
use ndarray::linalg;
use ndarray::{s, Array1, Array2, Array3, ArrayD, ArrayView1, ArrayView4, Axis, Ix2};
use ndarray::{s, Array1, Array2, Array3, ArrayD, ArrayView1, ArrayView4, Axis, Ix1, Ix2};
use std::iter;

/// Number of ULPS used to de-duplicate grid values in [`Grid::evolve_info`].
Expand Down Expand Up @@ -980,8 +980,6 @@ pub(crate) fn evolve_slice_with_many(
let dim: Vec<_> = infos.iter().map(|info| info.x0.len()).collect();

for subgrids_oc in grid.subgrids().axis_iter(Axis(1)) {
assert_eq!(infos[0].x0.len(), infos[1].x0.len());

let mut tables = vec![ArrayD::zeros(dim.clone()); channels0.len()];

for (subgrids_o, channel1) in subgrids_oc.axis_iter(Axis(1)).zip(grid.channels()) {
Expand Down Expand Up @@ -1020,19 +1018,13 @@ pub(crate) fn evolve_slice_with_many(
}
}

let mut tmp = Array2::zeros((last_x1[0].len(), infos[1].x0.len()));

for (pids1, factor) in channel1
.entry()
.iter()
.map(|(pids1, factor)| ([pids1[0], pids1[1]], factor))
{
for (pids1, factor) in channel1.entry() {
for (fk_table, ops) in
channels0
.iter()
.zip(tables.iter_mut())
.filter_map(|(pids0, fk_table)| {
izip!(pids0, &pids1, &pids01, &eko_slices)
izip!(pids0, pids1, &pids01, &eko_slices)
.map(|(&pid0, &pid1, pids, slices)| {
pids.iter().zip(slices).find_map(|(&(p0, p1), op)| {
((p0 == pid0) && (p1 == pid1)).then_some(op)
Expand All @@ -1043,20 +1035,22 @@ pub(crate) fn evolve_slice_with_many(
.map(|ops| (fk_table, ops))
})
{
general_tensor_mul(*factor, &array, &ops, &mut tmp, fk_table);
general_tensor_mul(*factor, &array, &ops, fk_table);
}
}
}

// TODO: generalize this for arbitrary scales and x values
let mut node_values = vec![NodeValues::UseThese(vec![infos[0].fac0])];

for info in infos {
node_values.push(NodeValues::UseThese(info.x0.clone()));
}

sub_fk_tables.extend(tables.into_iter().map(|table| {
PackedQ1X2SubgridV1::new(
PackedArray::from(table.insert_axis(Axis(0)).view()),
// TODO: generalize this for arbitrary scales and x values
vec![
NodeValues::UseThese(vec![infos[0].fac0]),
NodeValues::UseThese(infos[0].x0.clone()),
NodeValues::UseThese(infos[1].x0.clone()),
],
node_values.clone(),
)
.into()
}));
Expand All @@ -1067,8 +1061,8 @@ pub(crate) fn evolve_slice_with_many(
.into_shape((1, grid.bin_info().bins(), channels0.len()))
.unwrap(),
channels0
.iter()
.map(|c| channel![c[0], c[1], 1.0])
.into_iter()
.map(|c| Channel::new(vec![(c, 1.0)]))
.collect(),
))
}
Expand All @@ -1077,16 +1071,25 @@ fn general_tensor_mul(
factor: f64,
array: &ArrayD<f64>,
ops: &[&Array2<f64>],
tmp: &mut Array2<f64>,
fk_table: &mut ArrayD<f64>,
) {
// TODO: generalize this to n dimensions
assert_eq!(array.shape().len(), 2);
let array = array.view().into_dimensionality::<Ix2>().unwrap();
let mut fk_table = fk_table.view_mut().into_dimensionality::<Ix2>().unwrap();

// tmp = array * ops[1]^T
linalg::general_mat_mul(1.0, &array, &ops[1].t(), 0.0, tmp);
// fk_table += factor * ops[0] * tmp
linalg::general_mat_mul(factor, ops[0], &tmp, 1.0, &mut fk_table);
match array.shape().len() {
1 => {
let array = array.view().into_dimensionality::<Ix1>().unwrap();
let mut fk_table = fk_table.view_mut().into_dimensionality::<Ix1>().unwrap();
fk_table.scaled_add(factor, &ops[0].dot(&array));
}
2 => {
let array = array.view().into_dimensionality::<Ix2>().unwrap();
let mut fk_table = fk_table.view_mut().into_dimensionality::<Ix2>().unwrap();

let mut tmp = Array2::zeros((array.shape()[0], ops[1].shape()[0]));
// 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, &mut fk_table);
}
// TODO: generalize this to n dimensions
_ => unimplemented!(),
}
}
39 changes: 20 additions & 19 deletions pineappl/src/grid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -1275,15 +1275,16 @@ impl Grid {
operator.dim(),
);

let view = operator.view();

let (subgrids, channels) = if self.convolutions()[0] != Convolution::None
&& self.convolutions()[1] != Convolution::None
{
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, alphas_table)
}?;
let views = vec![operator.view(); self.convolutions().len()];
let infos = vec![info.clone(); self.convolutions().len()];
let (subgrids, channels) = evolution::evolve_slice_with_many(
self,
&views,
&infos,
order_mask,
(xi.0, xi.1, 1.0),
alphas_table,
)?;

let rhs = Self {
subgrids,
Expand Down Expand Up @@ -1734,7 +1735,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.0, 1.0],
vec![Convolution::UnpolPDF(2212), Convolution::UnpolPDF(2212)],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand All @@ -1755,7 +1756,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.0, 0.25, 0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand All @@ -1778,7 +1779,7 @@ mod tests {
vec![Order::new(1, 2, 0, 0, 0), Order::new(1, 2, 0, 1, 0)],
vec![0.0, 0.25, 0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand Down Expand Up @@ -1806,7 +1807,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.0, 0.25, 0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand All @@ -1832,7 +1833,7 @@ mod tests {
],
vec![0.0, 0.25, 0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand Down Expand Up @@ -1865,7 +1866,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.0, 0.25, 0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand All @@ -1884,7 +1885,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.0, 0.25, 0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand Down Expand Up @@ -1914,7 +1915,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.0, 0.25, 0.5],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand All @@ -1937,7 +1938,7 @@ mod tests {
vec![Order::new(0, 2, 0, 0, 0)],
vec![0.5, 0.75, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand Down Expand Up @@ -1972,7 +1973,7 @@ mod tests {
}],
vec![0.0, 1.0],
vec![Convolution::UnpolPDF(2212); 2],
v0::default_interps(),
v0::default_interps(2),
vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2],
Scales {
ren: ScaleFuncForm::Scale(0),
Expand Down
6 changes: 3 additions & 3 deletions pineappl/src/lagrange_subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ mod tests {

#[test]
fn fill_zero() {
let interps = v0::default_interps();
let interps = v0::default_interps(2);
let mut subgrid = LagrangeSubgridV2::new(&interps);

subgrid.fill(&interps, &[1000.0, 0.5, 0.5], 0.0);
Expand All @@ -148,7 +148,7 @@ mod tests {

#[test]
fn fill_outside_range() {
let interps = v0::default_interps();
let interps = v0::default_interps(2);
let mut subgrid = LagrangeSubgridV2::new(&interps);

subgrid.fill(&interps, &[1000.0, 1e-10, 0.5], 0.0);
Expand All @@ -169,7 +169,7 @@ mod tests {

#[test]
fn fill() {
let interps = v0::default_interps();
let interps = v0::default_interps(2);
let mut subgrid = LagrangeSubgridV2::new(&interps);

subgrid.fill(&interps, &[1000.0, 0.5, 0.5], 1.0);
Expand Down
2 changes: 1 addition & 1 deletion pineappl/src/packed_subgrid.rs
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ mod tests {
PackedArray::new(vec![0, 0, 0]),
vec![NodeValues::UseThese(Vec::new()); 3],
);
subgrid.fill(&v0::default_interps(), &[0.0; 3], 0.0);
subgrid.fill(&v0::default_interps(2), &[0.0; 3], 0.0);
}

#[test]
Expand Down
Loading

0 comments on commit ae4b9a1

Please sign in to comment.