Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize the APPLgrid exporter #285

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Generalize the APPLgrid exporter #285

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ca8dd1c
Generalize the APPLgrid exporter
janw20 Apr 30, 2024
503c12c
Increase `ulps` and remove `epsilon` in `approx_eq!`
janw20 May 3, 2024
2159c87
Remove leftover comment
janw20 May 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
94 changes: 61 additions & 33 deletions pineappl_cli/src/export/applgrid.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
use anyhow::{anyhow, bail, Result};
use cxx::{let_cxx_string, UniquePtr};
use float_cmp::approx_eq;
use itertools::Itertools;
use ndarray::{s, Axis};
use pineappl::grid::{Grid, Order};
use pineappl::subgrid::{Mu2, Subgrid, SubgridParams};
Expand All @@ -18,32 +19,61 @@ fn reconstruct_subgrid_params(grid: &Grid, order: usize, bin: usize) -> Result<S
.subgrids()
.slice(s![order, bin, ..])
.iter()
.map(|subgrid| {
.flat_map(|subgrid| {
subgrid
.mu2_grid()
.iter()
.into_iter()
.map(|&Mu2 { ren, fac }| {
if !approx_eq!(f64, ren, fac, ulps = 128) {
bail!("subgrid has mur2 != muf2, which APPLgrid does not support");
}

Ok(fac)
})
.collect::<Result<Vec<_>>>()
.collect_vec()
})
.collect::<Result<_>>()?;
let mut mu2_grid: Vec<_> = mu2_grid.into_iter().flatten().collect();
mu2_grid.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = 128));
let mu2_grid = mu2_grid.as_slice();

if let &[fac] = mu2_grid {
result.set_q2_bins(1);
result.set_q2_max(fac);
result.set_q2_min(fac);
.collect::<Result<Vec<_>>>()?
.into_iter()
.sorted_unstable_by(|a, b| a.total_cmp(b))
.dedup_by(|&a, &b| approx_eq!(f64, a, b, ulps = 128))
.collect_vec();

let x_grid = grid
.subgrids()
.slice(s![order, bin, ..])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This probably doesn't work quite as expected. APPLgrid stores different channels in a single appl_igrid, which has common parameter values NQ2, NQ2min, NQ2max, NQorder, Nx, xmin, xmax, xorder. This means that:

  • if several channels have different SubgridParam/ExtraSubgridParam values we should error out, unless this is a DIS grid in which we ignore the second dimension. In PineAPPL the second dimension can be either x1 or x2
  • the flat_map should therefore probably be a simple map and we must make sure that for all channels the entries are the same
  • the same applies also for mu_grid

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The OPTIMIZE_SUBGRID_TYPE optimization can change the ends of the x and mu2 grid points, right? So I would check if the grid points are equal modulo some endpoints, and then use the union as the x points that will be passed to APPLgrid?

Copy link
Contributor

@cschwan cschwan May 6, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would address one point of it, yes. It'll get interesting if you've got a single channel with different minimums/maximums.

.iter()
.flat_map(|subgrid| {
[
subgrid.x1_grid().into_owned(),
subgrid.x2_grid().into_owned(),
]
.into_iter()
})
.flatten()
.sorted_unstable_by(|a, b| a.total_cmp(b))
.dedup_by(|&a, &b| approx_eq!(f64, a, b, ulps = 128))
.collect_vec();

result.set_q2_bins(mu2_grid.len());
result.set_q2_min(*mu2_grid.first().unwrap());
result.set_q2_max(*mu2_grid.last().unwrap());
if mu2_grid.len() == 1 {
result.set_q2_order(0);
} else {
// TODO: reconstruct the interpolation order. for now leave it as the default
// result.set_q2_order(...);
}

result.set_x_bins(x_grid.len());
result.set_x_min(*x_grid.first().unwrap());
result.set_x_max(*x_grid.last().unwrap());
if x_grid.len() == 1 {
result.set_x_order(0);
} else {
// TODO: reconstruct the interpolation order. for now leave it as the default
// result.set_x_order(...);
}

// TODO: implement the general case
Ok(result)
}

Expand All @@ -55,15 +85,9 @@ pub fn convert_into_applgrid(
let bin_info = grid.bin_info();
let dim = bin_info.dimensions();

if dim > 1 {
bail!(
"grid has {} dimensions, but APPLgrid only supports one-dimensional distributions",
dim
);
}

if bin_info.slices().len() != 1 {
bail!("grid has non-consecutive bin limits, which APPLgrid does not support");
let integer_bin_limits = dim > 1 || bin_info.slices().len() > 1;
if integer_bin_limits {
println!("Info: APPLgrid only supports one-dimensional consecutive bin limits, so the bin limits will be set to [0, {}] in steps of 1, corresponding to the ordering given by `pineappl read --bins <GRID>`.", bin_info.bins());
}

let lumis = grid.lumi().len();
Expand Down Expand Up @@ -100,11 +124,15 @@ pub fn convert_into_applgrid(
ffi::make_lumi_pdf(id, &combinations).into_raw();

let limits = &bin_info.limits();
let limits: Vec<_> = limits
.iter()
.map(|vec| vec[0].0)
.chain(limits.last().map(|vec| vec[0].1))
.collect();
let limits = if integer_bin_limits {
(0..=limits.len()).map(|x| x as f64).collect::<Vec<_>>()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here there's very likely a normalization factor missing. If you change the bin sizes you must correct the cross sections.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, thanks. Does APPLgrid also divide by the bin widths in the end?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

AFAICR yes!

} else {
limits
.iter()
.map(|vec| vec[0].0)
.chain(limits.last().map(|vec| vec[0].1))
.collect_vec()
};

let order_mask = Order::create_mask(grid.orders(), 3, 0, false);
let orders_with_mask: Vec<_> = grid
Expand Down Expand Up @@ -209,12 +237,12 @@ pub fn convert_into_applgrid(
.map(|&x1| {
appl_x1
.iter()
.position(|&x| approx_eq!(f64, x, x1, ulps = 128))
.position(|&x| approx_eq!(f64, x, x1, ulps = 512))
.map_or_else(
|| {
Err(anyhow!(
"momentum fraction x1 = {} not found in APPLgrid",
x1
"momentum fraction x1 = {} not found in APPLgrid, the closest is {}",
x1, appl_x1.iter().min_by(|&a, &b| (a - x1).abs().total_cmp(&(b - x1).abs())).unwrap()
))
},
|idx| Ok(idx.try_into().unwrap()),
Expand All @@ -226,12 +254,12 @@ pub fn convert_into_applgrid(
.map(|&x2| {
appl_x2
.iter()
.position(|&x| approx_eq!(f64, x, x2, ulps = 128))
.position(|&x| approx_eq!(f64, x, x2, ulps = 512))
.map_or_else(
|| {
Err(anyhow!(
"momentum fraction x2 = {} not found in APPLgrid",
x2
"momentum fraction x2 = {} not found in APPLgrid, the closest is {}",
x2, appl_x2.iter().min_by(|&a, &b| (a - x2).abs().total_cmp(&(b - x2).abs())).unwrap()
))
},
|idx| Ok(idx.try_into().unwrap()),
Expand Down
Loading