diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 2a089d1e..192888b4 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -35,6 +35,7 @@ jobs: # `-C link-dead-code` is needed to prevent 'warning: XX functions have mismatched data' warnings RUSTFLAGS: '-Cinstrument-coverage -Clink-dead-code' run: | + export NEOPDF_DATA_PATH=/usr/local/share/LHAPDF # we need stderr, but we can't run test twice because it'll regenerate/modify the binaries which interferes with `llvm-cov` cargo test --features=applgrid,evolve,fastnlo,fktable --no-fail-fast 2> >(tee stderr 1>&2) # from https://stackoverflow.com/a/51141872/812178 diff --git a/CHANGELOG.md b/CHANGELOG.md index df97b4e6..834d6604 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- added support for the [NeoPDF](https://qcdlab.github.io/neopdf/) library as + an alternative to LHAPDF in the CLI; using `--backend=neopdf` one can use it. + ## [1.4.0] - 12/05/2026 Starting with this version, PineAPPL has an official logo! diff --git a/Cargo.lock b/Cargo.lock index 1ff40b9e..155aa2a8 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3,10 +3,10 @@ version = 4 [[package]] -name = "adler" -version = "1.0.2" +name = "adler2" +version = "2.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f26201604c87b1e01bd3d98f8d5d9a8fcbb815e8cedb41ffccbeb4bf593a35fe" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" [[package]] name = "aho-corasick" @@ -246,7 +246,7 @@ checksum = "af491d569909a7e4dee0ad7db7f5341fef5c614d5b8ec8cf765732aba3cff681" dependencies = [ "serde", "termcolor", - "unicode-width", + "unicode-width 0.2.2", ] [[package]] @@ -255,6 +255,17 @@ version = "1.0.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1d07550c9036bf2ae0c684c4297d503f838287c83c53686d05370d0e139ae570" +[[package]] +name = "console" +version = "0.16.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d64e8af5551369d19cf50138de61f1c42074ab970f74e99be916646777f8fc87" +dependencies = [ + "encode_unicode", + "libc", + "windows-sys 0.61.2", +] + [[package]] name = "cpufeatures" version = "0.2.17" @@ -445,6 +456,12 @@ version = "0.3.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "780955b8b195a21ab8e4ac6b60dd1dbdcec1dc6c51c0617964b08c81785e12c9" +[[package]] +name = "dyn-clone" +version = "1.0.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d0881ea181b1df73ff77ffaaf9c7544ecc11e82fba9b5f27b262a3c73a332555" + [[package]] name = "either" version = "1.15.0" @@ -485,7 +502,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "39cab71617ae0d63f51a36d69f866391735b51691dbda63cf6f96d042b63efeb" dependencies = [ "libc", - "windows-sys 0.52.0", + "windows-sys 0.59.0", ] [[package]] @@ -512,9 +529,9 @@ checksum = "5baebc0774151f905a1a2cc41989300b1e6fbb29aff0ceffa1064fdd3088d582" [[package]] name = "flate2" -version = "1.0.28" +version = "1.1.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "46303f565772937ffe1d394a4fac6f411c6013172fadde9dcdb1e147a086940e" +checksum = "843fba2746e448b37e26a819579957415c8cef339bf08564fe8b7ddbd959573c" dependencies = [ "crc32fast", "miniz_oxide", @@ -667,6 +684,17 @@ dependencies = [ "hashbrown", ] +[[package]] +name = "indicatif" +version = "0.18.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "25470f23803092da7d239834776d653104d551bc4d7eacaf31e6837854b8e9eb" +dependencies = [ + "console", + "portable-atomic", + "unit-prefix", +] + [[package]] name = "is-terminal" version = "0.4.17" @@ -675,14 +703,14 @@ checksum = "3640c1c38b8e4e43584d8df18be5fc6b0aa314ce6ebf51b53313d4306cca8e46" dependencies = [ "hermit-abi", "libc", - "windows-sys 0.52.0", + "windows-sys 0.59.0", ] [[package]] name = "itertools" -version = "0.14.0" +version = "0.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2b192c782037fadd9cfa75548310488aabdbf3d2da73885b31bd0abd03351285" +checksum = "413ee7dfc52ee1a4949ceeb7dbc8a33f2d6c088194d9f922fb8318faf1f01186" dependencies = [ "either", ] @@ -782,11 +810,12 @@ checksum = "f8ca58f447f06ed17d5fc4043ce1b10dd205e060fb3ce5b979b8ed8e59ff3f79" [[package]] name = "miniz_oxide" -version = "0.7.1" +version = "0.8.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e7810e0be55b428ada41041c41f32c9f1a42817901b4ccf45fa3d4b6561e74c7" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" dependencies = [ - "adler", + "adler2", + "simd-adler32", ] [[package]] @@ -818,6 +847,66 @@ dependencies = [ "zip", ] +[[package]] +name = "neopdf" +version = "0.3.3" +source = "git+https://github.com/QCDLab/neopdf.git#0adbddd49ae3312e42f9200a5850b538b1bfb354" +dependencies = [ + "bincode", + "flate2", + "git-version", + "indicatif", + "itertools", + "lz4_flex", + "ndarray", + "neopdf_legacy", + "ninterp", + "rayon", + "regex", + "serde", + "serde_yaml", + "tar", + "tempfile", + "thiserror", + "ureq", +] + +[[package]] +name = "neopdf_legacy" +version = "0.3.3" +source = "git+https://github.com/QCDLab/neopdf.git#0adbddd49ae3312e42f9200a5850b538b1bfb354" +dependencies = [ + "bincode", + "flate2", + "git-version", + "indicatif", + "itertools", + "lz4_flex", + "ndarray", + "ninterp", + "rayon", + "regex", + "serde", + "serde_yaml", + "tar", + "tempfile", + "thiserror", + "ureq", +] + +[[package]] +name = "ninterp" +version = "0.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "67588d36f7a821f32edf2ad9699804b65f65b835634ccc267f1c2bcdaca1acf4" +dependencies = [ + "dyn-clone", + "itertools", + "ndarray", + "num-traits", + "thiserror", +] + [[package]] name = "num-bigint" version = "0.4.6" @@ -997,6 +1086,7 @@ dependencies = [ "managed-lhapdf", "ndarray", "ndarray-npy", + "neopdf", "pineappl", "pineappl_applgrid", "pineappl_fastnlo", @@ -1098,7 +1188,7 @@ dependencies = [ "is-terminal", "lazy_static", "term", - "unicode-width", + "unicode-width 0.1.11", ] [[package]] @@ -1251,6 +1341,18 @@ dependencies = [ "thiserror", ] +[[package]] +name = "regex" +version = "1.12.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e10754a14b9137dd7b1e3e5b0493cc9171fdd105e0ab477f51b72e7f3ac0e276" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + [[package]] name = "regex-automata" version = "0.4.14" @@ -1304,7 +1406,7 @@ dependencies = [ "errno", "libc", "linux-raw-sys", - "windows-sys 0.52.0", + "windows-sys 0.59.0", ] [[package]] @@ -1438,6 +1540,12 @@ version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" +[[package]] +name = "simd-adler32" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "703d5c7ef118737c72f1af64ad2f6f8c5e1921f818cdcb97b8fe6fc69bf66214" + [[package]] name = "smallvec" version = "1.15.1" @@ -1606,6 +1714,18 @@ version = "0.1.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e51733f11c9c4f72aa0c160008246859e340b00807569a0da0e7a1079b27ba85" +[[package]] +name = "unicode-width" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b4ac048d71ede7ee76d585517add45da530660ef4390e49b098733c6e897f254" + +[[package]] +name = "unit-prefix" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "81e544489bf3d8ef66c953931f56617f423cd4b5494be343d9b9d3dda037b9a3" + [[package]] name = "unsafe-libyaml" version = "0.2.11" @@ -1721,7 +1841,7 @@ version = "0.1.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" dependencies = [ - "windows-sys 0.48.0", + "windows-sys 0.59.0", ] [[package]] @@ -1730,6 +1850,12 @@ version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + [[package]] name = "windows-sys" version = "0.48.0" @@ -1748,6 +1874,24 @@ dependencies = [ "windows-targets 0.52.6", ] +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets 0.52.6", +] + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + [[package]] name = "windows-targets" version = "0.48.5" diff --git a/Cargo.toml b/Cargo.toml index 4e8902d5..98c8e926 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -35,11 +35,12 @@ enum_dispatch = "0.3.7" flate2 = "1.0.22" float-cmp = { default-features = false, version = "0.10.0" } git-version = "0.3.5" -itertools = "0.14.0" +itertools = "0.13.0" lhapdf = { package = "managed-lhapdf", version = "0.4.1" } lz4_flex = "0.13.0" ndarray = { version = "0.17.2" } ndarray-npy = { default-features = false, features = ["npz"], version = "0.10.0" } +neopdf = { git = "https://github.com/QCDLab/neopdf.git" } num-complex = "0.4.4" numpy = "0.28.0" pkg-config = "0.3.26" diff --git a/README.md b/README.md index f3d5e2e8..65c8f80a 100644 --- a/README.md +++ b/README.md @@ -45,7 +45,8 @@ If you use PineAPPL, please cite By using PineAPPL, you're probably also using - [APPLgrid], -- [fastNLO] and/or +- [fastNLO], +- [NeoPDF] and/or - [LHAPDF]. If that is the case, please cite these accordingly. @@ -53,6 +54,7 @@ If that is the case, please cite these accordingly. [APPLgrid]: https://applgrid.hepforge.org [fastNLO]: https://fastnlo.hepforge.org [LHAPDF]: https://lhapdf.hepforge.org +[NeoPDF]: https://qcdlab.github.io/neopdf/ [zenodo DOI]: https://zenodo.org/badge/latestdoi/248306479 [paper]: https://inspirehep.net/literature/1814432 diff --git a/pineappl_cli/Cargo.toml b/pineappl_cli/Cargo.toml index 65e8041e..854cca95 100644 --- a/pineappl_cli/Cargo.toml +++ b/pineappl_cli/Cargo.toml @@ -30,6 +30,7 @@ lhapdf.workspace = true lz4_flex = { optional = true, workspace = true } ndarray.workspace = true ndarray-npy = { optional = true, workspace = true } +neopdf.workspace = true pineappl.workspace = true pineappl_applgrid = { optional = true, workspace = true } pineappl_fastnlo = { optional = true, workspace = true } diff --git a/pineappl_cli/src/analyze.rs b/pineappl_cli/src/analyze.rs index cc5f8738..90abfedc 100644 --- a/pineappl_cli/src/analyze.rs +++ b/pineappl_cli/src/analyze.rs @@ -64,7 +64,8 @@ pub struct CkfOpts { impl Subcommand for CkfOpts { fn run(&self, cfg: &GlobalConfiguration) -> Result { let grid = helpers::read_grid(&self.input)?; - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs, cfg.pdf_backend)?; let orders_den = if self.orders_den.is_empty() { grid.orders() @@ -84,7 +85,7 @@ impl Subcommand for CkfOpts { .map(|lumi| { let mut lumi_mask = vec![false; grid.channels().len()]; lumi_mask[lumi] = true; - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs.conv_types, @@ -101,7 +102,7 @@ impl Subcommand for CkfOpts { .map(|lumi| { let mut lumi_mask = vec![false; grid.channels().len()]; lumi_mask[lumi] = true; - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs.conv_types, diff --git a/pineappl_cli/src/channels.rs b/pineappl_cli/src/channels.rs index 0ba04375..247c54dd 100644 --- a/pineappl_cli/src/channels.rs +++ b/pineappl_cli/src/channels.rs @@ -64,7 +64,8 @@ pub struct Opts { impl Subcommand for Opts { fn run(&self, cfg: &GlobalConfiguration) -> Result { let grid = helpers::read_grid(&self.input)?; - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs, cfg.pdf_backend)?; let mut channels: Vec<_> = self.channels.iter().cloned().flatten().collect(); channels.sort_unstable(); @@ -90,7 +91,7 @@ impl Subcommand for Opts { .map(|channel| { let mut channel_mask = vec![false; grid.channels().len()]; channel_mask[channel] = true; - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs.conv_types, diff --git a/pineappl_cli/src/convolve.rs b/pineappl_cli/src/convolve.rs index f6ceae11..78728562 100644 --- a/pineappl_cli/src/convolve.rs +++ b/pineappl_cli/src/convolve.rs @@ -58,10 +58,11 @@ pub struct Opts { impl Subcommand for Opts { fn run(&self, cfg: &GlobalConfiguration) -> Result { let grid = helpers::read_grid(&self.input)?; - let mut conv_funs_0 = helpers::create_conv_funs(&self.conv_funs[0])?; + let mut conv_funs_0 = + helpers::create_conv_funs_with_backend(&self.conv_funs[0], cfg.pdf_backend)?; let bins: Vec<_> = self.bins.iter().cloned().flatten().collect(); - let results = helpers::convolve_scales( + let results = helpers::convolve_scales_with_backend( &grid, &mut conv_funs_0, &self.conv_funs[0].conv_types, @@ -91,8 +92,9 @@ impl Subcommand for Opts { .iter() .flat_map(|conv_funs| { let conv_types = &conv_funs.conv_types; - let mut conv_funs = helpers::create_conv_funs(conv_funs).unwrap(); - helpers::convolve( + let mut conv_funs = + helpers::create_conv_funs_with_backend(conv_funs, cfg.pdf_backend).unwrap(); + helpers::convolve_with_backend( &grid, &mut conv_funs, conv_types, diff --git a/pineappl_cli/src/diff.rs b/pineappl_cli/src/diff.rs index 151549be..5b1e7b05 100644 --- a/pineappl_cli/src/diff.rs +++ b/pineappl_cli/src/diff.rs @@ -120,7 +120,8 @@ impl Subcommand for Opts { bail!("channels differ"); } - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs, cfg.pdf_backend)?; let mut table = helpers::create_table(); let mut title = Row::empty(); @@ -140,7 +141,7 @@ impl Subcommand for Opts { table.set_titles(title); - let results1 = helpers::convolve( + let results1 = helpers::convolve_with_backend( &grid1, &mut conv_funs, &self.conv_funs.conv_types, @@ -151,7 +152,7 @@ impl Subcommand for Opts { ConvoluteMode::Normal, cfg, ); - let results2 = helpers::convolve( + let results2 = helpers::convolve_with_backend( &grid2, &mut conv_funs, &self.conv_funs.conv_types, @@ -204,7 +205,7 @@ impl Subcommand for Opts { let order_results1: Vec> = orders .iter() .map(|&order| { - helpers::convolve( + helpers::convolve_with_backend( &grid1, &mut conv_funs, &self.conv_funs.conv_types, @@ -220,7 +221,7 @@ impl Subcommand for Opts { let order_results2: Vec> = orders .iter() .map(|&order| { - helpers::convolve( + helpers::convolve_with_backend( &grid2, &mut conv_funs, &self.conv_funs.conv_types, diff --git a/pineappl_cli/src/evolve/mod.rs b/pineappl_cli/src/evolve/mod.rs index d3adf6ff..023819db 100644 --- a/pineappl_cli/src/evolve/mod.rs +++ b/pineappl_cli/src/evolve/mod.rs @@ -2,10 +2,10 @@ mod eko; use super::helpers::{self, ConvFuns, ConvoluteMode}; +use super::pdf_backend::PdfBackend; use super::{GlobalConfiguration, Subcommand}; use anyhow::{Result, anyhow}; use clap::{Parser, ValueHint}; -use lhapdf::Pdf; use pineappl::fk_table::FkTable; use pineappl::grid::Grid; use std::path::{Path, PathBuf}; @@ -59,8 +59,9 @@ impl Subcommand for Opts { use prettytable::row; let grid = helpers::read_grid(&self.input)?; - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?; - let results = helpers::convolve_scales( + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs, cfg.pdf_backend)?; + let results = helpers::convolve_scales_with_backend( &grid, &mut conv_funs, &self.conv_funs.conv_types, @@ -76,14 +77,14 @@ impl Subcommand for Opts { let fk_table = evolve_grid( &grid, &eko_paths, - &conv_funs[cfg.use_alphas_from], + &*conv_funs[cfg.use_alphas_from], &self.orders, self.xir, self.xif, self.xia, )?; - let evolved_results = helpers::convolve_scales( + let evolved_results = helpers::convolve_scales_with_backend( fk_table.grid(), &mut conv_funs, &self.conv_funs.conv_types, @@ -137,7 +138,7 @@ impl Subcommand for Opts { fn evolve_grid( grid: &Grid, ekos: &[&Path], - use_alphas_from: &Pdf, + use_alphas_from: &dyn PdfBackend, orders: &[(u8, u8)], xir: f64, xif: f64, @@ -171,7 +172,7 @@ fn evolve_grid( fn evolve_grid( _: &Grid, _: &[&Path], - _: &Pdf, + _: &dyn PdfBackend, _: &[(u8, u8)], _: f64, _: f64, diff --git a/pineappl_cli/src/helpers.rs b/pineappl_cli/src/helpers.rs index f618af41..c0bcbb76 100644 --- a/pineappl_cli/src/helpers.rs +++ b/pineappl_cli/src/helpers.rs @@ -1,7 +1,8 @@ use super::GlobalConfiguration; +use super::pdf_backend::{self, Backend, ForcePositive, PdfBackend, PdfSetBackend}; use anyhow::{Context as _, Error, Result, anyhow, bail}; use itertools::Itertools as _; -use lhapdf::{Pdf, PdfSet}; +use lhapdf::Pdf; use pineappl::boc::{ScaleFuncForm, Scales}; use pineappl::convolutions::{Conv, ConvType, ConvolutionCache}; use pineappl::grid::Grid; @@ -128,13 +129,14 @@ impl FromStr for ConvFuns { } } -#[derive(Clone, Copy)] +#[derive(Copy, Clone)] pub enum ConvoluteMode { Asymmetry, Integrated, Normal, } +/// Creates convolution functions using LHAPDF backend (legacy). pub fn create_conv_funs(funs: &ConvFuns) -> Result> { Ok(funs .lhapdf_names @@ -153,32 +155,38 @@ pub fn create_conv_funs(funs: &ConvFuns) -> Result> { .collect::>()?) } +/// Creates convolution functions using the specified backend. +pub fn create_conv_funs_with_backend( + funs: &ConvFuns, + backend: Backend, +) -> Result>> { + funs.lhapdf_names + .iter() + .zip(&funs.members) + .map(|(name, member)| { + let member = member.unwrap_or(0); + pdf_backend::create_pdf(name, member, backend) + }) + .collect() +} + +/// Creates convolution functions for a PDF set using the specified backend. +/// +/// Returns a tuple of (`PdfSetBackend`, Vec>>). pub fn create_conv_funs_for_set( funs: &ConvFuns, index_of_set: usize, -) -> Result<(PdfSet, Vec>)> { + backend: Backend, +) -> Result<(Box, Vec>>)> { let setname = &funs.lhapdf_names[index_of_set]; - let set = setname.parse().map_or_else( - |_| Ok::<_, Error>(PdfSet::new(setname)?), - |lhaid| { - Ok(PdfSet::new( - &lhapdf::lookup_pdf(lhaid) - .map(|(set, _)| set) - .ok_or_else(|| { - anyhow!("no convolution function for LHAID = `{lhaid}` found") - })?, - )?) - }, - )?; + let set = pdf_backend::create_pdf_set(setname, backend)?; - let conv_funs = set - .mk_pdfs()? + let set_members = set.mk_pdfs()?; + let conv_funs = set_members .into_iter() - .map(|conv_fun| { - // TODO: do not create objects that are getting overwritten in any case - let mut conv_funs = create_conv_funs(funs)?; - conv_funs[index_of_set] = conv_fun; - + .map(|member_pdf| { + let mut conv_funs = create_conv_funs_with_backend(funs, backend)?; + conv_funs[index_of_set] = member_pdf; Ok::<_, Error>(conv_funs) }) .collect::>()?; @@ -248,6 +256,7 @@ pub fn labels_and_units(grid: &Grid, integrated: bool) -> (Vec<(String, &str)>, ) } +/// Performs convolution with scale variations using LHAPDF backend (legacy). pub fn convolve_scales( grid: &Grid, conv_funs: &mut [Pdf], @@ -328,6 +337,125 @@ pub fn convolve_scales( let mut cache = ConvolutionCache::new(convolutions, xfx, &mut alphas_funs[cfg.use_alphas_from]); let mut results = grid.convolve(&mut cache, &orders, bins, channels, scales); + match mode { + ConvoluteMode::Asymmetry => { + let bin_count = grid.bwfl().len(); + assert!((bins.is_empty() || (bins.len() == bin_count)) && (bin_count % 2 == 0)); + + results + .iter() + .skip((bin_count / 2) * scales.len()) + .zip( + results + .chunks_exact(scales.len()) + .take(bin_count / 2) + .rev() + .flatten(), + ) + .map(|(pos, neg)| (pos - neg) / (pos + neg)) + .collect() + } + ConvoluteMode::Integrated => { + results + .iter_mut() + .zip( + grid.bwfl() + .normalizations() + .into_iter() + .enumerate() + .filter(|(index, _)| bins.is_empty() || bins.contains(index)) + .flat_map(|(_, norm)| iter::repeat(norm).take(scales.len())), + ) + .for_each(|(value, norm)| *value *= norm); + + results + } + ConvoluteMode::Normal => results, + } +} + +/// Performs convolution with scale variations using the backend abstraction. +#[allow(clippy::too_many_arguments)] +pub fn convolve_scales_with_backend( + grid: &Grid, + conv_funs: &mut [Box], + conv_types: &[ConvType], + orders: &[(u8, u8)], + bins: &[usize], + channels: &[bool], + scales: &[(f64, f64, f64)], + mode: ConvoluteMode, + cfg: &GlobalConfiguration, +) -> Vec { + let orders: Vec<_> = grid + .orders() + .iter() + .map(|order| { + orders.is_empty() + || orders + .iter() + .any(|other| (order.alphas == other.0) && (order.alpha == other.1)) + }) + .collect(); + + if cfg.force_positive { + for fun in conv_funs.iter_mut() { + fun.set_force_positive(ForcePositive::ClipNegative); + } + } + + // TODO: promote this to an error + assert!( + cfg.use_alphas_from < conv_funs.len(), + "expected `use_alphas_from` to be an integer within `[0, {})`, but got `{}`", + conv_funs.len(), + cfg.use_alphas_from + ); + + let x_min_max: Vec<_> = conv_funs + .iter_mut() + .map(|fun| (fun.x_min(), fun.x_max())) + .collect(); + + let mut funs: Vec<_> = conv_funs + .iter() + .zip(&x_min_max) + .map(|(fun, &(x_min, x_max))| { + move |id: i32, x: f64, q2: f64| { + if !cfg.allow_extrapolation && (x < x_min || x > x_max) { + 0.0 + } else { + fun.xfx_q2(id, x, q2) + } + } + }) + .collect(); + + let xfx: Vec<_> = funs + .iter_mut() + .map(|fun| fun as &mut dyn FnMut(i32, f64, f64) -> f64) + .collect(); + + let mut alphas_funs: Vec<_> = conv_funs + .iter() + .map(|fun| { + let fun_ref = fun.as_ref(); + move |q2: f64| fun_ref.alphas_q2(q2) + }) + .collect(); + + let convolutions: Vec<_> = conv_funs + .iter() + .zip(conv_types) + .map(|(fun, &conv_type)| { + let pid = fun.particle_id(); + Conv::new(conv_type, pid) + }) + .collect(); + + let mut cache = ConvolutionCache::new(convolutions, xfx, &mut alphas_funs[cfg.use_alphas_from]); + let mut results = grid.convolve(&mut cache, &orders, bins, channels, scales); + match mode { ConvoluteMode::Asymmetry => { let bin_count = grid.bwfl().len(); @@ -383,6 +511,7 @@ pub fn scales_vector(grid: &Grid, scales: usize) -> &[(f64, f64, f64)] { } } +/// Performs convolution using LHAPDF backend (legacy). pub fn convolve( grid: &Grid, conv_funs: &mut [Pdf], @@ -407,6 +536,32 @@ pub fn convolve( ) } +/// Performs convolution using the backend abstraction. +#[allow(clippy::too_many_arguments)] +pub fn convolve_with_backend( + grid: &Grid, + conv_funs: &mut [Box], + conv_types: &[ConvType], + orders: &[(u8, u8)], + bins: &[usize], + lumis: &[bool], + scales: usize, + mode: ConvoluteMode, + cfg: &GlobalConfiguration, +) -> Vec { + convolve_scales_with_backend( + grid, + conv_funs, + conv_types, + orders, + bins, + lumis, + scales_vector(grid, scales), + mode, + cfg, + ) +} + pub fn convolve_limits(grid: &Grid, bins: &[usize], mode: ConvoluteMode) -> Vec> { let limits: Vec<_> = grid .bwfl() diff --git a/pineappl_cli/src/lib.rs b/pineappl_cli/src/lib.rs index 50abbb09..98c50987 100644 --- a/pineappl_cli/src/lib.rs +++ b/pineappl_cli/src/lib.rs @@ -11,6 +11,7 @@ mod helpers; mod import; mod merge; mod orders; +pub mod pdf_backend; mod plot; mod pull; mod read; @@ -39,6 +40,9 @@ pub struct GlobalConfiguration { /// Choose the PDF/FF set for the strong coupling. #[arg(default_value = "0", long, value_name = "IDX")] pub use_alphas_from: usize, + /// Select the PDF interpolation backend: 'lhapdf' or 'neopdf'. + #[arg(default_value = "lhapdf", long, value_name = "BACKEND")] + pub pdf_backend: pdf_backend::Backend, } /// TODO. diff --git a/pineappl_cli/src/orders.rs b/pineappl_cli/src/orders.rs index 7000d738..efc42a87 100644 --- a/pineappl_cli/src/orders.rs +++ b/pineappl_cli/src/orders.rs @@ -41,7 +41,8 @@ pub struct Opts { impl Subcommand for Opts { fn run(&self, cfg: &GlobalConfiguration) -> Result { let grid = helpers::read_grid(&self.input)?; - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs, cfg.pdf_backend)?; let mut orders: Vec<_> = grid .orders() @@ -63,7 +64,7 @@ impl Subcommand for Opts { let results: Vec> = orders .iter() .map(|order| { - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs.conv_types, diff --git a/pineappl_cli/src/pdf_backend.rs b/pineappl_cli/src/pdf_backend.rs new file mode 100644 index 00000000..b652315e --- /dev/null +++ b/pineappl_cli/src/pdf_backend.rs @@ -0,0 +1,432 @@ +//! PDF backend abstraction layer. +//! +//! This module provides a unified interface for different PDF interpolation backends, +//! currently supporting `LHAPDF` and `NeoPDF`. It allows runtime selection of the backend +//! and provides type-safe access to PDF metadata. + +use anyhow::{Context, Result, anyhow}; +// use std::fmt; +use std::str::FromStr; + +pub use neopdf::uncertainty::CL_1_SIGMA; + +// /// The type of PDF set (space-like for PDFs, time-like for FFs). +// #[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +// pub enum SetType { +// /// Space-like PDF (standard parton distribution function). +// #[default] +// SpaceLike, +// /// Time-like PDF (fragmentation function). +// TimeLike, +// } + +// impl fmt::Display for SetType { +// fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { +// match self { +// Self::SpaceLike => write!(f, "SpaceLike"), +// Self::TimeLike => write!(f, "TimeLike"), +// } +// } +// } + +/// Method for handling negative PDF values. +#[derive(Clone, Copy, Debug, Default)] +pub enum ForcePositive { + /// No clipping - return values as-is. + #[default] + None, + /// Clip negative values to zero. + ClipNegative, +} + +/// Available PDF backends. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum Backend { + /// `LHAPDF` backend (C++ library with Rust bindings). + #[default] + Lhapdf, + /// `NeoPDF` backend (pure Rust implementation). + Neopdf, +} + +impl FromStr for Backend { + type Err = anyhow::Error; + + fn from_str(s: &str) -> Result { + match s.to_lowercase().as_str() { + "lhapdf" => Ok(Self::Lhapdf), + "neopdf" => Ok(Self::Neopdf), + _ => Err(anyhow!( + "unknown PDF backend '{}'; must be 'lhapdf' or 'neopdf'", + s + )), + } + } +} + +// impl fmt::Display for Backend { +// fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { +// match self { +// Self::Lhapdf => write!(f, "lhapdf"), +// Self::Neopdf => write!(f, "neopdf"), +// } +// } +// } + +/// Unified PDF backend interface. +/// +/// This trait provides a common interface for PDF interpolation backends, +/// allowing the CLI to work with different implementations transparently. +pub trait PdfBackend: Send { + /// Evaluates xf(x, Q^2) for a given flavor. + /// + /// # Arguments + /// * `id` - The Monte Carlo PDG flavor ID. + /// * `x` - The momentum fraction. + /// * `q2` - The squared energy scale. + /// + /// # Returns + /// The PDF value xf(x, Q^2). + /// + /// # TODO + /// Extend to support `NeoPDF` multi-parameters interpolation. + fn xfx_q2(&self, id: i32, x: f64, q2: f64) -> f64; + + /// Evaluates the strong coupling constant `alpha_s(Q^2)`. + fn alphas_q2(&self, q2: f64) -> f64; + + /// Returns the minimum valid x value. + fn x_min(&mut self) -> f64; + + /// Returns the maximum valid x value. + fn x_max(&mut self) -> f64; + + /// Returns the hadron particle ID (PDG code). + fn particle_id(&self) -> i32; + + // /// Returns whether this is a polarized PDF. + // fn is_polarized(&self) -> bool; + + // /// Returns the PDF set type (space-like or time-like). + // fn set_type(&self) -> SetType; + + /// Sets the method for handling negative PDF values. + fn set_force_positive(&mut self, method: ForcePositive); + + // /// Returns the error type of the PDF set (e.g., "replicas", "hessian"). + // fn error_type(&self) -> String; +} + +/// Unified PDF set interface for uncertainty calculations. +pub trait PdfSetBackend { + // /// Returns the number of members in this set. + // fn num_members(&self) -> usize; + + /// Creates all PDF members in the set. + fn mk_pdfs(&self) -> Result>>; + + /// Calculates the uncertainty for a set of values. + fn uncertainty(&self, values: &[f64], cl: f64, alternative: bool) -> Result; +} + +/// Uncertainty information from a PDF set. +#[derive(Clone, Debug)] +pub struct Uncertainty { + /// Central value. + pub central: f64, + /// Negative error (absolute value). + pub errminus: f64, + /// Positive error (absolute value). + pub errplus: f64, +} + +// ============================================================================ +// LHAPDF Backend Implementation +// ============================================================================ + +/// LHAPDF backend wrapper. +pub struct LhapdfPdf { + pdf: lhapdf::Pdf, +} + +impl LhapdfPdf { + /// Creates a new LHAPDF PDF by set name and member index. + pub fn with_setname_and_member(setname: &str, member: i32) -> Result { + Ok(Self { + pdf: lhapdf::Pdf::with_setname_and_member(setname, member).with_context(|| { + format!("failed to load LHAPDF set '{setname}' member {member}") + })?, + }) + } + + /// Creates a new LHAPDF PDF by LHAID. + pub fn with_lhaid(lhaid: i32) -> Result { + Ok(Self { + pdf: lhapdf::Pdf::with_lhaid(lhaid) + .with_context(|| format!("failed to load LHAPDF with LHAID {lhaid}"))?, + }) + } +} + +impl PdfBackend for LhapdfPdf { + fn xfx_q2(&self, id: i32, x: f64, q2: f64) -> f64 { + self.pdf.xfx_q2(id, x, q2) + } + + fn alphas_q2(&self, q2: f64) -> f64 { + self.pdf.alphas_q2(q2) + } + + fn x_min(&mut self) -> f64 { + self.pdf.x_min() + } + + fn x_max(&mut self) -> f64 { + self.pdf.x_max() + } + + fn particle_id(&self) -> i32 { + self.pdf + .set() + .entry("Particle") + .map_or(Ok(2212), |s| s.parse()) + .unwrap_or(2212) + } + + // fn is_polarized(&self) -> bool { + // self.pdf + // .set() + // .entry("Polarized") + // .is_some_and(|s| s.eq_ignore_ascii_case("true")) + // } + + // fn set_type(&self) -> SetType { + // self.pdf + // .set() + // .entry("SetType") + // .map_or(SetType::SpaceLike, |s| { + // if s.eq_ignore_ascii_case("timelike") { + // SetType::TimeLike + // } else { + // SetType::SpaceLike + // } + // }) + // } + + fn set_force_positive(&mut self, method: ForcePositive) { + match method { + ForcePositive::None => self.pdf.set_force_positive(0), + ForcePositive::ClipNegative => self.pdf.set_force_positive(1), + } + } + + // fn error_type(&self) -> String { + // self.pdf + // .set() + // .entry("ErrorType") + // .map(|s| s.to_owned()) + // .unwrap_or_default() + // } +} + +/// LHAPDF set wrapper. +pub struct LhapdfSet { + set: lhapdf::PdfSet, +} + +impl LhapdfSet { + /// Creates a new LHAPDF set by name. + pub fn new(setname: &str) -> Result { + Ok(Self { + set: lhapdf::PdfSet::new(setname) + .with_context(|| format!("failed to load LHAPDF set '{setname}'"))?, + }) + } +} + +impl PdfSetBackend for LhapdfSet { + // fn num_members(&self) -> usize { + // self.set + // .entry("NumMembers") + // .and_then(|s| s.parse().ok()) + // .unwrap_or(1) + // } + + fn mk_pdfs(&self) -> Result>> { + let pdfs = self.set.mk_pdfs()?; + Ok(pdfs + .into_iter() + .map(|pdf| Box::new(LhapdfPdf { pdf }) as Box) + .collect()) + } + + fn uncertainty(&self, values: &[f64], cl: f64, alternative: bool) -> Result { + let unc = self.set.uncertainty(values, cl, alternative)?; + Ok(Uncertainty { + central: unc.central, + errminus: unc.errminus, + errplus: unc.errplus, + }) + } +} + +// ============================================================================ +// NeoPDF Backend Implementation +// ============================================================================ + +/// `NeoPDF` backend wrapper. +pub struct NeopdfPdf { + pdf: neopdf::pdf::PDF, +} + +impl NeopdfPdf { + /// Loads a NeoPDF member by set name and member index. + pub fn load(pdf_name: &str, member: usize) -> Self { + Self { + pdf: neopdf::pdf::PDF::load(pdf_name, member), + } + } +} + +impl PdfBackend for NeopdfPdf { + fn xfx_q2(&self, id: i32, x: f64, q2: f64) -> f64 { + self.pdf.xfxq2(id, &[x, q2]) + } + + fn alphas_q2(&self, q2: f64) -> f64 { + self.pdf.alphas_q2(q2) + } + + fn x_min(&mut self) -> f64 { + self.pdf.metadata().x_min + } + + fn x_max(&mut self) -> f64 { + self.pdf.metadata().x_max + } + + fn particle_id(&self) -> i32 { + self.pdf.metadata().hadron_pid + } + + // fn is_polarized(&self) -> bool { + // self.pdf.metadata().polarised + // } + + // fn set_type(&self) -> SetType { + // match self.pdf.metadata().set_type { + // neopdf::metadata::SetType::SpaceLike => SetType::SpaceLike, + // neopdf::metadata::SetType::TimeLike => SetType::TimeLike, + // } + // } + + fn set_force_positive(&mut self, method: ForcePositive) { + match method { + ForcePositive::None => { + self.pdf + .set_force_positive(neopdf::gridpdf::ForcePositive::NoClipping); + } + ForcePositive::ClipNegative => { + self.pdf + .set_force_positive(neopdf::gridpdf::ForcePositive::ClipNegative); + } + } + } + + // fn error_type(&self) -> String { + // self.pdf.metadata().error_type.clone() + // } +} + +/// `NeoPDF` set wrapper. +pub struct NeopdfSet { + pdf_name: String, + // num_members: usize, + error_type: String, +} + +impl NeopdfSet { + /// Creates a new NeoPDF set by name. + pub fn new(pdf_name: &str) -> Self { + // Load member 0 to get metadata + let pdf0 = neopdf::pdf::PDF::load(pdf_name, 0); + let metadata = pdf0.metadata(); + + Self { + pdf_name: pdf_name.to_owned(), + // num_members: metadata.num_members as usize, + error_type: metadata.error_type.clone(), + } + } +} + +impl PdfSetBackend for NeopdfSet { + // fn num_members(&self) -> usize { + // self.num_members + // } + + fn mk_pdfs(&self) -> Result>> { + let pdfs = neopdf::pdf::PDF::load_pdfs(&self.pdf_name); + Ok(pdfs + .into_iter() + .map(|pdf| Box::new(NeopdfPdf { pdf }) as Box) + .collect()) + } + + fn uncertainty(&self, values: &[f64], cl: f64, alternative: bool) -> Result { + let unc = neopdf::uncertainty::uncertainty( + values, + &self.error_type, + neopdf::uncertainty::CL_1_SIGMA, + cl, + alternative, + ) + .map_err(|e| anyhow::anyhow!(e))?; + + Ok(Uncertainty { + central: unc.central, + errminus: unc.errminus, + errplus: unc.errplus, + }) + } +} + +// ============================================================================ +// Factory Functions +// ============================================================================ + +/// Creates a single PDF from a set name and member index. +pub fn create_pdf(name: &str, member: usize, backend: Backend) -> Result> { + match backend { + Backend::Lhapdf => { + if let Ok(lhaid) = name.parse::() { + Ok(Box::new(LhapdfPdf::with_lhaid(lhaid)?)) + } else { + Ok(Box::new(LhapdfPdf::with_setname_and_member( + name, + member.try_into().unwrap(), + )?)) + } + } + // NOTE: `NeoPDF` doesn't support reading LHAID yet. + Backend::Neopdf => Ok(Box::new(NeopdfPdf::load(name, member))), + } +} + +/// Creates a PDF set for uncertainty calculations. +pub fn create_pdf_set(name: &str, backend: Backend) -> Result> { + match backend { + Backend::Lhapdf => { + // Try parsing as LHAID first + let setname = if let Ok(lhaid) = name.parse::() { + lhapdf::lookup_pdf(lhaid) + .map(|(set, _)| set) + .ok_or_else(|| anyhow!("no convolution function for LHAID = `{lhaid}` found"))? + } else { + name.to_owned() + }; + Ok(Box::new(LhapdfSet::new(&setname)?)) + } + Backend::Neopdf => Ok(Box::new(NeopdfSet::new(name))), + } +} diff --git a/pineappl_cli/src/plot.rs b/pineappl_cli/src/plot.rs index d661a016..6a0e5325 100644 --- a/pineappl_cli/src/plot.rs +++ b/pineappl_cli/src/plot.rs @@ -63,7 +63,8 @@ impl Subcommand for Opts { }; let grid = helpers::read_grid(&self.input)?; - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs[0])?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs[0], cfg.pdf_backend)?; let slices = grid.bwfl().slices(); let mut data_string = String::new(); @@ -98,7 +99,7 @@ impl Subcommand for Opts { { let bins: Vec<_> = slice.collect(); - let results = helpers::convolve( + let results = helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs[0].conv_types, @@ -125,7 +126,7 @@ impl Subcommand for Opts { }) .collect(); - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs[0].conv_types, @@ -159,9 +160,10 @@ impl Subcommand for Opts { .map(|conv_funs| { if self.no_conv_fun_unc { let conv_types = &conv_funs.conv_types; - let mut conv_funs = helpers::create_conv_funs(conv_funs)?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(conv_funs, cfg.pdf_backend)?; - let results = helpers::convolve( + let results = helpers::convolve_with_backend( &grid, &mut conv_funs, conv_types, @@ -178,12 +180,13 @@ impl Subcommand for Opts { let (set, funs) = helpers::create_conv_funs_for_set( conv_funs, self.conv_fun_uncert_from, + cfg.pdf_backend, )?; let pdf_results: Vec<_> = funs .into_par_iter() .flat_map(|mut funs| { - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut funs, &conv_funs.conv_types, @@ -279,7 +282,7 @@ impl Subcommand for Opts { channel_mask[channel] = true; ( map_format_channel(&grid.channels()[channel], &grid), - helpers::convolve( + helpers::convolve_with_backend( &grid, &mut conv_funs, &self.conv_funs[0].conv_types, diff --git a/pineappl_cli/src/pull.rs b/pineappl_cli/src/pull.rs index 0b40586d..09f1635e 100644 --- a/pineappl_cli/src/pull.rs +++ b/pineappl_cli/src/pull.rs @@ -1,8 +1,8 @@ use super::helpers::{self, ConvFuns, ConvoluteMode}; +use super::pdf_backend::{PdfBackend, PdfSetBackend}; use super::{GlobalConfiguration, Subcommand}; use anyhow::{Error, Result}; use clap::{Parser, ValueHint}; -use lhapdf::{Pdf, PdfSet}; use prettytable::{Row, cell}; use rayon::{ThreadPoolBuilder, prelude::*}; use std::num::NonZeroUsize; @@ -53,9 +53,9 @@ impl Subcommand for Opts { let grid = helpers::read_grid(&self.input)?; let (set1, mut conv_funs1) = - helpers::create_conv_funs_for_set(&self.conv_funs1, self.pull_from)?; + helpers::create_conv_funs_for_set(&self.conv_funs1, self.pull_from, cfg.pdf_backend)?; let (set2, mut conv_funs2) = - helpers::create_conv_funs_for_set(&self.conv_funs2, self.pull_from)?; + helpers::create_conv_funs_for_set(&self.conv_funs2, self.pull_from, cfg.pdf_backend)?; ThreadPoolBuilder::new() .num_threads(self.threads) @@ -67,7 +67,7 @@ impl Subcommand for Opts { let results1: Vec<_> = conv_funs1 .par_iter_mut() .map(|funs| { - Ok::<_, Error>(helpers::convolve( + Ok::<_, Error>(helpers::convolve_with_backend( &grid, funs, &self.conv_funs1.conv_types, @@ -86,7 +86,7 @@ impl Subcommand for Opts { let results2: Vec<_> = conv_funs2 .par_iter_mut() .map(|funs| { - Ok::<_, Error>(helpers::convolve( + Ok::<_, Error>(helpers::convolve_with_backend( &grid, funs, &self.conv_funs2.conv_types, @@ -141,79 +141,81 @@ impl Subcommand for Opts { (diff / unc1.hypot(unc2), unc1, unc2) }; - let channel_results = - |conv_funs: &ConvFuns, pdfset: &mut [Vec], set: &PdfSet| -> Vec { - if let Some(member) = conv_funs.members[self.pull_from] { - (0..grid.channels().len()) - .map(|channel| { - let mut channel_mask = vec![false; grid.channels().len()]; - channel_mask[channel] = true; - match helpers::convolve( - &grid, - &mut pdfset[member], - &conv_funs.conv_types, - &self.orders, - &[bin], - &channel_mask, - 1, - ConvoluteMode::Normal, - cfg, - ) - .as_slice() - { - [value] => *value, - _ => unreachable!(), - } - }) - .collect() - } else { - let results: Vec<_> = pdfset - .iter_mut() - .flat_map(|fun| { - (0..grid.channels().len()) - .map(|channel| { - let mut channel_mask = vec![false; grid.channels().len()]; - channel_mask[channel] = true; - match helpers::convolve( - &grid, - fun, - &conv_funs.conv_types, - &self.orders, - &[bin], - &channel_mask, - 1, - ConvoluteMode::Normal, - cfg, - ) - .as_slice() - { - [value] => *value, - _ => unreachable!(), - } - }) - .collect::>() - }) - .collect(); - - (0..grid.channels().len()) - .map(|channel| { - let central: Vec<_> = results - .iter() - .skip(channel) - .step_by(grid.channels().len()) - .copied() - .collect(); - set.uncertainty(¢ral, self.cl, false).unwrap().central - }) - .collect() - } - }; + let channel_results = |conv_funs: &ConvFuns, + pdfset: &mut [Vec>], + set: &dyn PdfSetBackend| + -> Vec { + if let Some(member) = conv_funs.members[self.pull_from] { + (0..grid.channels().len()) + .map(|channel| { + let mut channel_mask = vec![false; grid.channels().len()]; + channel_mask[channel] = true; + match helpers::convolve_with_backend( + &grid, + &mut pdfset[member], + &conv_funs.conv_types, + &self.orders, + &[bin], + &channel_mask, + 1, + ConvoluteMode::Normal, + cfg, + ) + .as_slice() + { + [value] => *value, + _ => unreachable!(), + } + }) + .collect() + } else { + let results: Vec<_> = pdfset + .iter_mut() + .flat_map(|fun| { + (0..grid.channels().len()) + .map(|channel| { + let mut channel_mask = vec![false; grid.channels().len()]; + channel_mask[channel] = true; + match helpers::convolve_with_backend( + &grid, + fun, + &conv_funs.conv_types, + &self.orders, + &[bin], + &channel_mask, + 1, + ConvoluteMode::Normal, + cfg, + ) + .as_slice() + { + [value] => *value, + _ => unreachable!(), + } + }) + .collect::>() + }) + .collect(); + + (0..grid.channels().len()) + .map(|channel| { + let central: Vec<_> = results + .iter() + .skip(channel) + .step_by(grid.channels().len()) + .copied() + .collect(); + set.uncertainty(¢ral, self.cl, false).unwrap().central + }) + .collect() + } + }; let mut pull_tuples = if self.limit == 0 { Vec::new() } else { - let channel_results1 = channel_results(&self.conv_funs1, &mut conv_funs1, &set1); - let channel_results2 = channel_results(&self.conv_funs2, &mut conv_funs2, &set2); + let channel_results1 = channel_results(&self.conv_funs1, &mut conv_funs1, &*set1); + let channel_results2 = channel_results(&self.conv_funs2, &mut conv_funs2, &*set2); let pull_tuples: Vec<_> = channel_results2 .iter() diff --git a/pineappl_cli/src/uncert.rs b/pineappl_cli/src/uncert.rs index f116c496..52632724 100644 --- a/pineappl_cli/src/uncert.rs +++ b/pineappl_cli/src/uncert.rs @@ -66,7 +66,7 @@ pub struct Opts { #[command(flatten)] group: Group, /// Confidence level in per cent, for convolution function uncertainties. - #[arg(default_value_t = lhapdf::CL_1_SIGMA, long)] + #[arg(default_value_t = super::pdf_backend::CL_1_SIGMA, long)] cl: f64, /// Show integrated numbers (without bin widths) instead of differential ones. #[arg(long, short)] @@ -94,7 +94,8 @@ pub struct Opts { impl Subcommand for Opts { fn run(&self, cfg: &GlobalConfiguration) -> Result { let grid = helpers::read_grid(&self.input)?; - let mut conv_funs = helpers::create_conv_funs(&self.conv_funs)?; + let mut conv_funs = + helpers::create_conv_funs_with_backend(&self.conv_funs, cfg.pdf_backend)?; let limits = helpers::convolve_limits( &grid, @@ -116,11 +117,12 @@ impl Subcommand for Opts { .conv_fun .iter() .map(|&index| { - let (set, funs) = helpers::create_conv_funs_for_set(&self.conv_funs, index)?; + let (set, funs) = + helpers::create_conv_funs_for_set(&self.conv_funs, index, cfg.pdf_backend)?; let results: Vec<_> = funs .into_par_iter() .map(|mut funs| { - Ok::<_, Error>(helpers::convolve( + Ok::<_, Error>(helpers::convolve_with_backend( &grid, &mut funs, &self.conv_funs.conv_types, @@ -159,7 +161,7 @@ impl Subcommand for Opts { .max() .unwrap_or(1); let scale_tuples = helpers::scales_vector(&grid, scales_max); - let scale_results = helpers::convolve_scales( + let scale_results = helpers::convolve_scales_with_backend( &grid, &mut conv_funs, &self.conv_funs.conv_types, diff --git a/pineappl_cli/tests/convolve.rs b/pineappl_cli/tests/convolve.rs index 95c66fd9..28ad066e 100644 --- a/pineappl_cli/tests/convolve.rs +++ b/pineappl_cli/tests/convolve.rs @@ -536,6 +536,86 @@ fn issue_334() { .stdout(NO_CHANNELS_GRID_STR); } +#[test] +fn neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "convolve", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(str::contains(DEFAULT_STR)); +} + +#[test] +fn integrated_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "convolve", + "--integrated", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(INTEGRATED_STR); +} + +#[test] +fn three_pdfs_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "convolve", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed/0", + "NNPDF31_nlo_as_0118_luxqed/1", + "NNPDF31_nlo_as_0118_luxqed/2", + ]) + .assert() + .success() + .stdout(THREE_PDFS_STR); +} + +#[test] +fn bins_13567_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "convolve", + "--bins=1,3,5-7", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(BINS_13567_STR); +} + +#[test] +fn force_positive_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--force-positive", + "--pdf-backend=neopdf", + "convolve", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(FORCE_POSITIVE_STR); +} + #[test] fn lagrange_subgrid_v1() { Command::cargo_bin("pineappl") diff --git a/pineappl_cli/tests/main.rs b/pineappl_cli/tests/main.rs index 2663622e..2ddc5869 100644 --- a/pineappl_cli/tests/main.rs +++ b/pineappl_cli/tests/main.rs @@ -29,6 +29,7 @@ Options: --force-positive Forces negative PDF values to zero --allow-extrapolation Allow extrapolation of PDFs outside their region of validity --use-alphas-from Choose the PDF/FF set for the strong coupling [default: 0] + --pdf-backend Select the PDF interpolation backend: 'lhapdf' or 'neopdf' [default: lhapdf] -h, --help Print help -V, --version Print version "; diff --git a/pineappl_cli/tests/uncert.rs b/pineappl_cli/tests/uncert.rs index dfcd0067..30ac3c39 100644 --- a/pineappl_cli/tests/uncert.rs +++ b/pineappl_cli/tests/uncert.rs @@ -1,6 +1,7 @@ #![expect(missing_docs, reason = "non-public items will not be documented")] use assert_cmd::Command; +use predicates::str; use std::num::NonZeroUsize; use std::thread; @@ -204,6 +205,28 @@ const SCALE_ENV_9_STR: &str = "b etal dsig/detal 9pt-svar (env) 7 4 4.5 2.7517266e1 -5.36 5.22 "; +const NEOPDF_CONV_FUN_STR: &str = "-+----+----+-----------+-----------+---------+--------- +0 2 2.25 7.5459110e2 7.5461655e2 -1.14 1.14 +1 2.25 2.5 6.9028342e2 6.9027941e2 -1.16 1.16 +2 2.5 2.75 6.0025198e2 6.0022595e2 -1.18 1.18 +3 2.75 3 4.8552235e2 4.8548211e2 -1.22 1.22 +4 3 3.25 3.6195456e2 3.6191001e2 -1.27 1.27 +5 3.25 3.5 2.4586691e2 2.4582640e2 -1.35 1.35 +6 3.5 4 1.1586851e2 1.1584074e2 -1.51 1.51 +7 4 4.5 2.7517266e1 2.7504644e1 -2.77 2.77 +"; + +const NEOPDF_CONV_FUN_CL_90_STR: &str = "-+----+----+-----------+-----------+---------+--------- +0 2 2.25 7.5459110e2 7.5461655e2 -1.87 1.87 +1 2.25 2.5 6.9028342e2 6.9027941e2 -1.90 1.90 +2 2.5 2.75 6.0025198e2 6.0022595e2 -1.95 1.95 +3 2.75 3 4.8552235e2 4.8548211e2 -2.00 2.00 +4 3 3.25 3.6195456e2 3.6191001e2 -2.08 2.08 +5 3.25 3.5 2.4586691e2 2.4582640e2 -2.22 2.22 +6 3.5 4 1.1586851e2 1.1584074e2 -2.48 2.48 +7 4 4.5 2.7517266e1 2.7504644e1 -4.56 4.56 +"; + #[test] fn help() { Command::cargo_bin("pineappl") @@ -437,3 +460,88 @@ fn scale_env_9() { .success() .stdout(SCALE_ENV_9_STR); } + +#[test] +fn conv_fun_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "uncert", + "--conv-fun", + "--threads=1", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(str::contains(NEOPDF_CONV_FUN_STR)); +} + +#[test] +fn conv_fun_cl_90_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "uncert", + "--conv-fun", + "--cl=90", + "--threads=1", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(str::contains(NEOPDF_CONV_FUN_CL_90_STR)); +} + +#[test] +fn scale_env_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "uncert", + "--scale-env", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(str::contains(SCALE_ENV_STR)); +} + +#[test] +fn conv_fun_orders_a2_as1a2_neopdf_beckend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "uncert", + "--conv-fun=0", + "--orders=a2,as1a2", + "--threads=1", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(ORDERS_A2_AS1A2_STR); +} + +#[test] +fn scale_abs_9_neopdf_backend() { + Command::cargo_bin("pineappl") + .unwrap() + .args([ + "--pdf-backend=neopdf", + "uncert", + "--scale-abs=9", + "../test-data/LHCB_WP_7TEV_opt.pineappl.lz4", + "NNPDF31_nlo_as_0118_luxqed", + ]) + .assert() + .success() + .stdout(SCALE_ABS_9_STR); +}