diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 9a1b4d8a5..bee3e7adc 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -2,6 +2,7 @@ # See https://pre-commit.com/hooks.html for more hooks ci: autofix_prs: false + skip: [fmt-eko, fmt-ekore] repos: - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.5.0 @@ -42,6 +43,24 @@ repos: args: ["--add-ignore=D107,D105"] additional_dependencies: - toml + - repo: local + hooks: + - id: fmt-eko + name: fmt-eko + description: Format eko files with cargo fmt. + entry: cargo fmt --manifest-path crates/eko/Cargo.toml -- + language: system + files: ^crates/eko/.*\.rs$ + args: [] + - repo: local + hooks: + - id: fmt-ekore + name: fmt-ekore + description: Format ekore files with cargo fmt. + entry: cargo fmt --manifest-path crates/ekore/Cargo.toml -- + language: system + files: ^crates/ekore/.*\.rs$ + args: [] - repo: https://github.com/pre-commit/pre-commit rev: v3.6.0 hooks: diff --git a/crates/.gitignore b/crates/.gitignore new file mode 100644 index 000000000..ea8c4bf7f --- /dev/null +++ b/crates/.gitignore @@ -0,0 +1 @@ +/target diff --git a/crates/Cargo.lock b/crates/Cargo.lock new file mode 100644 index 000000000..b6d5d084d --- /dev/null +++ b/crates/Cargo.lock @@ -0,0 +1,204 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "ahash" +version = "0.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2c99f64d1e06488f620f932677e24bc6e2897582980441ae90a671415bd7ec2f" +dependencies = [ + "cfg-if", + "once_cell", + "version_check", +] + +[[package]] +name = "allocator-api2" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0942ffc6dcaadf03badf6e6a2d0228460359d5e34b57ccdc720b7382dfbd5ec5" + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "eko" +version = "0.1.0" +dependencies = [ + "ekore", + "katexit", + "num", +] + +[[package]] +name = "ekore" +version = "0.1.0" +dependencies = [ + "float-cmp", + "hashbrown", + "katexit", + "num", +] + +[[package]] +name = "float-cmp" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "98de4bbd547a563b716d8dfa9aad1cb19bfab00f4fa09a6a4ed21dbcf44ce9c4" +dependencies = [ + "num-traits", +] + +[[package]] +name = "hashbrown" +version = "0.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2c6201b9ff9fd90a5a3bac2e56a830d0caa509576f0e503818ee82c181b3437a" +dependencies = [ + "ahash", + "allocator-api2", +] + +[[package]] +name = "katexit" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb1304c448ce2c207c2298a34bc476ce7ae47f63c23fa2b498583b26be9bc88c" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "num" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b05180d69e3da0e530ba2a1dae5110317e49e3b7f3d41be227dc5f92e49ee7af" +dependencies = [ + "num-bigint", + "num-complex", + "num-integer", + "num-iter", + "num-rational", + "num-traits", +] + +[[package]] +name = "num-bigint" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f93ab6289c7b344a8a9f60f88d80aa20032336fe78da341afc91c8a2341fc75f" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0d21255c828d6f128a1e41534206671e8c3ea0c62f32291e808dc82cff17d" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-iter" +version = "0.1.43" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7d03e6c028c5dc5cac6e2dec0efda81fc887605bb3d884578bb6d6bf7514e252" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-bigint", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f30b0abd723be7e2ffca1272140fac1a2f084c77ec3e123c192b66af1ee9e6c2" +dependencies = [ + "autocfg", +] + +[[package]] +name = "once_cell" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" + +[[package]] +name = "proc-macro2" +version = "1.0.66" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "18fb31db3f9bddb2ea821cde30a9f70117e3f119938b5ee630b7403aa6e2ead9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "50f3b39ccfb720540debaa0164757101c08ecb8d326b15358ce76a62c7e85965" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "unicode-ident" +version = "1.0.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "301abaae475aa91687eb82514b328ab47a211a533026cb25fc3e519b86adfc3c" + +[[package]] +name = "version_check" +version = "0.9.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "49874b5167b65d7193b8aba1567f5c7d93d001cafc34600cee003eda787e483f" diff --git a/crates/Cargo.toml b/crates/Cargo.toml new file mode 100644 index 000000000..8a69f7202 --- /dev/null +++ b/crates/Cargo.toml @@ -0,0 +1,3 @@ +[workspace] + +members = ["eko", "ekore"] diff --git a/crates/eko/Cargo.toml b/crates/eko/Cargo.toml new file mode 100644 index 000000000..50aee6504 --- /dev/null +++ b/crates/eko/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "eko" +version = "0.1.0" +edition = "2021" + +[lib] +name = "ekors" +crate-type = ["cdylib"] + +[dependencies] +num = "0.4.1" +katexit = "0.1.4" +ekore = { version = "0.1.0", path = "../ekore" } diff --git a/crates/eko/pyproject.toml b/crates/eko/pyproject.toml new file mode 100644 index 000000000..2fff64504 --- /dev/null +++ b/crates/eko/pyproject.toml @@ -0,0 +1,16 @@ +[build-system] +requires = ["maturin>=1.1,<2.0"] +build-backend = "maturin" + +[project] +name = "ekors" +requires-python = ">=3.8" +classifiers = [ + "Programming Language :: Rust", + "Programming Language :: Python :: Implementation :: CPython", + "Programming Language :: Python :: Implementation :: PyPy", +] +dependencies = ["cffi"] + +[tool.maturin] +bindings = "cffi" diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs new file mode 100644 index 000000000..e6b3e44cc --- /dev/null +++ b/crates/eko/src/lib.rs @@ -0,0 +1,188 @@ +//! Interface to the eko Python package. + +use ekore; +use ekore::harmonics::cache::Cache; +use std::ffi::c_void; + +mod mellin; + +/// QCD intergration kernel inside quad. +#[no_mangle] +pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { + let args = *(rargs as *mut QuadQCDargs); + let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0); + // prepare gamma + let path = mellin::TalbotPath::new(u, args.logx, is_singlet); + let jac = path.jac() * path.prefactor(); + let mut c = Cache::new(path.n()); + let mut re = Vec::::new(); + let mut im = Vec::::new(); + if is_singlet { + let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd( + args.order_qcd, + &mut c, + args.nf, + ); + for k in 0..args.order_qcd { + for l in 0..2 { + for m in 0..2 { + re.push(res[k][l][m].re); + im.push(res[k][l][m].im); + } + } + } + } else { + let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd( + args.order_qcd, + args.mode0, + &mut c, + args.nf, + ); + for j in 0..args.order_qcd { + re.push(res[j].re); + im.push(res[j].im); + } + } + // pass on + (args.py)( + re.as_ptr(), + im.as_ptr(), + c.n.re, + c.n.im, + jac.re, + jac.im, + args.order_qcd, + is_singlet, + args.mode0, + args.mode1, + args.nf, + args.is_log, + args.logx, + args.areas, + args.areas_x, + args.areas_y, + args.L, + args.method_num, + args.as1, + args.as0, + args.ev_op_iterations, + args.ev_op_max_order_qcd, + args.sv_mode_num, + args.is_threshold, + ) +} + +/// Python callback signature +type PyQuadKerQCDT = unsafe extern "C" fn( + *const f64, + *const f64, + f64, + f64, + f64, + f64, + usize, + bool, + u16, + u16, + u8, + bool, + f64, + *const f64, + u8, + u8, + f64, + u8, + f64, + f64, + u8, + u8, + u8, + bool, +) -> f64; + +/// Additional integration parameters +#[allow(non_snake_case)] +#[repr(C)] +#[derive(Clone, Copy)] +pub struct QuadQCDargs { + pub order_qcd: usize, + pub mode0: u16, + pub mode1: u16, + pub is_polarized: bool, + pub is_time_like: bool, + pub nf: u8, + pub py: PyQuadKerQCDT, + pub is_log: bool, + pub logx: f64, + pub areas: *const f64, + pub areas_x: u8, + pub areas_y: u8, + pub L: f64, + pub method_num: u8, + pub as1: f64, + pub as0: f64, + pub ev_op_iterations: u8, + pub ev_op_max_order_qcd: u8, + pub sv_mode_num: u8, + pub is_threshold: bool, +} + +/// empty placeholder function for python callback +pub unsafe extern "C" fn my_py( + _re_gamma: *const f64, + _im_gamma: *const f64, + _re_n: f64, + _im_n: f64, + _re_jac: f64, + _im_jac: f64, + _order_qcd: usize, + _is_singlet: bool, + _mode0: u16, + _mode1: u16, + _nf: u8, + _is_log: bool, + _logx: f64, + _areas: *const f64, + _areas_x: u8, + _areas_y: u8, + _l: f64, + _method_num: u8, + _as1: f64, + _as0: f64, + _ev_op_iterations: u8, + _ev_op_max_order_qcd: u8, + _sv_mode_num: u8, + _is_threshold: bool, +) -> f64 { + 0. +} + +/// Return empty additional arguments. +/// +/// This is required to make the arguments part of the API, otherwise it won't be added to the compiled +/// package (since it does not appear in the signature of `quad_ker_qcd`). +#[no_mangle] +pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { + QuadQCDargs { + order_qcd: 0, + mode0: 0, + mode1: 0, + is_polarized: false, + is_time_like: false, + nf: 0, + py: my_py, + is_log: true, + logx: 0., + areas: [].as_ptr(), + areas_x: 0, + areas_y: 0, + L: 0., + method_num: 0, + as1: 0., + as0: 0., + ev_op_iterations: 0, + ev_op_max_order_qcd: 0, + sv_mode_num: 0, + is_threshold: false, + } +} diff --git a/crates/eko/src/mellin.rs b/crates/eko/src/mellin.rs new file mode 100644 index 000000000..5afc96f5f --- /dev/null +++ b/crates/eko/src/mellin.rs @@ -0,0 +1,77 @@ +//! Tools for Mellin inversion. +//! +//! We provide all necessary toold to deal with the +//! [inverse Mellin transformation](https://en.wikipedia.org/wiki/Mellin_inversion_theorem). + +use num::complex::Complex; +use std::f64::consts::PI; + +#[cfg_attr(doc, katexit::katexit)] +/// Talbot inversion path. +/// +/// Implements the algorithm presented in [\[Abate\]](crate::bib::Abate). +/// $p_{\text{Talbot}}(t) = o + r \cdot ( \theta \cot(\theta) + i\theta)$ with $\theta = \pi(2t-1)$ +/// The default values for the parameters $r,o$ are given by $r = 1/2, o = 0$ for +/// the non-singlet integrals and by $r = \frac{2}{5} \frac{16}{1 - \ln(x)}, o = 1$ +/// for the singlet sector. Note that the non-singlet kernels evolve poles only up to +/// $N=0$ whereas the singlet kernels have poles up to $N=1$. +pub struct TalbotPath { + /// integration variable + t: f64, + + /// bending variable + r: f64, + + /// real offset + o: f64, +} + +impl TalbotPath { + /// Auxilary angle. + fn theta(&self) -> f64 { + PI * (2.0 * self.t - 1.0) + } + + /// Constructor from parameters. + pub fn new(t: f64, logx: f64, is_singlet: bool) -> Self { + Self { + t, + r: 0.4 * 16.0 / (1.0 - logx), + o: if is_singlet { 1. } else { 0. }, + } + } + + /// Mellin N. + pub fn n(&self) -> Complex { + let theta = self.theta(); + // treat singular point separately + let re = if 0.5 == self.t { + 1.0 + } else { + theta / theta.tan() + }; + let im = theta; + self.o + self.r * Complex:: { re, im } + } + + /// Transformation jacobian. + pub fn jac(&self) -> Complex { + let theta = self.theta(); + // treat singular point separately + let re = if 0.5 == self.t { + 0.0 + } else { + 1.0 / theta.tan() - theta / theta.sin().powf(2.) + }; + let im = 1.0; + self.r * PI * 2.0 * Complex:: { re, im } + } + + /// Mellin inversion prefactor. + pub fn prefactor(&self) -> Complex { + Complex:: { + re: 0., + im: -1. / PI, + } + } +} diff --git a/crates/ekore/Cargo.toml b/crates/ekore/Cargo.toml new file mode 100644 index 000000000..ef8499dba --- /dev/null +++ b/crates/ekore/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "ekore" +version = "0.1.0" +edition = "2021" +description = "EKO expressions" +license = "GPL-3.0-or-later" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +num = "0.4.1" +float-cmp = "0.9.0" +katexit = "0.1.4" +hashbrown = "0.14" diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib new file mode 100644 index 000000000..3c63b69e7 --- /dev/null +++ b/crates/ekore/refs.bib @@ -0,0 +1,52 @@ +@article{Moch2004pa, + author = "Moch, S. and Vermaseren, J. A. M. and Vogt, A.", + title = "{The Three loop splitting functions in QCD: The + Nonsinglet case}", + journal = "Nucl. Phys.", + volume = "B688", + year = "2004", + pages = "101-134", + doi = "10.1016/j.nuclphysb.2004.03.030", + eprint = "hep-ph/0403192", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "DESY-04-047, SFB-CPP-04-09, NIKHEF-04-001", + SLACcitation = "%%CITATION = HEP-PH/0403192;%%" +} +@article{Vogt2004mw, + author = "Vogt, A. and Moch, S. and Vermaseren, J. A. M.", + title = "{The Three-loop splitting functions in QCD: The Singlet + case}", + journal = "Nucl. Phys.", + volume = "B691", + year = "2004", + pages = "129-181", + doi = "10.1016/j.nuclphysb.2004.04.024", + eprint = "hep-ph/0404111", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "NIKHEF-04-004, DESY-04-060, SFB-CPP-04-12", + SLACcitation = "%%CITATION = HEP-PH/0404111;%%" +} +@article{KOLBIG1972221, + title = "Programs for computing the logarithm of the gamma function, and the digamma function, for complex argument", + journal = "Computer Physics Communications", + volume = "4", + number = "2", + pages = "221 - 226", + year = "1972", + issn = "0010-4655", + doi = "https://doi.org/10.1016/0010-4655(72)90012-4", + url = "http://www.sciencedirect.com/science/article/pii/0010465572900124", + author = "K.S. Kölbig" +} +@article{Abate, + author = {Abate, J. and Valkó, P.}, + year = {2004}, + month = {06}, + pages = {979 - 993}, + title = {Multi‐precision Laplace transform inversion}, + volume = {60}, + journal = {International Journal for Numerical Methods in Engineering}, + doi = {10.1002/nme.995} +} diff --git a/crates/ekore/src/anomalous_dimensions.rs b/crates/ekore/src/anomalous_dimensions.rs new file mode 100644 index 000000000..6a395d982 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions.rs @@ -0,0 +1,3 @@ +//! The anomalous dimensions for DGLAP evolution. + +pub mod unpolarized; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized.rs b/crates/ekore/src/anomalous_dimensions/unpolarized.rs new file mode 100644 index 000000000..19bb654d2 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized.rs @@ -0,0 +1,3 @@ +//! The unpolarized anomalous dimensions for space-like and time-like kinematics. + +pub mod spacelike; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs new file mode 100644 index 000000000..c2d0fb0d7 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -0,0 +1,26 @@ +//! The unpolarized, space-like anomalous dimensions at various couplings power. + +use crate::harmonics::cache::Cache; +use num::complex::Complex; +use num::Zero; +pub mod as1; + +/// Compute the tower of the non-singlet anomalous dimensions. +pub fn gamma_ns_qcd(order_qcd: usize, _mode: u16, c: &mut Cache, nf: u8) -> Vec> { + let mut gamma_ns = vec![Complex::::zero(); order_qcd]; + gamma_ns[0] = as1::gamma_ns(c, nf); + gamma_ns +} + +/// Compute the tower of the singlet anomalous dimension matrices. +pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Complex; 2]; 2]> { + let mut gamma_S = vec![ + [ + [Complex::::zero(), Complex::::zero()], + [Complex::::zero(), Complex::::zero()] + ]; + order_qcd + ]; + gamma_S[0] = as1::gamma_singlet(c, nf); + gamma_S +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs new file mode 100644 index 000000000..71e285ea9 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -0,0 +1,116 @@ +//! LO QCD + +use num::complex::Complex; + +use crate::constants::{CA, CF, TR}; +use crate::harmonics::cache::{Cache, K}; + +/// Compute the non-singlet anomalous dimension. +/// +/// Implements Eq. (3.4) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. +pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n; + let S1 = c.get(K::S1); + let gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0)); + CF * gamma +} + +/// Compute the quark-gluon anomalous dimension. +/// +/// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw). +pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { + let N = c.n; + let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N + 2.0)); + 2.0 * TR * 2.0 * (nf as f64) * gamma +} + +/// Compute the gluon-quark anomalous dimension. +/// +/// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw). +pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n; + let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N - 1.0)); + 2.0 * CF * gamma +} + +/// Compute the gluon-gluon anomalous dimension. +/// +/// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw). +pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { + let N = c.n; + let S1 = c.get(K::S1); + let gamma = S1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0); + CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * (nf as f64) +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let gamma_qq = gamma_ns(c, nf); + [ + [gamma_qq, gamma_qg(c, nf)], + [gamma_gq(c, nf), gamma_gg(c, nf)], + ] +} + +#[cfg(test)] +mod tests { + use crate::cmplx; + use crate::{anomalous_dimensions::unpolarized::spacelike::as1::*, harmonics::cache::Cache}; + use float_cmp::assert_approx_eq; + use num::complex::Complex; + const NF: u8 = 5; + + #[test] + fn number_conservation() { + const N: Complex = cmplx![1., 0.]; + let mut c = Cache::new(N); + let me = gamma_ns(&mut c, NF); + assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12); + assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + } + + #[test] + fn quark_momentum_conservation() { + const N: Complex = cmplx![2., 0.]; + let mut c = Cache::new(N); + let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF); + assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12); + assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + } + + #[test] + fn gluon_momentum_conservation() { + const N: Complex = cmplx![2., 0.]; + let mut c = Cache::new(N); + let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF); + assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12); + assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + } + + #[test] + fn gamma_qg_() { + const N: Complex = cmplx![1., 0.]; + let mut c = Cache::new(N); + let me = gamma_qg(&mut c, NF); + assert_approx_eq!(f64, me.re, -20. / 3.0, ulps = 32); + assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + } + + #[test] + fn gamma_gq_() { + const N: Complex = cmplx![0., 1.]; + let mut c = Cache::new(N); + let me = gamma_gq(&mut c, NF); + assert_approx_eq!(f64, me.re, 4. / 3.0, ulps = 32); + assert_approx_eq!(f64, me.im, -4. / 3.0, ulps = 32); + } + + #[test] + fn gamma_gg_() { + const N: Complex = cmplx![0., 1.]; + let mut c = Cache::new(N); + let me = gamma_gg(&mut c, NF); + assert_approx_eq!(f64, me.re, 5.195725159621, ulps = 32, epsilon = 1e-11); + assert_approx_eq!(f64, me.im, 10.52008856962, ulps = 32, epsilon = 1e-11); + } +} diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs new file mode 100644 index 000000000..c62f16d02 --- /dev/null +++ b/crates/ekore/src/bib.rs @@ -0,0 +1,45 @@ +//! List of References (autogenerated on 2023-11-29T18:07:41.518255). + +/// The Three loop splitting functions in QCD: The Nonsinglet case +/// +/// Moch, S. and Vermaseren, J. A. M. and Vogt, A. +/// +/// Published in: Nucl. Phys. B688 (2004), 101-134 +/// +/// e-Print: [hep-ph/0403192](https://arxiv.org/abs/hep-ph/0403192) +/// +/// DOI: [10.1016/j.nuclphysb.2004.03.030](https:dx.doi.org/10.1016/j.nuclphysb.2004.03.030) +pub fn Moch2004pa() {} + +/// The Three-loop splitting functions in QCD: The Singlet case +/// +/// Vogt, A. and Moch, S. and Vermaseren, J. A. M. +/// +/// Published in: Nucl. Phys. B691 (2004), 129-181 +/// +/// e-Print: [hep-ph/0404111](https://arxiv.org/abs/hep-ph/0404111) +/// +/// DOI: [10.1016/j.nuclphysb.2004.04.024](https:dx.doi.org/10.1016/j.nuclphysb.2004.04.024) +pub fn Vogt2004mw() {} + +/// Programs for computing the logarithm of the gamma function, and the digamma function, for complex argument +/// +/// K.S. Kölbig +/// +/// Published in: Computer Physics Communications 4 (1972), 221 - 226 +/// +/// +/// +/// DOI: [https://doi.org/10.1016/0010-4655(72)90012-4](https://doi.org/10.1016/0010-4655(72)90012-4) +pub fn KOLBIG1972221() {} + +/// Multi‐precision Laplace transform inversion +/// +/// Abate, J. and Valkó, P. +/// +/// Published in: International Journal for Numerical Methods in Engineering 60 (2004), 979 - 993 +/// +/// +/// +/// DOI: [10.1002/nme.995](https:dx.doi.org/10.1002/nme.995) +pub fn Abate() {} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs new file mode 100644 index 000000000..3cdd13a15 --- /dev/null +++ b/crates/ekore/src/constants.rs @@ -0,0 +1,25 @@ +//! Global constants. + +#[cfg_attr(doc, katexit::katexit)] +/// The number of colors. +/// +/// Defaults to $N_C = 3$. +pub const NC: u8 = 3; + +#[cfg_attr(doc, katexit::katexit)] +/// The normalization of fundamental generators. +/// +/// Defaults to $T_R = 1/2$. +pub const TR: f64 = 1.0 / 2.0; + +#[cfg_attr(doc, katexit::katexit)] +/// Second Casimir constant in the adjoint representation. +/// +/// Defaults to $C_A = N_C = 3$. +pub const CA: f64 = NC as f64; + +#[cfg_attr(doc, katexit::katexit)] +/// Second Casimir constant in the fundamental representation. +/// +/// Defaults to $C_F = \frac{N_C^2-1}{2N_C} = 4/3$. +pub const CF: f64 = ((NC * NC - 1) as f64) / ((2 * NC) as f64); diff --git a/crates/ekore/src/harmonics.rs b/crates/ekore/src/harmonics.rs new file mode 100644 index 000000000..b829060fc --- /dev/null +++ b/crates/ekore/src/harmonics.rs @@ -0,0 +1,5 @@ +//! Tools to compute harmonic sums and related special functions. + +pub mod cache; +pub mod polygamma; +mod w1; diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs new file mode 100644 index 000000000..6ad72171d --- /dev/null +++ b/crates/ekore/src/harmonics/cache.rs @@ -0,0 +1,48 @@ +//! Cache harmonic sums for given Mellin N. + +use hashbrown::HashMap; +use num::complex::Complex; + +use crate::harmonics::w1; + +#[cfg_attr(doc, katexit::katexit)] +/// List of available elements. +#[derive(Debug, PartialEq, Eq, Hash)] +pub enum K { + /// $S_1(N)$ + S1, +} + +/// Hold all cached values. +pub struct Cache { + /// Mellin N + pub n: Complex, + /// Mapping + m: HashMap>, +} + +impl Cache { + /// Initialize new, empty Cache at given Mellin N. + pub fn new(n: Complex) -> Self { + Self { + n, + m: HashMap::new(), + } + } + + /// Retrieve an element. + pub fn get(&mut self, k: K) -> Complex { + let val = self.m.get(&k); + // already there? + if val.is_some() { + return *val.unwrap(); + } + // compute new + let val = match k { + K::S1 => w1::S1(self.n), + }; + // insert + self.m.insert(k, val); + val + } +} diff --git a/crates/ekore/src/harmonics/polygamma.rs b/crates/ekore/src/harmonics/polygamma.rs new file mode 100644 index 000000000..efbdae147 --- /dev/null +++ b/crates/ekore/src/harmonics/polygamma.rs @@ -0,0 +1,210 @@ +//! Tools to compute [polygamma functions](https://en.wikipedia.org/wiki/Polygamma_function). + +use num::{complex::Complex, Zero}; +use std::f64::consts::PI; + +#[cfg_attr(doc, katexit::katexit)] +/// Compute the polygamma functions $\psi_k(z)$. +/// +/// Reimplementation of ``WPSIPG`` (C317) in [CERNlib](http://cernlib.web.cern.ch/cernlib/) given by [[KOLBIG1972221]][crate::bib::KOLBIG1972221]. +/// +/// TODO: introduce back errors +pub fn cern_polygamma(Z: Complex, K: usize) -> Complex { + // const DELTA: f64 = 5e-13; + const R1: f64 = 1.0; + const HF: f64 = R1 / 2.0; + let C1: f64 = PI.powi(2); + let C2: f64 = 2. * PI.powi(3); + let C3: f64 = 2. * PI.powi(4); + let C4: f64 = 8. * PI.powi(5); + + // SGN is originally indexed 0:4 -> no shift + const SGN: [i8; 5] = [-1, 1, -1, 1, -1]; + // FCT is originally indexed -1:4 -> shift +1 + const FCT: [u8; 6] = [0, 1, 1, 2, 6, 24]; + + // C is originally indexed 1:6 x 0:4 -> swap indices and shift new last -1 + const C: [[f64; 6]; 5] = [ + [ + 8.33333333333333333e-2, + -8.33333333333333333e-3, + 3.96825396825396825e-3, + -4.16666666666666667e-3, + 7.57575757575757576e-3, + -2.10927960927960928e-2, + ], + [ + 1.66666666666666667e-1, + -3.33333333333333333e-2, + 2.38095238095238095e-2, + -3.33333333333333333e-2, + 7.57575757575757576e-2, + -2.53113553113553114e-1, + ], + [ + 5.00000000000000000e-1, + -1.66666666666666667e-1, + 1.66666666666666667e-1, + -3.00000000000000000e-1, + 8.33333333333333333e-1, + -3.29047619047619048e+0, + ], + [ + 2.00000000000000000e+0, + -1.00000000000000000e+0, + 1.33333333333333333e+0, + -3.00000000000000000e+0, + 1.00000000000000000e+1, + -4.60666666666666667e+1, + ], + [10., -7., 12., -33., 130., -691.], + ]; + let mut U = Z; + let mut X = U.re; + let A = X.abs(); + // if (K < 0 || K > 4) + // throw InvalidPolygammaOrder("Order K has to be in [0:4]"); + let A_as_int = A as i64; + // if (fabs(U.imag()) < DELTA && fabs(X+A_as_int) < DELTA) + // throw InvalidPolygammaArgument("Argument Z equals non-positive integer"); + let K1 = K + 1; + if X < 0. { + U = -U; + } + let mut V = U; + let mut H = Complex::zero(); + if A < 15. { + H = 1. / V.powu(K1 as u32); + for _ in 1..(14 - A_as_int + 1) { + V = V + 1.; + H = H + 1. / V.powu(K1 as u32); + } + V = V + 1.; + } + let mut R = 1. / V.powu(2); + let mut P = R * C[K][6 - 1]; + for i in (1..=5).rev() + // (int i = 5; i>1-1; i--) + { + P = R * (C[K][i - 1] + P); + } + H = (SGN[K] as f64) + * ((FCT[K + 1] as f64) * H + + (V * ((FCT[K] as f64) + P) + HF * (FCT[K + 1] as f64)) / V.powu(K1 as u32)); + if 0 == K { + H = H + V.ln(); + } + if X < 0. { + V = PI * U; + X = V.re; + let Y = V.im; + let A = X.sin(); + let B = X.cos(); + let T = Y.tanh(); + P = Complex:: { re: B, im: -A * T } / Complex:: { re: A, im: B * T }; + if 0 == K { + H = H + 1. / U + PI * P; + } else if 1 == K { + H = -H + 1. / U.powu(2) + C1 * (P.powu(2) + 1.); + } else if 2 == K { + H = H + 2. / U.powu(3) + C2 * P * (P.powu(2) + 1.); + } else if 3 == K { + R = P.powu(2); + H = -H + 6. / U.powu(4) + C3 * ((3. * R + 4.) * R + 1.); + } else if 4 == K { + R = P.powu(2); + H = H + 24. / U.powu(5) + C4 * P * ((3. * R + 5.) * R + 2.); + } + } + H +} + +#[cfg(test)] +mod tests { + use crate::cmplx; + use crate::harmonics::polygamma::cern_polygamma; + use float_cmp::assert_approx_eq; + use num::complex::Complex; + + #[test] + fn fortran() { + const ZS: [Complex; 9] = [ + cmplx![1.0, 0.], + cmplx![2.0, 0.], + cmplx![3.0, 0.], + cmplx![0., 1.], + cmplx![-1., 1.], + cmplx![-2., 1.], + cmplx![-1., 2.], + cmplx![-2., 2.], + cmplx![-3., 2.], + ]; + const KS: [usize; 5] = [0, 1, 2, 3, 4]; + + const FORTRAN_REF: [[Complex; 9]; 5] = [ + [ + cmplx![-0.5772156649015332, 0.], + cmplx![0.42278433509846636, 0.], + cmplx![0.9227843350984672, 0.], + cmplx![0.09465032062247669, 2.0766740474685816], + cmplx![0.5946503206224767, 2.5766740474685816], + cmplx![0.9946503206224772, 2.7766740474685814], + cmplx![0.9145915153739776, 2.2208072826422303], + cmplx![1.1645915153739772, 2.47080728264223], + cmplx![1.395360746143208, 2.624653436488384], + ], + [ + cmplx![1.6449340668482264, 0.], + cmplx![0.6449340668482264, 0.], + cmplx![0.3949340668482264, 0.], + cmplx![-0.5369999033772361, -0.7942335427593189], + cmplx![-0.5369999033772362, -0.2942335427593189], + cmplx![-0.4169999033772362, -0.13423354275931887], + cmplx![-0.24506883785905695, -0.3178255501472297], + cmplx![-0.24506883785905695, -0.19282555014722969], + cmplx![-0.21548303904248892, -0.12181963298746643], + ], + [ + cmplx![-2.404113806319188, 0.], + cmplx![-0.40411380631918853, 0.], + cmplx![-0.15411380631918858, 0.], + cmplx![0.3685529315879351, -1.233347149654934], + cmplx![-0.13144706841206477, -0.7333471496549337], + cmplx![-0.09944706841206462, -0.5573471496549336], + cmplx![0.03902435405364951, -0.15743252404131272], + cmplx![-0.02347564594635048, -0.09493252404131272], + cmplx![-0.031668636387861625, -0.053057239562477945], + ], + [ + cmplx![6.49393940226683, 0.], + cmplx![0.49393940226682925, 0.], + cmplx![0.11893940226682913, 0.], + cmplx![4.4771255510465044, -0.31728657866196064], + cmplx![2.9771255510464925, -0.3172865786619599], + cmplx![2.909925551046492, -0.08688657866195917], + cmplx![0.12301766661068443, -0.05523068481179527], + cmplx![0.02926766661068438, -0.055230684811795216], + cmplx![0.004268541930176011, -0.03002148345329936], + ], + [ + cmplx![-24.88626612344089, 0.], + cmplx![-0.8862661234408784, 0.], + cmplx![-0.13626612344087824, 0.], + cmplx![3.2795081690440493, 21.41938794863803], + cmplx![0.2795081690440445, 18.419387948637894], + cmplx![-0.012331830955960252, 18.734267948637896], + cmplx![0.14223316576854003, 0.10023607930398608], + cmplx![0.04848316576854002, 0.006486079303986134], + cmplx![0.009893695996688708, 0.014372034600746361], + ], + ]; + for kit in KS.iter().enumerate() { + for zit in ZS.iter().enumerate() { + let fref = FORTRAN_REF[kit.0][zit.0]; + let me = cern_polygamma(*zit.1, *kit.1); + assert_approx_eq!(f64, me.re, fref.re, ulps = 32); + assert_approx_eq!(f64, me.im, fref.im, ulps = 32); + } + } + } +} diff --git a/crates/ekore/src/harmonics/w1.rs b/crates/ekore/src/harmonics/w1.rs new file mode 100644 index 000000000..2b861cb30 --- /dev/null +++ b/crates/ekore/src/harmonics/w1.rs @@ -0,0 +1,11 @@ +use num::complex::Complex; + +use crate::harmonics::polygamma::cern_polygamma; + +/// Compute the harmonic sum $S_1(N)$. +/// +/// $$S_1(N) = \sum\limits_{j=1}^N \frac 1 j = \psi_0(N+1)+\gamma_E$$ +/// with $\psi_0(N)$ the digamma function and $\gamma_E$ the Euler-Mascheroni constant. +pub fn S1(N: Complex) -> Complex { + cern_polygamma(N + 1.0, 0) + 0.5772156649015328606065120 +} diff --git a/crates/ekore/src/lib.rs b/crates/ekore/src/lib.rs new file mode 100644 index 000000000..af39644cc --- /dev/null +++ b/crates/ekore/src/lib.rs @@ -0,0 +1,10 @@ +//! Library for anomalous dimension in DGLAP and transition matrix elements. + +// Let's stick to the original names which often come from FORTRAN, where such convention do not exists +#![allow(non_snake_case)] + +pub mod anomalous_dimensions; +pub mod bib; +mod constants; +pub mod harmonics; +pub mod util; diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs new file mode 100644 index 000000000..5c4e34f4d --- /dev/null +++ b/crates/ekore/src/util.rs @@ -0,0 +1,9 @@ +//! Helper utilities. + +/// Shorthand complex number contructor. +#[macro_export] +macro_rules! cmplx { + ($re:expr, $im:expr) => { + Complex::new($re, $im) + }; +} diff --git a/crates/make_bib.py b/crates/make_bib.py new file mode 100644 index 000000000..1b15d663d --- /dev/null +++ b/crates/make_bib.py @@ -0,0 +1,67 @@ +"""Parse bibtex file to rust crate.""" +import datetime +import pathlib +import re + +import bibtexparser + +# path to bib file +BIBFILE = pathlib.Path(__file__).parent / "ekore" / "refs.bib" + +# A single entry +ENTRY = """/// {title} +/// +/// {author} +/// +/// {publication} +/// +/// {eprint} +/// +/// {doi}""" +# Combine publication information +PUB = """Published in: {journal} {volume} ({year}), {pages}""" + + +def clean_nl(t: str) -> str: + """Shrink whitespace to one.""" + return re.sub("\n\\s+", " ", t) + + +# collect output +out = "//! List of References (autogenerated on {now}).\n\n".format( + now=datetime.datetime.now().isoformat() +) +# Parse .bib file +bib_str = pathlib.Path(BIBFILE).read_text("utf8") +bib_database = bibtexparser.parse_string(bib_str) +# iterate on all elements +for el in bib_database.entries: + title = re.sub(r"^\{(.+)\}$", r"\1", clean_nl(el.fields_dict["title"].value)) + author = el.fields_dict["author"].value + publication = PUB.format( + journal=el.fields_dict["journal"].value, + volume=el.fields_dict["volume"].value, + year=el.fields_dict["year"].value, + pages=el.fields_dict["pages"].value, + ) + eprint = "" + if ( + "eprint" in el.fields_dict + and "archivePrefix" in el.fields_dict + and el.fields_dict["archivePrefix"].value == "arXiv" + ): + ar = el.fields_dict["eprint"].value + eprint = f"e-Print: [{ar}](https://arxiv.org/abs/{ar})" + doi = "" + if "doi" in el.fields_dict: + ddoi = el.fields_dict["doi"].value + url = ddoi if "http" in ddoi else f"https:dx.doi.org/{ddoi}" + doi = f"DOI: [{ddoi}]({url})" + doc = ENTRY.format( + title=title, author=author, publication=publication, eprint=eprint, doi=doi + ) + out += doc + "\n" + f"pub fn {el.key}() {{}} \n\n" + +# Print to terminal - the user can dump to the relevant file +out = re.sub(r" +\n", "\n", out.strip()) +print(out) diff --git a/poetry.lock b/poetry.lock index 3ca57fb0e..997523b2d 100644 --- a/poetry.lock +++ b/poetry.lock @@ -661,13 +661,13 @@ files = [ [[package]] name = "exceptiongroup" -version = "1.1.3" +version = "1.1.1" description = "Backport of PEP 654 (exception groups)" optional = false python-versions = ">=3.7" files = [ - {file = "exceptiongroup-1.1.3-py3-none-any.whl", hash = "sha256:343280667a4585d195ca1cf9cef84a4e178c4b6cf2274caef9859782b567d5e3"}, - {file = "exceptiongroup-1.1.3.tar.gz", hash = "sha256:097acd85d473d75af5bb98e41b61ff7fe35efe6675e4f9370ec6ec5126d160e9"}, + {file = "exceptiongroup-1.1.1-py3-none-any.whl", hash = "sha256:232c37c63e4f682982c8b6459f33a8981039e5fb8756b2074364e5055c498c9e"}, + {file = "exceptiongroup-1.1.1.tar.gz", hash = "sha256:d484c3090ba2889ae2928419117447a14daf3c1231d5e30d0aae34f354f01785"}, ] [package.extras] @@ -945,13 +945,13 @@ files = [ [[package]] name = "ipykernel" -version = "6.25.1" +version = "6.23.2" description = "IPython Kernel for Jupyter" optional = false python-versions = ">=3.8" files = [ - {file = "ipykernel-6.25.1-py3-none-any.whl", hash = "sha256:c8a2430b357073b37c76c21c52184db42f6b4b0e438e1eb7df3c4440d120497c"}, - {file = "ipykernel-6.25.1.tar.gz", hash = "sha256:050391364c0977e768e354bdb60cbbfbee7cbb943b1af1618382021136ffd42f"}, + {file = "ipykernel-6.23.2-py3-none-any.whl", hash = "sha256:7ccb6e2d32fd958c21453db494c914f3474908a2fdefd99ab548a5375b548d1f"}, + {file = "ipykernel-6.23.2.tar.gz", hash = "sha256:fcfb67c5b504aa1bfcda1c5b3716636239e0f7b9290958f1c558c79b4c0e7ed5"}, ] [package.dependencies] @@ -1569,13 +1569,13 @@ test = ["flaky", "ipykernel (>=6.19.3)", "ipython", "ipywidgets", "nbconvert (>= [[package]] name = "nbconvert" -version = "7.7.4" +version = "7.5.0" description = "Converting Jupyter Notebooks" optional = false python-versions = ">=3.8" files = [ - {file = "nbconvert-7.7.4-py3-none-any.whl", hash = "sha256:ace26f4386d08eb5c55833596a942048c5502a95e05590cb523826a749a40a37"}, - {file = "nbconvert-7.7.4.tar.gz", hash = "sha256:1113d039fa3fc3a846ffa5a3b0a019e85aaa94c566a09fa0c400fb7638e46087"}, + {file = "nbconvert-7.5.0-py3-none-any.whl", hash = "sha256:852e44392d5650ef217a5ce3a8050747051d4e6ba75f0574cb5435049ee6c0d9"}, + {file = "nbconvert-7.5.0.tar.gz", hash = "sha256:f78fd22fd2410b960d5d9bcecf3e1d6c7bdc5fec2c865964c84aa4e74e6e88da"}, ] [package.dependencies] @@ -3078,13 +3078,13 @@ zstd = ["zstandard (>=0.18.0)"] [[package]] name = "virtualenv" -version = "20.24.3" +version = "20.23.0" description = "Virtual Python Environment builder" optional = false python-versions = ">=3.7" files = [ - {file = "virtualenv-20.24.3-py3-none-any.whl", hash = "sha256:95a6e9398b4967fbcb5fef2acec5efaf9aa4972049d9ae41f95e0972a683fd02"}, - {file = "virtualenv-20.24.3.tar.gz", hash = "sha256:e5c3b4ce817b0b328af041506a2a299418c98747c4b1e68cb7527e74ced23efc"}, + {file = "virtualenv-20.23.0-py3-none-any.whl", hash = "sha256:6abec7670e5802a528357fdc75b26b9f57d5d92f29c5462ba0fbe45feacc685e"}, + {file = "virtualenv-20.23.0.tar.gz", hash = "sha256:a85caa554ced0c0afbd0d638e7e2d7b5f92d23478d05d17a76daeac8f279f924"}, ] [package.dependencies] diff --git a/pyproject.toml.patch b/pyproject.toml.patch new file mode 100644 index 000000000..65b5daa2c --- /dev/null +++ b/pyproject.toml.patch @@ -0,0 +1,39 @@ +diff --git a/pyproject.toml b/pyproject.toml +index 31be6cb0..b4ec7c95 100644 +--- a/pyproject.toml ++++ b/pyproject.toml +@@ -1,6 +1,20 @@ + [build-system] +-requires = ["poetry-core>=1.0.0", "poetry-dynamic-versioning"] +-build-backend = "poetry.core.masonry.api" ++requires = ["maturin>=0.14.8,<0.15", "poetry-dynamic-versioning"] ++build-backend = "maturin" ++ ++[project] ++name = "eko" ++requires-python = ">=3.8" ++classifiers = [ ++ "Programming Language :: Rust", ++ "Programming Language :: Python :: Implementation :: CPython", ++ "Programming Language :: Python :: Implementation :: PyPy", ++] ++version = "0.0.0" ++ ++[tool.maturin] ++manifest-path = "crates/eko/Cargo.toml" ++python-packages = ["ekomark", "ekobox"] + + [tool.poetry] + name = "eko" +@@ -124,6 +138,11 @@ asv-publish = "asv publish --config benchmarks/asv.conf.json" + asv-show = "asv show --config benchmarks/asv.conf.json" + asv-clean = { "shell" = "rm -rf benchmarks/env benchmarks/html benchmarks/results" } + asv = ["asv-run", "asv-publish", "asv-preview"] ++compile = "pip install -e crates/eko/" ++rdocs = "cargo doc --workspace --manifest-path crates/Cargo.toml --no-deps" ++rdocs-view = "xdg-open crates/target/doc/ekors/index.html" ++rdocs-clean = "rm -rf crates/target/doc/" ++rtest = "cargo test --workspace --manifest-path crates/Cargo.toml" + + [tool.pytest.ini_options] + testpaths = ['tests/', 'benchmarks/'] diff --git a/rustify.sh b/rustify.sh new file mode 100755 index 000000000..05078d0b7 --- /dev/null +++ b/rustify.sh @@ -0,0 +1,5 @@ +#!/usr/bin/bash + +patch -p1 `. + """ + +-import functools + import logging + import os + import time + from multiprocessing import Pool + ++import ekors + import numba as nb + import numpy as np +-from scipy import integrate ++from scipy import LowLevelCallable, integrate + + import ekore.anomalous_dimensions.polarized.space_like as ad_ps + import ekore.anomalous_dimensions.unpolarized.space_like as ad_us +@@ -28,92 +28,10 @@ from ..kernels import utils + from ..kernels import valence_qed as qed_v + from ..matchings import Segment + from ..member import OpMember ++from .quad_ker import cb_quad_ker_qcd + + logger = logging.getLogger(__name__) + +- +-@nb.njit(cache=True) +-def select_singlet_element(ker, mode0, mode1): +- """Select element of the singlet matrix. +- +- Parameters +- ---------- +- ker : numpy.ndarray +- singlet integration kernel +- mode0 : int +- id for first sector element +- mode1 : int +- id for second sector element +- +- Returns +- ------- +- complex +- singlet integration kernel element +- """ +- k = 0 if mode0 == 100 else 1 +- l = 0 if mode1 == 100 else 1 +- return ker[k, l] +- +- +-@nb.njit(cache=True) +-def select_QEDsinglet_element(ker, mode0, mode1): +- """Select element of the QEDsinglet matrix. +- +- Parameters +- ---------- +- ker : numpy.ndarray +- QEDsinglet integration kernel +- mode0 : int +- id for first sector element +- mode1 : int +- id for second sector element +- Returns +- ------- +- ker : complex +- QEDsinglet integration kernel element +- """ +- if mode0 == 21: +- index1 = 0 +- elif mode0 == 22: +- index1 = 1 +- elif mode0 == 100: +- index1 = 2 +- else: +- index1 = 3 +- if mode1 == 21: +- index2 = 0 +- elif mode1 == 22: +- index2 = 1 +- elif mode1 == 100: +- index2 = 2 +- else: +- index2 = 3 +- return ker[index1, index2] +- +- +-@nb.njit(cache=True) +-def select_QEDvalence_element(ker, mode0, mode1): +- """ +- Select element of the QEDvalence matrix. +- +- Parameters +- ---------- +- ker : numpy.ndarray +- QEDvalence integration kernel +- mode0 : int +- id for first sector element +- mode1 : int +- id for second sector element +- Returns +- ------- +- ker : complex +- QEDvalence integration kernel element +- """ +- index1 = 0 if mode0 == 10200 else 1 +- index2 = 0 if mode1 == 10200 else 1 +- return ker[index1, index2] +- +- + spec = [ + ("is_singlet", nb.boolean), + ("is_QEDsinglet", nb.boolean), +@@ -185,403 +103,6 @@ class QuadKerBase: + return self.path.prefactor * pj * self.path.jac + + +-@nb.njit(cache=True) +-def quad_ker( +- u, +- order, +- mode0, +- mode1, +- method, +- is_log, +- logx, +- areas, +- as_list, +- mu2_from, +- mu2_to, +- a_half, +- alphaem_running, +- nf, +- L, +- ev_op_iterations, +- ev_op_max_order, +- sv_mode, +- is_threshold, +- n3lo_ad_variation, +- is_polarized, +- is_time_like, +-): +- """Raw evolution kernel inside quad. +- +- Parameters +- ---------- +- u : float +- quad argument +- order : int +- perturbation order +- mode0: int +- pid for first sector element +- mode1 : int +- pid for second sector element +- method : str +- method +- is_log : boolean +- is a logarithmic interpolation +- logx : float +- Mellin inversion point +- areas : tuple +- basis function configuration +- as1 : float +- target coupling value +- as0 : float +- initial coupling value +- mu2_from : float +- initial value of mu2 +- mu2_from : float +- final value of mu2 +- aem_list : list +- list of electromagnetic coupling values +- alphaem_running : bool +- whether alphaem is running or not +- nf : int +- number of active flavors +- L : float +- logarithm of the squared ratio of factorization and renormalization scale +- ev_op_iterations : int +- number of evolution steps +- ev_op_max_order : int +- perturbative expansion order of U +- sv_mode: int, `enum.IntEnum` +- scale variation mode, see `eko.scale_variations.Modes` +- is_threshold : boolean +- is this an intermediate threshold operator? +- n3lo_ad_variation : tuple +- |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` +- is_polarized : boolean +- is polarized evolution ? +- is_time_like : boolean +- is time-like evolution ? +- +- Returns +- ------- +- float +- evaluated integration kernel +- """ +- ker_base = QuadKerBase(u, is_log, logx, mode0) +- integrand = ker_base.integrand(areas) +- if integrand == 0.0: +- return 0.0 +- if order[1] == 0: +- ker = quad_ker_qcd( +- ker_base, +- order, +- mode0, +- mode1, +- method, +- as_list[-1], +- as_list[0], +- nf, +- L, +- ev_op_iterations, +- ev_op_max_order, +- sv_mode, +- is_threshold, +- is_polarized, +- is_time_like, +- n3lo_ad_variation, +- ) +- else: +- ker = quad_ker_qed( +- ker_base, +- order, +- mode0, +- mode1, +- method, +- as_list, +- mu2_from, +- mu2_to, +- a_half, +- alphaem_running, +- nf, +- L, +- ev_op_iterations, +- ev_op_max_order, +- sv_mode, +- is_threshold, +- n3lo_ad_variation, +- ) +- +- # recombine everything +- return np.real(ker * integrand) +- +- +-@nb.njit(cache=True) +-def quad_ker_qcd( +- ker_base, +- order, +- mode0, +- mode1, +- method, +- as1, +- as0, +- nf, +- L, +- ev_op_iterations, +- ev_op_max_order, +- sv_mode, +- is_threshold, +- is_polarized, +- is_time_like, +- n3lo_ad_variation, +-): +- """Raw evolution kernel inside quad. +- +- Parameters +- ---------- +- quad_ker : float +- quad argument +- order : int +- perturbation order +- mode0: int +- pid for first sector element +- mode1 : int +- pid for second sector element +- method : str +- method +- as1 : float +- target coupling value +- as0 : float +- initial coupling value +- nf : int +- number of active flavors +- L : float +- logarithm of the squared ratio of factorization and renormalization scale +- ev_op_iterations : int +- number of evolution steps +- ev_op_max_order : int +- perturbative expansion order of U +- sv_mode: int, `enum.IntEnum` +- scale variation mode, see `eko.scale_variations.Modes` +- is_threshold : boolean +- is this an itermediate threshold operator? +- n3lo_ad_variation : tuple +- |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` +- +- Returns +- ------- +- float +- evaluated integration kernel +- """ +- # compute the actual evolution kernel for pure QCD +- if ker_base.is_singlet: +- if is_polarized: +- if is_time_like: +- raise NotImplementedError("Polarized, time-like is not implemented") +- else: +- gamma_singlet = ad_ps.gamma_singlet(order, ker_base.n, nf) +- else: +- if is_time_like: +- gamma_singlet = ad_ut.gamma_singlet(order, ker_base.n, nf) +- else: +- gamma_singlet = ad_us.gamma_singlet( +- order, ker_base.n, nf, n3lo_ad_variation +- ) +- # scale var exponentiated is directly applied on gamma +- if sv_mode == sv.Modes.exponentiated: +- gamma_singlet = sv.exponentiated.gamma_variation( +- gamma_singlet, order, nf, L +- ) +- ker = s.dispatcher( +- order, +- method, +- gamma_singlet, +- as1, +- as0, +- nf, +- ev_op_iterations, +- ev_op_max_order, +- ) +- # scale var expanded is applied on the kernel +- if sv_mode == sv.Modes.expanded and not is_threshold: +- ker = np.ascontiguousarray( +- sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L, dim=2) +- ) @ np.ascontiguousarray(ker) +- ker = select_singlet_element(ker, mode0, mode1) +- else: +- if is_polarized: +- if is_time_like: +- raise NotImplementedError("Polarized, time-like is not implemented") +- else: +- gamma_ns = ad_ps.gamma_ns(order, mode0, ker_base.n, nf) +- else: +- if is_time_like: +- gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf) +- else: +- gamma_ns = ad_us.gamma_ns(order, mode0, ker_base.n, nf) +- if sv_mode == sv.Modes.exponentiated: +- gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) +- ker = ns.dispatcher( +- order, +- method, +- gamma_ns, +- as1, +- as0, +- nf, +- ev_op_iterations, +- ) +- if sv_mode == sv.Modes.expanded and not is_threshold: +- ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker +- return ker +- +- +-@nb.njit(cache=True) +-def quad_ker_qed( +- ker_base, +- order, +- mode0, +- mode1, +- method, +- as_list, +- mu2_from, +- mu2_to, +- a_half, +- alphaem_running, +- nf, +- L, +- ev_op_iterations, +- ev_op_max_order, +- sv_mode, +- is_threshold, +- n3lo_ad_variation, +-): +- """Raw evolution kernel inside quad. +- +- Parameters +- ---------- +- ker_base : QuadKerBase +- quad argument +- order : int +- perturbation order +- mode0: int +- pid for first sector element +- mode1 : int +- pid for second sector element +- method : str +- method +- as1 : float +- target coupling value +- as0 : float +- initial coupling value +- mu2_from : float +- initial value of mu2 +- mu2_from : float +- final value of mu2 +- aem_list : list +- list of electromagnetic coupling values +- alphaem_running : bool +- whether alphaem is running or not +- nf : int +- number of active flavors +- L : float +- logarithm of the squared ratio of factorization and renormalization scale +- ev_op_iterations : int +- number of evolution steps +- ev_op_max_order : int +- perturbative expansion order of U +- sv_mode: int, `enum.IntEnum` +- scale variation mode, see `eko.scale_variations.Modes` +- is_threshold : boolean +- is this an itermediate threshold operator? +- n3lo_ad_variation : tuple +- |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` +- +- Returns +- ------- +- float +- evaluated integration kernel +- """ +- # compute the actual evolution kernel for QEDxQCD +- if ker_base.is_QEDsinglet: +- gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation) +- # scale var exponentiated is directly applied on gamma +- if sv_mode == sv.Modes.exponentiated: +- gamma_s = sv.exponentiated.gamma_variation_qed( +- gamma_s, order, nf, L, alphaem_running +- ) +- ker = qed_s.dispatcher( +- order, +- method, +- gamma_s, +- as_list, +- a_half, +- nf, +- ev_op_iterations, +- ev_op_max_order, +- ) +- # scale var expanded is applied on the kernel +- # TODO : in this way a_half[-1][1] is the aem value computed in +- # the middle point of the last step. Instead we want aem computed in mu2_final. +- # However the distance between the two is very small and affects only the running aem +- if sv_mode == sv.Modes.expanded and not is_threshold: +- ker = np.ascontiguousarray( +- sv.expanded.singlet_variation_qed( +- gamma_s, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +- ) +- ) @ np.ascontiguousarray(ker) +- ker = select_QEDsinglet_element(ker, mode0, mode1) +- elif ker_base.is_QEDvalence: +- gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) +- # scale var exponentiated is directly applied on gamma +- if sv_mode == sv.Modes.exponentiated: +- gamma_v = sv.exponentiated.gamma_variation_qed( +- gamma_v, order, nf, L, alphaem_running +- ) +- ker = qed_v.dispatcher( +- order, +- method, +- gamma_v, +- as_list, +- a_half, +- nf, +- ev_op_iterations, +- ev_op_max_order, +- ) +- # scale var expanded is applied on the kernel +- if sv_mode == sv.Modes.expanded and not is_threshold: +- ker = np.ascontiguousarray( +- sv.expanded.valence_variation_qed( +- gamma_v, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +- ) +- ) @ np.ascontiguousarray(ker) +- ker = select_QEDvalence_element(ker, mode0, mode1) +- else: +- gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) +- # scale var exponentiated is directly applied on gamma +- if sv_mode == sv.Modes.exponentiated: +- gamma_ns = sv.exponentiated.gamma_variation_qed( +- gamma_ns, order, nf, L, alphaem_running +- ) +- ker = qed_ns.dispatcher( +- order, +- method, +- gamma_ns, +- as_list, +- a_half[:, 1], +- alphaem_running, +- nf, +- ev_op_iterations, +- mu2_from, +- mu2_to, +- ) +- if sv_mode == sv.Modes.expanded and not is_threshold: +- ker = ( +- sv.expanded.non_singlet_variation_qed( +- gamma_ns, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +- ) +- * ker +- ) +- return ker +- +- + class Operator(sv.ModeMixin): + """Internal representation of a single EKO. + +@@ -784,49 +305,6 @@ class Operator(sv.ModeMixin): + labels.extend(br.singlet_unified_labels) + return labels + +- def quad_ker(self, label, logx, areas): +- """Return partially initialized integrand function. +- +- Parameters +- ---------- +- label: tuple +- operator element pids +- logx: float +- Mellin inversion point +- areas : tuple +- basis function configuration +- +- Returns +- ------- +- functools.partial +- partially initialized integration kernel +- +- """ +- return functools.partial( +- quad_ker, +- order=self.order, +- mode0=label[0], +- mode1=label[1], +- method=self.config["method"], +- is_log=self.int_disp.log, +- logx=logx, +- areas=areas, +- as_list=self.as_list, +- mu2_from=self.q2_from, +- mu2_to=self.q2_to, +- a_half=self.a_half_list, +- alphaem_running=self.alphaem_running, +- nf=self.nf, +- L=np.log(self.xif2), +- ev_op_iterations=self.config["ev_op_iterations"], +- ev_op_max_order=tuple(self.config["ev_op_max_order"]), +- sv_mode=self.sv_mode, +- is_threshold=self.is_threshold, +- n3lo_ad_variation=self.config["n3lo_ad_variation"], +- is_polarized=self.config["polarized"], +- is_time_like=self.config["time_like"], +- ) +- + def initialize_op_members(self): + """Init all operators with the identity or zeros.""" + eye = OpMember( +@@ -849,10 +327,7 @@ class Operator(sv.ModeMixin): + else: + self.op_members[n] = zero.copy() + +- def run_op_integration( +- self, +- log_grid, +- ): ++ def run_op_integration(self, log_grid): + """Run the integration for each grid point. + + Parameters +@@ -867,18 +342,56 @@ class Operator(sv.ModeMixin): + """ + column = [] + k, logx = log_grid ++ # call(!) self.labels only once ++ labels = self.labels + start_time = time.perf_counter() ++ # start preparing C arguments ++ cfg = ekors.lib.empty_qcd_args() ++ cfg.order_qcd = self.order[0] ++ cfg.is_polarized = self.config["polarized"] ++ cfg.is_time_like = self.config["time_like"] ++ cfg.nf = self.nf ++ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qcd.address) ++ cfg.is_log = self.int_disp.log ++ cfg.logx = logx ++ cfg.L = np.log(self.xif2) ++ cfg.method_num = 1 ++ cfg.as1 = self.as_list[1] ++ cfg.as0 = self.as_list[0] ++ cfg.ev_op_iterations = self.config["ev_op_iterations"] ++ cfg.ev_op_max_order_qcd = self.config["ev_op_max_order"][0] ++ cfg.sv_mode_num = 1 ++ cfg.is_threshold = self.is_threshold ++ + # iterate basis functions + for l, bf in enumerate(self.int_disp): + if k == l and l == self.grid_size - 1: + continue ++ # add emtpy labels with 0s ++ if bf.is_below_x(np.exp(logx)): ++ column.append({label: (0.0, 0.0) for label in labels}) ++ continue + temp_dict = {} ++ # prepare areas for C ++ curareas = bf.areas_representation ++ areas_len = curareas.shape[0] * curareas.shape[1] ++ # force the variable in scope ++ areas_ffi = ekors.ffi.new( ++ f"double[{areas_len}]", curareas.flatten().tolist() ++ ) ++ cfg.areas = areas_ffi ++ cfg.areas_x = curareas.shape[0] ++ cfg.areas_y = curareas.shape[1] + # iterate sectors +- for label in self.labels: ++ for label in labels: ++ cfg.mode0 = label[0] ++ cfg.mode1 = label[1] ++ # construct the low level object ++ func = LowLevelCallable( ++ ekors.lib.rust_quad_ker_qcd, ekors.ffi.addressof(cfg) ++ ) + res = integrate.quad( +- self.quad_ker( +- label=label, logx=logx, areas=bf.areas_representation +- ), ++ func, + 0.5, + 1.0 - self._mellin_cut, + epsabs=1e-12, diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py new file mode 100644 index 000000000..9303f8aa9 --- /dev/null +++ b/src/eko/evolution_operator/quad_ker.py @@ -0,0 +1,367 @@ +r"""Integration kernels.""" + +import logging + +import numba as nb +import numpy as np +from scipy import LowLevelCallable, integrate + +from .. import basis_rotation as br +from .. import interpolation +from .. import scale_variations as sv +from ..kernels import non_singlet as ns +from ..kernels import singlet as s +from ..kernels import utils +from ..matchings import Segment +from ..member import OpMember + +logger = logging.getLogger(__name__) + + +@nb.njit(cache=True) +def select_singlet_element(ker, mode0, mode1): + """Select element of the singlet matrix. + + Parameters + ---------- + ker : numpy.ndarray + singlet integration kernel + mode0 : int + id for first sector element + mode1 : int + id for second sector element + + Returns + ------- + complex + singlet integration kernel element + """ + k = 0 if mode0 == 100 else 1 + l = 0 if mode1 == 100 else 1 + return ker[k, l] + + +@nb.cfunc( + nb.types.double( + nb.types.CPointer(nb.types.double), # re_gamma_ns_raw + nb.types.CPointer(nb.types.double), # im_gamma_ns_raw + nb.types.double, # re_n + nb.types.double, # im_n + nb.types.double, # re_jac + nb.types.double, # im_jac + nb.types.uint, # order_qcd + nb.types.bool_, # is_singlet + nb.types.uint, # mode0 + nb.types.uint, # mode1 + nb.types.uint, # nf + nb.types.bool_, # is_log + nb.types.double, # logx + nb.types.CPointer(nb.types.double), # areas_raw + nb.types.uint, # areas_x + nb.types.uint, # areas_y + nb.types.double, # L + nb.types.uint, # method_num + nb.types.double, # as1 + nb.types.double, # as0 + nb.types.uint, # ev_op_iterations + nb.types.uint, # ev_op_max_order_qcd + nb.types.uint, # sv_mode_num + nb.types.bool_, # is_threshold + ), + cache=True, +) +def cb_quad_ker_qcd( + re_gamma_raw, + im_gamma_raw, + re_n, + im_n, + re_jac, + im_jac, + order_qcd, + is_singlet, + mode0, + mode1, + nf, + is_log, + logx, + areas_raw, + areas_x, + areas_y, + L, + _method_num, + as1, + as0, + ev_op_iterations, + ev_op_max_order_qcd, + _sv_mode_num, + is_threshold, +): + """C Callback inside integration kernel.""" + # recover complex variables + n = re_n + im_n * 1j + jac = re_jac + im_jac * 1j + # combute basis functions + areas = nb.carray(areas_raw, (areas_x, areas_y)) + pj = interpolation.evaluate_grid(n, is_log, logx, areas) + # TODO recover parameters + method = "iterate-expanded" + sv_mode = sv.Modes.exponentiated + order = (order_qcd, 0) + ev_op_max_order = (ev_op_max_order_qcd, 0) + if is_singlet: + # reconstruct singlet matrices + re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, 2, 2)) + im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, 2, 2)) + gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation( + gamma_singlet, order, nf, L + ) + # construct eko + ker = s.dispatcher( + order, + method, + gamma_singlet, + as1, + as0, + nf, + ev_op_iterations, + ev_op_max_order, + ) + # scale var expanded is applied on the kernel + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = np.ascontiguousarray( + sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L, dim=2) + ) @ np.ascontiguousarray(ker) + ker = select_singlet_element(ker, mode0, mode1) + else: + # construct non-singlet matrices + re_gamma_ns = nb.carray(re_gamma_raw, order_qcd) + im_gamma_ns = nb.carray(im_gamma_raw, order_qcd) + gamma_ns = re_gamma_ns + im_gamma_ns * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) + # construct eko + ker = ns.dispatcher( + order, + method, + gamma_ns, + as1, + as0, + nf, + ev_op_iterations, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker + # recombine everything + res = ker * pj * jac + return np.real(res) + + +# from ..kernels import singlet_qed as qed_s +# from ..kernels import non_singlet_qed as qed_ns +# from ..kernels import valence_qed as qed_v + +# @nb.njit(cache=True) +# def select_QEDsinglet_element(ker, mode0, mode1): +# """Select element of the QEDsinglet matrix. + +# Parameters +# ---------- +# ker : numpy.ndarray +# QEDsinglet integration kernel +# mode0 : int +# id for first sector element +# mode1 : int +# id for second sector element +# Returns +# ------- +# ker : complex +# QEDsinglet integration kernel element +# """ +# if mode0 == 21: +# index1 = 0 +# elif mode0 == 22: +# index1 = 1 +# elif mode0 == 100: +# index1 = 2 +# else: +# index1 = 3 +# if mode1 == 21: +# index2 = 0 +# elif mode1 == 22: +# index2 = 1 +# elif mode1 == 100: +# index2 = 2 +# else: +# index2 = 3 +# return ker[index1, index2] + + +# @nb.njit(cache=True) +# def select_QEDvalence_element(ker, mode0, mode1): +# """ +# Select element of the QEDvalence matrix. + +# Parameters +# ---------- +# ker : numpy.ndarray +# QEDvalence integration kernel +# mode0 : int +# id for first sector element +# mode1 : int +# id for second sector element +# Returns +# ------- +# ker : complex +# QEDvalence integration kernel element +# """ +# index1 = 0 if mode0 == 10200 else 1 +# index2 = 0 if mode1 == 10200 else 1 +# return ker[index1, index2] + + +# @nb.njit(cache=True) +# def quad_ker_qed( +# ker_base, +# order, +# mode0, +# mode1, +# method, +# as_list, +# mu2_from, +# mu2_to, +# a_half, +# alphaem_running, +# nf, +# L, +# ev_op_iterations, +# ev_op_max_order, +# sv_mode, +# is_threshold, +# ): +# """Raw evolution kernel inside quad. + +# Parameters +# ---------- +# ker_base : QuadKerBase +# quad argument +# order : int +# perturbation order +# mode0: int +# pid for first sector element +# mode1 : int +# pid for second sector element +# method : str +# method +# as1 : float +# target coupling value +# as0 : float +# initial coupling value +# mu2_from : float +# initial value of mu2 +# mu2_from : float +# final value of mu2 +# aem_list : list +# list of electromagnetic coupling values +# alphaem_running : bool +# whether alphaem is running or not +# nf : int +# number of active flavors +# L : float +# logarithm of the squared ratio of factorization and renormalization scale +# ev_op_iterations : int +# number of evolution steps +# ev_op_max_order : int +# perturbative expansion order of U +# sv_mode: int, `enum.IntEnum` +# scale variation mode, see `eko.scale_variations.Modes` +# is_threshold : boolean +# is this an itermediate threshold operator? + +# Returns +# ------- +# float +# evaluated integration kernel +# """ +# # compute the actual evolution kernel for QEDxQCD +# if ker_base.is_QEDsinglet: +# gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf) +# # scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_s = sv.exponentiated.gamma_variation_qed( +# gamma_s, order, nf, L, alphaem_running +# ) +# ker = qed_s.dispatcher( +# order, +# method, +# gamma_s, +# as_list, +# a_half, +# nf, +# ev_op_iterations, +# ev_op_max_order, +# ) +# # scale var expanded is applied on the kernel +# # TODO : in this way a_half[-1][1] is the aem value computed in +# # the middle point of the last step. Instead we want aem computed in mu2_final. +# # However the distance between the two is very small and affects only the running aem +# if sv_mode == sv.Modes.expanded and not is_threshold: +# ker = np.ascontiguousarray( +# sv.expanded.singlet_variation_qed( +# gamma_s, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +# ) +# ) @ np.ascontiguousarray(ker) +# ker = select_QEDsinglet_element(ker, mode0, mode1) +# elif ker_base.is_QEDvalence: +# gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) +# # scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_v = sv.exponentiated.gamma_variation_qed( +# gamma_v, order, nf, L, alphaem_running +# ) +# ker = qed_v.dispatcher( +# order, +# method, +# gamma_v, +# as_list, +# a_half, +# nf, +# ev_op_iterations, +# ev_op_max_order, +# ) +# # scale var expanded is applied on the kernel +# if sv_mode == sv.Modes.expanded and not is_threshold: +# ker = np.ascontiguousarray( +# sv.expanded.valence_variation_qed( +# gamma_v, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +# ) +# ) @ np.ascontiguousarray(ker) +# ker = select_QEDvalence_element(ker, mode0, mode1) +# else: +# gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) +# # scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_ns = sv.exponentiated.gamma_variation_qed( +# gamma_ns, order, nf, L, alphaem_running +# ) +# ker = qed_ns.dispatcher( +# order, +# method, +# gamma_ns, +# as_list, +# a_half[:, 1], +# alphaem_running, +# nf, +# ev_op_iterations, +# mu2_from, +# mu2_to, +# ) +# if sv_mode == sv.Modes.expanded and not is_threshold: +# ker = ( +# sv.expanded.non_singlet_variation_qed( +# gamma_ns, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +# ) +# * ker +# ) +# return ker diff --git a/tests/eko/evolution_operator/test_init.py.patch b/tests/eko/evolution_operator/test_init.py.patch new file mode 100644 index 000000000..494f59a25 --- /dev/null +++ b/tests/eko/evolution_operator/test_init.py.patch @@ -0,0 +1,601 @@ +diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py +index 47d5e700..2315a2cb 100644 +--- a/tests/eko/evolution_operator/test_init.py ++++ b/tests/eko/evolution_operator/test_init.py +@@ -5,7 +5,6 @@ import pytest + import scipy.integrate + + import eko.runner.legacy +-import ekore.anomalous_dimensions.unpolarized.space_like as ad + from eko import basis_rotation as br + from eko import interpolation, mellin + from eko.evolution_operator import Operator, quad_ker +@@ -16,177 +15,176 @@ from eko.kernels import non_singlet_qed as qed_ns + from eko.kernels import singlet as s + from eko.matchings import Segment + +- +-def test_quad_ker_errors(): +- for mode0 in [br.non_singlet_pids_map["ns+"], 21]: +- with pytest.raises(NotImplementedError): +- quad_ker( +- u=0.3, +- order=(1, 0), +- mode0=mode0, +- mode1=0, +- method="", +- is_log=True, +- logx=0.1, +- areas=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], +- as_list=[2.0, 1.0], +- mu2_from=1.0, +- mu2_to=2.0, +- a_half=np.array([[1.5, 0.01]]), +- alphaem_running=False, +- nf=3, +- L=0, +- ev_op_iterations=1, +- ev_op_max_order=(1, 0), +- sv_mode=1, +- is_threshold=False, +- is_polarized=True, +- is_time_like=True, +- n3lo_ad_variation=(0, 0, 0, 0), +- ) ++# def test_quad_ker_errors(): ++# for mode0 in [br.non_singlet_pids_map["ns+"], 21]: ++# with pytest.raises(NotImplementedError): ++# quad_ker( ++# u=0.3, ++# order=(1, 0), ++# mode0=mode0, ++# mode1=0, ++# method="", ++# is_log=True, ++# logx=0.1, ++# areas=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], ++# as_list=[2.0, 1.0], ++# mu2_from=1.0, ++# mu2_to=2.0, ++# a_half=np.array([[1.5, 0.01]]), ++# alphaem_running=False, ++# nf=3, ++# L=0, ++# ev_op_iterations=1, ++# ev_op_max_order=(1, 0), ++# sv_mode=1, ++# is_threshold=False, ++# is_polarized=True, ++# is_time_like=True, ++# n3lo_ad_variation=(0, 0, 0, 0), ++# ) + + +-def test_quad_ker(monkeypatch): +- monkeypatch.setattr( +- mellin, "Talbot_path", lambda *args: 2 +- ) # N=2 is a safe evaluation point +- monkeypatch.setattr( +- mellin, "Talbot_jac", lambda *args: complex(0, np.pi) +- ) # negate mellin prefactor +- monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 1) +- monkeypatch.setattr(interpolation, "evaluate_Nx", lambda *args: 1) +- monkeypatch.setattr(ns, "dispatcher", lambda *args: 1.0) +- monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0) +- monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2)) +- params = [ +- ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0), +- ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.123, 1.0), +- ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0), +- ((1, 0), 100, 100, "", 0.123, 1.0), +- ((1, 0), 100, 21, "", 0.0, 0.0), +- ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0), +- ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0), +- ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0), +- ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0), +- ] +- for order, mode0, mode1, method, logx, res in params: +- for is_log in [True, False]: +- for t, p in [(False, False), (False, True), (True, False)]: +- res_ns = quad_ker( +- u=0, +- order=order, +- mode0=mode0, +- mode1=mode1, +- method=method, +- is_log=is_log, +- logx=logx, +- areas=np.zeros(3), +- as_list=[2.0, 1.0], +- mu2_from=1, +- mu2_to=2, +- a_half=np.array([[1.5, 0.01]]), +- alphaem_running=False, +- nf=3, +- L=0, +- ev_op_iterations=0, +- ev_op_max_order=(0, 0), +- sv_mode=1, +- is_threshold=False, +- is_polarized=p, +- is_time_like=t, +- n3lo_ad_variation=(0, 0, 0, 0), +- ) +- np.testing.assert_allclose(res_ns, res) +- for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]: +- for sv in [2, 3]: +- for polarized in [True, False]: +- res_sv = quad_ker( +- u=0, +- order=(1, 0), +- mode0=label[0], +- mode1=label[1], +- method="", +- is_log=True, +- logx=0.123, +- areas=np.zeros(3), +- as_list=[2.0, 1.0], +- mu2_from=1, +- mu2_to=2, +- a_half=np.array([[1.5, 0.01]]), +- alphaem_running=False, +- nf=3, +- L=0, +- ev_op_iterations=0, +- ev_op_max_order=(1, 0), +- sv_mode=sv, +- is_threshold=False, +- is_polarized=polarized, +- is_time_like=False, +- n3lo_ad_variation=(0, 0, 0, 0), +- ) +- np.testing.assert_allclose(res_sv, 1.0) +- for label in [ +- (100, 100), +- (21, 21), +- (22, 22), +- (101, 101), +- (10200, 10200), +- (10204, 10204), +- (10202, 0), +- ]: +- for sv in [2, 3]: +- res_sv = quad_ker( +- u=0, +- order=(1, 1), +- mode0=label[0], +- mode1=label[1], +- method="iterate-exact", +- is_log=True, +- logx=0.123, +- areas=np.zeros(3), +- as_list=[2.0, 1.0], +- mu2_from=1, +- mu2_to=2, +- a_half=np.array([[1.5, 0.01]]), +- alphaem_running=False, +- nf=3, +- L=0, +- ev_op_iterations=0, +- ev_op_max_order=(1, 0), +- sv_mode=sv, +- is_threshold=False, +- n3lo_ad_variation=(0, 0, 0, 0), +- is_polarized=False, +- is_time_like=False, +- ) +- np.testing.assert_allclose(res_sv, 1.0) ++# def test_quad_ker(monkeypatch): ++# monkeypatch.setattr( ++# mellin, "Talbot_path", lambda *args: 2 ++# ) # N=2 is a safe evaluation point ++# monkeypatch.setattr( ++# mellin, "Talbot_jac", lambda *args: complex(0, np.pi) ++# ) # negate mellin prefactor ++# monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 1) ++# monkeypatch.setattr(interpolation, "evaluate_Nx", lambda *args: 1) ++# monkeypatch.setattr(ns, "dispatcher", lambda *args: 1.0) ++# monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0) ++# monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2)) ++# params = [ ++# ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0), ++# ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.123, 1.0), ++# ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0), ++# ((1, 0), 100, 100, "", 0.123, 1.0), ++# ((1, 0), 100, 21, "", 0.0, 0.0), ++# ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0), ++# ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0), ++# ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0), ++# ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0), ++# ] ++# for order, mode0, mode1, method, logx, res in params: ++# for is_log in [True, False]: ++# for t, p in [(False, False), (False, True), (True, False)]: ++# res_ns = quad_ker( ++# u=0, ++# order=order, ++# mode0=mode0, ++# mode1=mode1, ++# method=method, ++# is_log=is_log, ++# logx=logx, ++# areas=np.zeros(3), ++# as_list=[2.0, 1.0], ++# mu2_from=1, ++# mu2_to=2, ++# a_half=np.array([[1.5, 0.01]]), ++# alphaem_running=False, ++# nf=3, ++# L=0, ++# ev_op_iterations=0, ++# ev_op_max_order=(0, 0), ++# sv_mode=1, ++# is_threshold=False, ++# is_polarized=p, ++# is_time_like=t, ++# n3lo_ad_variation=(0, 0, 0, 0), ++# ) ++# np.testing.assert_allclose(res_ns, res) ++# for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]: ++# for sv in [2, 3]: ++# for polarized in [True, False]: ++# res_sv = quad_ker( ++# u=0, ++# order=(1, 0), ++# mode0=label[0], ++# mode1=label[1], ++# method="", ++# is_log=True, ++# logx=0.123, ++# areas=np.zeros(3), ++# as_list=[2.0, 1.0], ++# mu2_from=1, ++# mu2_to=2, ++# a_half=np.array([[1.5, 0.01]]), ++# alphaem_running=False, ++# nf=3, ++# L=0, ++# ev_op_iterations=0, ++# ev_op_max_order=(1, 0), ++# sv_mode=sv, ++# is_threshold=False, ++# is_polarized=polarized, ++# is_time_like=False, ++# n3lo_ad_variation=(0, 0, 0, 0), ++# ) ++# np.testing.assert_allclose(res_sv, 1.0) ++# for label in [ ++# (100, 100), ++# (21, 21), ++# (22, 22), ++# (101, 101), ++# (10200, 10200), ++# (10204, 10204), ++# (10202, 0), ++# ]: ++# for sv in [2, 3]: ++# res_sv = quad_ker( ++# u=0, ++# order=(1, 1), ++# mode0=label[0], ++# mode1=label[1], ++# method="iterate-exact", ++# is_log=True, ++# logx=0.123, ++# areas=np.zeros(3), ++# as_list=[2.0, 1.0], ++# mu2_from=1, ++# mu2_to=2, ++# a_half=np.array([[1.5, 0.01]]), ++# alphaem_running=False, ++# nf=3, ++# L=0, ++# ev_op_iterations=0, ++# ev_op_max_order=(1, 0), ++# sv_mode=sv, ++# is_threshold=False, ++# n3lo_ad_variation=(0, 0, 0, 0), ++# is_polarized=False, ++# is_time_like=False, ++# ) ++# np.testing.assert_allclose(res_sv, 1.0) + +- monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 0) +- res_ns = quad_ker( +- u=0, +- order=(1, 0), +- mode0=br.non_singlet_pids_map["ns+"], +- mode1=0, +- method="", +- is_log=True, +- logx=0.0, +- areas=np.zeros(3), +- as_list=[2.0, 1.0], +- mu2_from=1, +- mu2_to=2, +- a_half=np.array([[1.5, 0.01]]), +- alphaem_running=False, +- nf=3, +- L=0, +- ev_op_iterations=0, +- ev_op_max_order=(0, 0), +- sv_mode=1, +- is_threshold=False, +- n3lo_ad_variation=(0, 0, 0, 0), +- is_polarized=False, +- is_time_like=False, +- ) +- np.testing.assert_allclose(res_ns, 0.0) ++# monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 0) ++# res_ns = quad_ker( ++# u=0, ++# order=(1, 0), ++# mode0=br.non_singlet_pids_map["ns+"], ++# mode1=0, ++# method="", ++# is_log=True, ++# logx=0.0, ++# areas=np.zeros(3), ++# as_list=[2.0, 1.0], ++# mu2_from=1, ++# mu2_to=2, ++# a_half=np.array([[1.5, 0.01]]), ++# alphaem_running=False, ++# nf=3, ++# L=0, ++# ev_op_iterations=0, ++# ev_op_max_order=(0, 0), ++# sv_mode=1, ++# is_threshold=False, ++# n3lo_ad_variation=(0, 0, 0, 0), ++# is_polarized=False, ++# is_time_like=False, ++# ) ++# np.testing.assert_allclose(res_ns, 0.0) + + + class FakeCoupling: +@@ -319,30 +317,30 @@ class TestOperator: + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, + ) + +- def test_compute_no_skip_sv( +- self, monkeypatch, theory_ffns, operator_card, tmp_path +- ): +- tcard: TheoryCard = theory_ffns(3) +- tcard.xif = 2.0 +- ocard: OperatorCard = operator_card +- ocard.configs.scvar_method = ScaleVariationsMethod.EXPANDED +- r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar") +- g = r.op_grid +- # setup objs +- o = Operator(g.config, g.managers, Segment(2.0, 2.0, 3)) +- # fake quad +- v = 0.1234 +- monkeypatch.setattr( +- scipy.integrate, "quad", lambda *args, v=v, **kwargs: (v, 0.56) +- ) +- o.compute() +- lx = len(ocard.xgrid.raw) +- res = np.full((lx, lx), v) +- res[-1, -1] = 1.0 +- # ns are all diagonal, so they start from an identity matrix +- for k in br.non_singlet_labels: +- assert k in o.op_members +- np.testing.assert_allclose(o.op_members[k].value, res, err_msg=k) ++ # def test_compute_no_skip_sv( ++ # self, monkeypatch, theory_ffns, operator_card, tmp_path ++ # ): ++ # tcard: TheoryCard = theory_ffns(3) ++ # tcard.xif = 2.0 ++ # ocard: OperatorCard = operator_card ++ # ocard.configs.scvar_method = ScaleVariationsMethod.EXPANDED ++ # r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar") ++ # g = r.op_grid ++ # # setup objs ++ # o = Operator(g.config, g.managers, Segment(2.0, 2.0, 3)) ++ # # fake quad ++ # v = 0.1234 ++ # monkeypatch.setattr( ++ # scipy.integrate, "quad", lambda *args, v=v, **kwargs: (v, 0.56) ++ # ) ++ # o.compute() ++ # lx = len(ocard.xgrid.raw) ++ # res = np.full((lx, lx), v) ++ # res[-1, -1] = 1.0 ++ # # ns are all diagonal, so they start from an identity matrix ++ # for k in br.non_singlet_labels: ++ # assert k in o.op_members ++ # np.testing.assert_allclose(o.op_members[k].value, res, err_msg=k) + + def test_compute(self, monkeypatch, theory_ffns, operator_card, tmp_path): + tcard: TheoryCard = theory_ffns(3) +@@ -395,97 +393,97 @@ class TestOperator: + ) + + +-def test_pegasus_path(): +- def quad_ker_pegasus( +- u, order, mode0, method, logx, areas, a1, a0, nf, ev_op_iterations +- ): +- # compute the mellin inversion as done in pegasus +- phi = 3 / 4 * np.pi +- c = 1.9 +- n = complex(c + u * np.exp(1j * phi)) +- gamma_ns = ad.gamma_ns(order, mode0, n, nf) +- ker = ns.dispatcher( +- order, +- method, +- gamma_ns, +- a1, +- a0, +- nf, +- ev_op_iterations, +- ) +- pj = interpolation.log_evaluate_Nx(n, logx, areas) +- return np.imag(np.exp(1j * phi) / np.pi * pj * ker) ++# def test_pegasus_path(): ++# def quad_ker_pegasus( ++# u, order, mode0, method, logx, areas, a1, a0, nf, ev_op_iterations ++# ): ++# # compute the mellin inversion as done in pegasus ++# phi = 3 / 4 * np.pi ++# c = 1.9 ++# n = complex(c + u * np.exp(1j * phi)) ++# gamma_ns = ad.gamma_ns(order, mode0, n, nf) ++# ker = ns.dispatcher( ++# order, ++# method, ++# gamma_ns, ++# a1, ++# a0, ++# nf, ++# ev_op_iterations, ++# ) ++# pj = interpolation.log_evaluate_Nx(n, logx, areas) ++# return np.imag(np.exp(1j * phi) / np.pi * pj * ker) + +- # It might be useful to test with a different function +- # monkeypatch.setattr(ns, "dispatcher", lambda x, *args: np.exp( - x ** 2 ) ) +- xgrid = np.geomspace(1e-7, 1, 10) +- int_disp = InterpolatorDispatcher(xgrid, 1, True) +- order = (2, 0) +- mode0 = br.non_singlet_pids_map["ns+"] +- mode1 = 0 +- method = "" +- logxs = np.log(int_disp.xgrid.raw) +- a1 = 1 +- a0 = 2 +- mu2_from = 1 +- mu2_to = 2**2 +- nf = 3 +- L = 0 +- ev_op_iterations = 10 +- for logx in logxs: +- for bf in int_disp: +- res_ns, _ = scipy.integrate.quad( +- quad_ker, +- 0.5, +- 1.0, +- args=( +- order, +- mode0, +- mode1, +- method, +- int_disp.log, +- logx, +- bf.areas_representation, +- [a0, a1], +- mu2_from, +- mu2_to, +- [[(a1 + a0) / 2, 0.00058]], +- False, +- nf, +- L, +- ev_op_iterations, +- 10, +- 0, +- False, +- (0, 0, 0, 0), +- False, +- False, +- ), +- epsabs=1e-12, +- epsrel=1e-5, +- limit=100, +- full_output=1, +- )[:2] ++# # It might be useful to test with a different function ++# # monkeypatch.setattr(ns, "dispatcher", lambda x, *args: np.exp( - x ** 2 ) ) ++# xgrid = np.geomspace(1e-7, 1, 10) ++# int_disp = InterpolatorDispatcher(xgrid, 1, True) ++# order = (2, 0) ++# mode0 = br.non_singlet_pids_map["ns+"] ++# mode1 = 0 ++# method = "" ++# logxs = np.log(int_disp.xgrid.raw) ++# a1 = 1 ++# a0 = 2 ++# mu2_from = 1 ++# mu2_to = 2**2 ++# nf = 3 ++# L = 0 ++# ev_op_iterations = 10 ++# for logx in logxs: ++# for bf in int_disp: ++# res_ns, _ = scipy.integrate.quad( ++# quad_ker, ++# 0.5, ++# 1.0, ++# args=( ++# order, ++# mode0, ++# mode1, ++# method, ++# int_disp.log, ++# logx, ++# bf.areas_representation, ++# [a0, a1], ++# mu2_from, ++# mu2_to, ++# [[(a1 + a0) / 2, 0.00058]], ++# False, ++# nf, ++# L, ++# ev_op_iterations, ++# 10, ++# 0, ++# False, ++# (0, 0, 0, 0), ++# False, ++# False, ++# ), ++# epsabs=1e-12, ++# epsrel=1e-5, ++# limit=100, ++# full_output=1, ++# )[:2] + +- res_test, _ = scipy.integrate.quad( +- quad_ker_pegasus, +- 0, +- np.inf, +- args=( +- order, +- mode0, +- method, +- logx, +- bf.areas_representation, +- a1, +- a0, +- nf, +- ev_op_iterations, +- ), +- epsabs=1e-12, +- epsrel=1e-5, +- limit=100, +- full_output=1, +- )[:2] ++# res_test, _ = scipy.integrate.quad( ++# quad_ker_pegasus, ++# 0, ++# np.inf, ++# args=( ++# order, ++# mode0, ++# method, ++# logx, ++# bf.areas_representation, ++# a1, ++# a0, ++# nf, ++# ev_op_iterations, ++# ), ++# epsabs=1e-12, ++# epsrel=1e-5, ++# limit=100, ++# full_output=1, ++# )[:2] + +- np.testing.assert_allclose(res_ns, res_test, rtol=2e-6) ++# np.testing.assert_allclose(res_ns, res_test, rtol=2e-6)