Skip to content

Commit

Permalink
Merge branch 'master' into eko-pids
Browse files Browse the repository at this point in the history 
* master:
  Add preliminary support for evolution basis
  Always optimize the luminosity function
  Use x1/x2 grid
  Use imu2 grid also in `LumiCache::xfx`
  Cache mu2/x1/x2 grids and use mu2 grid
  Add methods `clear` and `set_grids` to `LumiCache`
  Add new struct `LumiCache`
  Remove use of `RefCell`
  • Loading branch information
@alecandido
alecandido committed Oct 18, 2021
2 parents c422fb6 + 4e40743 commit 589c88d
Show file tree
Hide file tree
Showing 9 changed files with 485 additions and 41 deletions.
4 changes: 2 additions & 2 deletions pineappl/src/empty_subgrid.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ impl Subgrid for EmptySubgridV1 {
_: &[f64],
_: &[f64],
_: &[Mu2],
_: &dyn Fn(usize, usize, usize) -> f64,
_: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
0.0
}
Expand Down Expand Up @@ -65,7 +65,7 @@ mod tests {
#[test]
fn create_empty() {
let mut subgrid = EmptySubgridV1::default();
assert_eq!(subgrid.convolute(&[], &[], &[], &|_, _, _| 0.0), 0.0,);
assert_eq!(subgrid.convolute(&[], &[], &[], &mut |_, _, _| 0.0), 0.0,);
assert!(subgrid.is_empty());
subgrid.merge(&mut EmptySubgridV1::default().into(), false);
subgrid.scale(2.0);
Expand Down
145 changes: 123 additions & 22 deletions pineappl/src/grid.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,10 @@ use super::empty_subgrid::EmptySubgridV1;
use super::fk_table::FkTable;
use super::import_only_subgrid::ImportOnlySubgridV1;
use super::lagrange_subgrid::{LagrangeSparseSubgridV1, LagrangeSubgridV1, LagrangeSubgridV2};
use super::lumi::LumiEntry;
use super::lumi::{LumiCache, LumiEntry};
use super::lumi_entry;
use super::ntuple_subgrid::NtupleSubgridV1;
use super::pids;
use super::sparse_array3::SparseArray3;
use super::subgrid::{ExtraSubgridParams, Subgrid, SubgridEnum, SubgridParams};
use float_cmp::approx_eq;
Expand All @@ -19,7 +20,6 @@ use ndarray::{s, Array3, Array5, Dimension};
use rustc_hash::FxHashMap;
use serde::{Deserialize, Serialize};
use std::borrow::Cow;
use std::cell::RefCell;
use std::cmp::Ordering;
use std::collections::HashMap;
use std::convert::{TryFrom, TryInto};
Expand Down Expand Up @@ -316,6 +316,26 @@ impl Grid {
})
}

fn pdg_lumi(&self) -> Cow<[LumiEntry]> {
if let Some(key_values) = self.key_values() {
if let Some(lumi_id_types) = key_values.get("lumi_id_types") {
match lumi_id_types.as_str() {
"pdg_mc_ids" => {}
"evol" => {
return self
.lumi
.iter()
.map(|entry| LumiEntry::translate(entry, &pids::evol_to_pdg_mc_ids))
.collect();
}
_ => unimplemented!(),
}
}
}

Cow::Borrowed(self.lumi())
}

/// Performs a convolution of the contained subgrids with the given PDFs, `xfx1` for the first
/// parton and `xfx2` for the second parton, `alphas` for the evaluation of the strong
/// coupling. The parameters `order_mask` and `lumi_mask` can be used to selectively enable
Expand Down Expand Up @@ -347,9 +367,9 @@ impl Grid {
let bin_sizes = self.bin_info().normalizations();

// prepare pdf and alpha caches
let pdf_cache1 = RefCell::new(FxHashMap::default());
let pdf_cache2 = RefCell::new(FxHashMap::default());
let alphas_cache = RefCell::new(FxHashMap::default());
let mut pdf_cache1 = FxHashMap::default();
let mut pdf_cache2 = FxHashMap::default();
let mut alphas_cache = FxHashMap::default();
let mut last_xif = 0.0;

let (mut mu2_grid, mut x1_grid, mut x2_grid) = self
Expand All @@ -370,6 +390,7 @@ impl Grid {
.iter()
.map(|xi| xir_values.iter().position(|xir| xi.0 == *xir).unwrap())
.collect();
let self_lumi = self.pdg_lumi();

// iterate over the elements of `xi` and a corresponding index, but sorted using the
// factorisation value of `xi`
Expand All @@ -381,8 +402,8 @@ impl Grid {
{
// whenever the value `xif` changes we can clear the PDF cache
if xif != last_xif {
pdf_cache1.borrow_mut().clear();
pdf_cache2.borrow_mut().clear();
pdf_cache1.clear();
pdf_cache2.clear();
last_xif = xif;
}
// iterate subgrids
Expand All @@ -404,7 +425,7 @@ impl Grid {
None => continue,
};

let lumi_entry = &self.lumi[k];
let lumi_entry = &self_lumi[k];

let mut value = if subgrid.is_empty() {
0.0
Expand All @@ -417,20 +438,18 @@ impl Grid {

if mu2_grid_changed {
mu2_grid = new_mu2_grid;
alphas_cache.borrow_mut().clear();
alphas_cache.clear();
}

if mu2_grid_changed || (new_x1_grid != x1_grid) || (new_x2_grid != x2_grid) {
x1_grid = new_x1_grid;
x2_grid = new_x2_grid;
pdf_cache1.borrow_mut().clear();
pdf_cache2.borrow_mut().clear();
pdf_cache1.clear();
pdf_cache2.clear();
}

// pass convolute down, using caching
subgrid.convolute(&x1_grid, &x2_grid, &mu2_grid, &|ix1, ix2, imu2| {
let mut pdf_cache1 = pdf_cache1.borrow_mut();
let mut pdf_cache2 = pdf_cache2.borrow_mut();
subgrid.convolute(&x1_grid, &x2_grid, &mu2_grid, &mut |ix1, ix2, imu2| {
let x1 = x1_grid[ix1];
let x2 = x2_grid[ix2];
let mu2 = &mu2_grid[imu2];
Expand All @@ -454,7 +473,6 @@ impl Grid {
lumi += xfx1 * xfx2 * entry.2 / (x1 * x2);
}

let mut alphas_cache = alphas_cache.borrow_mut();
let alphas = alphas_cache
.entry(xir_values.len() * imu2 + xir_index)
.or_insert_with(|| alphas(xir * xir * mu2.ren));
Expand Down Expand Up @@ -483,6 +501,91 @@ impl Grid {
bins
}

/// TODO
///
/// Panics
///
/// TODO
pub fn convolute2(
&self,
lumi_cache: &mut LumiCache,
order_mask: &[bool],
bin_indices: &[usize],
lumi_mask: &[bool],
xi: &[(f64, f64)],
) -> Vec<f64> {
let bin_indices = if bin_indices.is_empty() {
(0..self.bin_limits.bins()).collect()
} else {
bin_indices.to_vec()
};
let mut bins: Vec<f64> = vec![0.0; bin_indices.len() * xi.len()];
let normalizations = self.bin_info().normalizations();
let self_lumi = self.pdg_lumi();

for (xi_index, &(xir, xif)) in xi.iter().enumerate() {
for ((ord, bin, lumi), subgrid) in self.subgrids.indexed_iter() {
let order = &self.orders[ord];

if ((order.logxir > 0) && (xir == 1.0)) || ((order.logxif > 0) && (xif == 1.0)) {
continue;
}

if (!order_mask.is_empty() && !order_mask[ord])
|| (!lumi_mask.is_empty() && !lumi_mask[lumi])
{
continue;
}

let bin_index = match bin_indices.iter().position(|&index| index == bin) {
Some(i) => i,
None => continue,
};

if subgrid.is_empty() {
continue;
}

let lumi_entry = &self_lumi[lumi];
let mu2_grid = subgrid.mu2_grid();
let x1_grid = subgrid.x1_grid();
let x2_grid = subgrid.x2_grid();

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

let mut value =
subgrid.convolute(&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 lumi_entry.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);
}

let alphas = lumi_cache.alphas(imu2);

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

if order.logxir > 0 {
value *= (xir * xir).ln().powi(order.logxir.try_into().unwrap());
}

if order.logxif > 0 {
value *= (xif * xif).ln().powi(order.logxif.try_into().unwrap());
}

bins[xi_index + xi.len() * bin_index] += value / normalizations[bin];
}
}

bins
}

/// Convolutes a single subgrid `(order, bin, lumi)` with the PDFs strong coupling given by
/// `xfx1`, `xfx2` and `alphas`. The convolution result is fully differentially, such that the
/// axes of the result correspond to the values given by the subgrid `q2`, `x1` and `x2` grid
Expand All @@ -504,9 +607,9 @@ impl Grid {
) -> Array3<f64> {
let normalization = self.bin_info().normalizations()[bin];

let pdf_cache1 = RefCell::new(FxHashMap::default());
let pdf_cache2 = RefCell::new(FxHashMap::default());
let alphas_cache = RefCell::new(FxHashMap::default());
let mut pdf_cache1 = FxHashMap::default();
let mut pdf_cache2 = FxHashMap::default();
let mut alphas_cache = FxHashMap::default();

let subgrid = &self.subgrids[[order, bin, lumi]];
let order = &self.orders[order];
Expand All @@ -527,8 +630,6 @@ impl Grid {
let mut array = Array3::zeros((mu2_grid.len(), x1_grid.len(), x2_grid.len()));

for ((imu2, ix1, ix2), value) in subgrid.iter() {
let mut pdf_cache1 = pdf_cache1.borrow_mut();
let mut pdf_cache2 = pdf_cache2.borrow_mut();
let x1 = x1_grid[ix1];
let x2 = x2_grid[ix2];
let mu2 = &mu2_grid[imu2];
Expand All @@ -552,7 +653,6 @@ impl Grid {
lumi += xfx1 * xfx2 * entry.2 / (x1 * x2);
}

let mut alphas_cache = alphas_cache.borrow_mut();
let alphas = alphas_cache
.entry(imu2)
.or_insert_with(|| alphas(xir * xir * mu2.ren));
Expand Down Expand Up @@ -937,9 +1037,10 @@ impl Grid {
.map_or(true, |map| map["initial_state_1"] == map["initial_state_2"])
{
self.symmetrize_lumi();
self.optimize_lumi();
}

self.optimize_lumi();

for subgrid in self.subgrids.iter_mut() {
if subgrid.is_empty() {
*subgrid = EmptySubgridV1::default().into();
Expand Down
10 changes: 6 additions & 4 deletions pineappl/src/import_only_subgrid.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ impl Subgrid for ImportOnlySubgridV1 {
_: &[f64],
_: &[f64],
_: &[Mu2],
lumi: &dyn Fn(usize, usize, usize) -> f64,
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.array
.indexed_iter()
Expand Down Expand Up @@ -251,7 +251,7 @@ impl Subgrid for ImportOnlySubgridV2 {
_: &[f64],
_: &[f64],
_: &[Mu2],
lumi: &dyn Fn(usize, usize, usize) -> f64,
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.array
.indexed_iter()
Expand Down Expand Up @@ -397,7 +397,8 @@ mod tests {
assert_eq!(grid1.iter().nth(3), Some(((0, 7, 1), &8.0)));

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

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

Expand Down Expand Up @@ -482,7 +483,8 @@ mod tests {
assert_eq!(grid1.iter().nth(3), Some(((0, 7, 1), &8.0)));

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

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

Expand Down
24 changes: 14 additions & 10 deletions pineappl/src/lagrange_subgrid.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ impl Subgrid for LagrangeSubgridV1 {
x1: &[f64],
x2: &[f64],
_: &[Mu2],
lumi: &dyn Fn(usize, usize, usize) -> f64,
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.grid.as_ref().map_or(0.0, |grid| {
grid.indexed_iter()
Expand Down Expand Up @@ -450,7 +450,7 @@ impl Subgrid for LagrangeSubgridV2 {
x1: &[f64],
x2: &[f64],
_: &[Mu2],
lumi: &dyn Fn(usize, usize, usize) -> f64,
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.grid.as_ref().map_or(0.0, |grid| {
grid.indexed_iter()
Expand Down Expand Up @@ -745,7 +745,7 @@ impl Subgrid for LagrangeSparseSubgridV1 {
x1: &[f64],
x2: &[f64],
_: &[Mu2],
lumi: &dyn Fn(usize, usize, usize) -> f64,
lumi: &mut dyn FnMut(usize, usize, usize) -> f64,
) -> f64 {
self.array
.indexed_iter()
Expand Down Expand Up @@ -958,7 +958,8 @@ mod tests {
let x2 = grid.x2_grid();
let mu2 = grid.mu2_grid();

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

let mut test = 0.0;

Expand Down Expand Up @@ -1008,13 +1009,14 @@ mod tests {
let x2 = grid1.x2_grid().into_owned();
let mu2 = grid1.mu2_grid().into_owned();

let reference = grid1.convolute(&x1, &x2, &mu2, &|ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));
let reference =
grid1.convolute(&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.convolute(&x1, &x2, &mu2, &|ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));
let merged = grid2.convolute(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));

assert!(approx_eq!(f64, reference, merged, ulps = 8));

Expand Down Expand Up @@ -1045,7 +1047,7 @@ mod tests {

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

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

assert!(approx_eq!(f64, 2.0 * reference, merged, ulps = 8));
}
Expand Down Expand Up @@ -1100,7 +1102,7 @@ mod tests {
let x2 = grid.x2_grid();
let mu2 = grid.mu2_grid();

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

assert_eq!(result, 0.0);
}
Expand Down Expand Up @@ -1158,8 +1160,10 @@ mod tests {
let sparse = LagrangeSparseSubgridV1::from(&dense);
assert!(!sparse.is_empty());

let reference = dense.convolute(&x1, &x2, &mu2, &|ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));
let converted = sparse.convolute(&x1, &x2, &mu2, &|ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));
let reference =
dense.convolute(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));
let converted =
sparse.convolute(&x1, &x2, &mu2, &mut |ix1, ix2, _| 1.0 / (x1[ix1] * x2[ix2]));

assert!(approx_eq!(f64, reference, converted, ulps = 8));
}
Expand Down
Loading

0 comments on commit 589c88d

Please sign in to comment.