Skip to content

Commit

Permalink
Migrate the remaining parts of plot to ConvFuns
Browse files Browse the repository at this point in the history
  • Loading branch information
@cschwan
cschwan committed Jun 27, 2024
1 parent c59f38e commit ddbf227
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 65 deletions.
115 changes: 78 additions & 37 deletions pineappl_cli/src/helpers.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,23 +97,6 @@ pub fn create_conv_funs_for_set(
Ok((set, conv_funs))
}

pub fn create_pdfset(pdfset: &str) -> Result<(PdfSet, Option<usize>)> {
let pdfset = pdfset.split_once('=').map_or(pdfset, |(name, _)| name);
let (pdfset, member) = pdfset
.rsplit_once('/')
.map_or((pdfset, None), |(set, member)| {
(set, Some(member.parse::<usize>().unwrap()))
});

Ok((
PdfSet::new(&pdfset.parse().map_or_else(
|_| pdfset.to_owned(),
|lhaid| lhapdf::lookup_pdf(lhaid).unwrap().0,
))?,
member,
))
}

pub fn read_grid(input: &Path) -> Result<Grid> {
Grid::read(File::open(input).context(format!("unable to open '{}'", input.display()))?)
.context(format!("unable to read '{}'", input.display()))
Expand Down Expand Up @@ -381,36 +364,94 @@ pub fn convolve_limits(grid: &Grid, bins: &[usize], mode: ConvoluteMode) -> Vec<

pub fn convolve_subgrid(
grid: &Grid,
lhapdf: &mut Pdf,
conv_funs: &mut [Pdf],
order: usize,
bin: usize,
lumi: usize,
cfg: &GlobalConfiguration,
) -> Array3<f64> {
// if the field 'Particle' is missing we assume it's a proton PDF
let pdf_pdg_id = lhapdf
.set()
.entry("Particle")
.map_or(Ok(2212), |string| string.parse::<i32>())
.unwrap();

if cfg.force_positive {
lhapdf.set_force_positive(1);
for fun in conv_funs.iter_mut() {
fun.set_force_positive(1);
}
}

let x_max = lhapdf.x_max();
let x_min = lhapdf.x_min();
let mut pdf = |id, x, q2| {
if !cfg.allow_extrapolation && (x < x_min || x > x_max) {
0.0
} else {
lhapdf.xfx_q2(id, x, q2)
match conv_funs {
[fun] => {
// there's only one convolution function from which we can use the strong coupling
assert_eq!(cfg.use_alphas_from, 0);

// if the field 'Particle' is missing we assume it's a proton PDF
let pdg_id = fun
.set()
.entry("Particle")
.map_or(Ok(2212), |string| string.parse::<i32>())
.unwrap();

let x_max = fun.x_max();
let x_min = fun.x_min();
let mut alphas = |q2| fun.alphas_q2(q2);
let mut fun = |id, x, q2| {
if !cfg.allow_extrapolation && (x < x_min || x > x_max) {
0.0
} else {
fun.xfx_q2(id, x, q2)
}
};

let mut cache = LumiCache::with_one(pdg_id, &mut fun, &mut alphas);

grid.convolve_subgrid(&mut cache, order, bin, lumi, 1.0, 1.0)
}
};
let mut alphas = |q2| lhapdf.alphas_q2(q2);
let mut cache = LumiCache::with_one(pdf_pdg_id, &mut pdf, &mut alphas);
[fun1, fun2] => {
let pdg_id1 = fun1
.set()
.entry("Particle")
.map_or(Ok(2212), |string| string.parse::<i32>())
.unwrap();

let pdg_id2 = fun2
.set()
.entry("Particle")
.map_or(Ok(2212), |string| string.parse::<i32>())
.unwrap();

let x_max1 = fun1.x_max();
let x_min1 = fun1.x_min();
let x_max2 = fun2.x_max();
let x_min2 = fun2.x_min();

let mut alphas = |q2| match cfg.use_alphas_from {
0 => fun1.alphas_q2(q2),
1 => fun2.alphas_q2(q2),
_ => panic!(
"expected `use_alphas_from` to be `0` or `1`, is {}",
cfg.use_alphas_from
),
};
let mut fun1 = |id, x, q2| {
if !cfg.allow_extrapolation && (x < x_min1 || x > x_max1) {
0.0
} else {
fun1.xfx_q2(id, x, q2)
}
};

let mut fun2 = |id, x, q2| {
if !cfg.allow_extrapolation && (x < x_min2 || x > x_max2) {
0.0
} else {
fun2.xfx_q2(id, x, q2)
}
};

let mut cache =
LumiCache::with_two(pdg_id1, &mut fun1, pdg_id2, &mut fun2, &mut alphas);

grid.convolve_subgrid(&mut cache, order, bin, lumi, 1.0, 1.0)
grid.convolve_subgrid(&mut cache, order, bin, lumi, 1.0, 1.0)
}
_ => unimplemented!(),
}
}

pub fn parse_integer_range(range: &str) -> Result<RangeInclusive<usize>> {
Expand Down
63 changes: 35 additions & 28 deletions pineappl_cli/src/plot.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,6 @@ use std::fmt::Write;
use std::num::NonZeroUsize;
use std::path::{Path, PathBuf};
use std::process::ExitCode;
use std::slice;
use std::thread;

/// Creates a matplotlib script plotting the contents of the grid.
Expand Down Expand Up @@ -529,19 +528,13 @@ impl Subcommand for Opts {
metadata = format_metadata(&vector),
);
} else {
let (pdfset1, pdfset2) = self
.conv_funs
.iter()
.map(|fun| {
assert_eq!(fun.lhapdf_names.len(), 1);
if let Some(member) = fun.members[0] {
format!("{}/{member}", fun.lhapdf_names[0])
} else {
fun.lhapdf_names[0].clone()
}
})
.collect_tuple()
.unwrap();
// TODO: enforce two arguments with clap
assert_eq!(self.conv_funs.len(), 2);

let (set1, mut conv_funs1) =
helpers::create_conv_funs_for_set(&self.conv_funs[0], self.conv_fun_uncert_from)?;
let (set2, mut conv_funs2) =
helpers::create_conv_funs_for_set(&self.conv_funs[1], self.conv_fun_uncert_from)?;
let (order, bin, channel) = self
.subgrid_pull
.iter()
Expand All @@ -552,17 +545,15 @@ impl Subcommand for Opts {
let cl = lhapdf::CL_1_SIGMA;
let grid = helpers::read_grid(&self.input)?;

let (set1, member1) = helpers::create_pdfset(&pdfset1)?;
let (set2, member2) = helpers::create_pdfset(&pdfset2)?;
let mut pdfset1 = set1.mk_pdfs()?;
let mut pdfset2 = set2.mk_pdfs()?;
let member1 = self.conv_funs[0].members[self.conv_fun_uncert_from];
let member2 = self.conv_funs[1].members[self.conv_fun_uncert_from];

let values1: Vec<_> = pdfset1
let values1: Vec<_> = conv_funs1
.par_iter_mut()
.map(|pdf| {
.map(|conv_funs| {
let values = helpers::convolve(
&grid,
slice::from_mut(pdf),
conv_funs,
&[],
&[bin],
&[],
Expand All @@ -574,12 +565,12 @@ impl Subcommand for Opts {
values[0]
})
.collect();
let values2: Vec<_> = pdfset2
let values2: Vec<_> = conv_funs2
.par_iter_mut()
.map(|pdf| {
.map(|conv_funs| {
let values = helpers::convolve(
&grid,
slice::from_mut(pdf),
conv_funs,
&[],
&[bin],
&[],
Expand Down Expand Up @@ -613,10 +604,26 @@ impl Subcommand for Opts {
unc1.hypot(unc2)
};

let res1 = helpers::convolve_subgrid(&grid, &mut pdfset1[0], order, bin, channel, cfg)
.sum_axis(Axis(0));
let res2 = helpers::convolve_subgrid(&grid, &mut pdfset2[0], order, bin, channel, cfg)
.sum_axis(Axis(0));
// TODO: if no member is given, the zeroth is used, but we should show the averaged
// result of all members instead
let res1 = helpers::convolve_subgrid(
&grid,
&mut conv_funs1[member1.unwrap_or(0)],
order,
bin,
channel,
cfg,
)
.sum_axis(Axis(0));
let res2 = helpers::convolve_subgrid(
&grid,
&mut conv_funs2[member2.unwrap_or(0)],
order,
bin,
channel,
cfg,
)
.sum_axis(Axis(0));

let subgrid = &grid.subgrids()[[order, bin, channel]];
//let q2 = subgrid.q2_grid();
Expand Down

0 comments on commit ddbf227

Please sign in to comment.