From 176284667f9a422da3b7767e57e314babe0d3a74 Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 16 Oct 2024 15:29:56 +0200 Subject: [PATCH 01/71] start implementation of aem1 --- .../unpolarized/spacelike.rs | 1 + .../unpolarized/spacelike/aem1.rs | 62 +++++++++++++++++++ crates/ekore/src/constants.rs | 29 +++++++++ 3 files changed, 92 insertions(+) create mode 100644 crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index dec73bf04..3ac896f93 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -4,6 +4,7 @@ use crate::constants::{PID_NSM, PID_NSP, PID_NSV}; use crate::harmonics::cache::Cache; use num::complex::Complex; use num::Zero; +pub mod aem1; pub mod as1; pub mod as2; pub mod as3; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs new file mode 100644 index 000000000..3a2b728fd --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -0,0 +1,62 @@ +//! |LO| |QED|. +use num::complex::Complex; + +use crate::constants::{ed2, eu2, uplike_flavors, CF, NC, TR}; +use crate::harmonics::cache::Cache; + +use crate::anomalous_dimensions::unpolarized::spacelike::as1; + +/// Compute the leading-order photon-quark anomalous dimension. +/// +/// Implements Eq. (2.5) of +pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex { + as1::gamma_gq(c, nf) / CF +} + +/// Compute the leading-order quark-photon anomalous dimension. +/// +/// Implements Eq. (2.5) of +pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { + as1::gamma_qg(c, nf) / TR * (NC as f64) +} + +/// Compute the leading-order photon-photon anomalous dimension. +/// +/// Implements Eq. (2.5) of +pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { + let nu = uplike_flavors(nf); + let nd = nf - nu; + (4.0 / 3.0 * (NC as f64) * ((nu as f64) * eu2 + (nd as f64) * ed2)).into() +} + +/// Compute the leading-order non-singlet QED anomalous dimension +/// +/// Implements Eq. (2.5) of +pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { + as1::gamma_ns(c, nf) / CF +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + 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_cmplx!(f64, me, Complex::zero(), 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_phq(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); + } +} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 8cb90eb1f..377cb1f77 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -1,4 +1,5 @@ //! Global constants. +use std::unimplemented; /// The number of colors. /// @@ -20,6 +21,16 @@ pub const CA: f64 = NC as f64; /// 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); +/// Up quark charge square. +/// +/// Defaults to $e_u^2 = 4./9$ +pub const eu2: f64 = 4. / 9.; + +/// Down quark charge square. +/// +/// Defaults to $e_d^2 = 1./9$ +pub const ed2: f64 = 1. / 9.; + /// Riemann zeta function at z = 2. /// /// $\zeta(2) = \pi^2 / 6$. @@ -41,3 +52,21 @@ pub const PID_NSM: u16 = 10201; /// non-singlet all-valence |PID|. pub const PID_NSV: u16 = 10200; + +/// compute the number of up flavors +pub fn uplike_flavors(nf: u8) -> u8 { + if nf > 6 { + unimplemented!("Selected nf is not implemented") + } + nf / 2 +} + +pub fn charge_combinations(nf: u8) -> Vec { + let nu = uplike_flavors(nf) as f64; + let nd = (nf as f64) - nu; + let e2avg = (nu * eu2 + nd * ed2) / (nf as f64); + let vue2m = nu / (nf as f64) * (eu2 - ed2); + let vde2m = nd / (nf as f64) * (eu2 - ed2); + let e2delta = vde2m - vue2m + e2avg; + vec![e2avg, vue2m, vde2m, e2delta] +} From 61220f8a3340a7274d81fdc314b729fbf9634f80 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 21 Oct 2024 11:26:08 +0200 Subject: [PATCH 02/71] complete aem1.py --- .../unpolarized/spacelike/aem1.rs | 61 ++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 3a2b728fd..91427622c 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -1,7 +1,8 @@ //! |LO| |QED|. use num::complex::Complex; +use num::Zero; -use crate::constants::{ed2, eu2, uplike_flavors, CF, NC, TR}; +use crate::constants::{charge_combinations, ed2, eu2, uplike_flavors, CF, NC, TR}; use crate::harmonics::cache::Cache; use crate::anomalous_dimensions::unpolarized::spacelike::as1; @@ -36,6 +37,64 @@ pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { as1::gamma_ns(c, nf) / CF } +/// Compute the leading-order singlet QED anomalous dimension matrix +/// +/// Implements Eq. (2.5) of +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + let cc = charge_combinations(nf); + let e2avg = cc[0]; + let vue2m = cc[1]; + let vde2m = cc[2]; + let e2delta = cc[3]; + + let gamma_ph_q = gamma_phq(c, nf); + let gamma_q_ph = gamma_qph(c, nf); + let gamma_nonsinglet = gamma_ns(c, nf); + + [ + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + Complex::::zero(), + gamma_phph(c, nf), + e2avg * gamma_ph_q, + vue2m * gamma_ph_q, + ], + [ + Complex::::zero(), + e2avg * gamma_q_ph, + e2avg * gamma_nonsinglet, + vue2m * gamma_nonsinglet, + ], + [ + Complex::::zero(), + vde2m * gamma_q_ph, + vde2m * gamma_nonsinglet, + e2delta * gamma_nonsinglet, + ], + ] +} + +/// Compute the leading-order valence QED anomalous dimension matrix +/// +/// Implements Eq. (2.5) of +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let cc = charge_combinations(nf); + let e2avg = cc[0]; + let vue2m = cc[1]; + let vde2m = cc[2]; + let e2delta = cc[3]; + + [ + [e2avg * gamma_ns(c, nf), vue2m * gamma_ns(c, nf)], + [vde2m * gamma_ns(c, nf), e2delta * gamma_ns(c, nf)], + ] +} + #[cfg(test)] mod tests { use super::*; From a40908b8a11b91a6a938868b81f159ff41a8844a Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 21 Oct 2024 11:35:08 +0200 Subject: [PATCH 03/71] unit tests --- .../unpolarized/spacelike/aem1.rs | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 91427622c..68e30694c 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -118,4 +118,21 @@ mod tests { let me = gamma_ns(&mut c, NF) + gamma_phq(&mut c, NF); assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); } + + #[test] + fn photon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let nu = uplike_flavors(nf); + let nd = nf - nu; + assert_approx_eq_cmplx!( + f64, + eu2 * gamma_qph(&mut c, nu) + ed2 * gamma_qph(&mut c, nd) + gamma_phph(&mut c, nf), + cmplx!(0., 0.), + epsilon = 2e-6 + ); + } + } } From f0edfae5134c029689c98a40fa0131a5e7a2dd11 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 21 Oct 2024 15:21:20 +0200 Subject: [PATCH 04/71] add struct for charge combinations --- .../unpolarized/spacelike/aem1.rs | 28 +++++++------------ crates/ekore/src/constants.rs | 17 +++++++++-- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 68e30694c..8d8d97af1 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -42,10 +42,6 @@ pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { /// Implements Eq. (2.5) of pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = charge_combinations(nf); - let e2avg = cc[0]; - let vue2m = cc[1]; - let vde2m = cc[2]; - let e2delta = cc[3]; let gamma_ph_q = gamma_phq(c, nf); let gamma_q_ph = gamma_qph(c, nf); @@ -61,20 +57,20 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { [ Complex::::zero(), gamma_phph(c, nf), - e2avg * gamma_ph_q, - vue2m * gamma_ph_q, + cc.e2avg * gamma_ph_q, + cc.vue2m * gamma_ph_q, ], [ Complex::::zero(), - e2avg * gamma_q_ph, - e2avg * gamma_nonsinglet, - vue2m * gamma_nonsinglet, + cc.e2avg * gamma_q_ph, + cc.e2avg * gamma_nonsinglet, + cc.vue2m * gamma_nonsinglet, ], [ Complex::::zero(), - vde2m * gamma_q_ph, - vde2m * gamma_nonsinglet, - e2delta * gamma_nonsinglet, + cc.vde2m * gamma_q_ph, + cc.vde2m * gamma_nonsinglet, + cc.e2delta * gamma_nonsinglet, ], ] } @@ -84,14 +80,10 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { /// Implements Eq. (2.5) of pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = charge_combinations(nf); - let e2avg = cc[0]; - let vue2m = cc[1]; - let vde2m = cc[2]; - let e2delta = cc[3]; [ - [e2avg * gamma_ns(c, nf), vue2m * gamma_ns(c, nf)], - [vde2m * gamma_ns(c, nf), e2delta * gamma_ns(c, nf)], + [cc.e2avg * gamma_ns(c, nf), cc.vue2m * gamma_ns(c, nf)], + [cc.vde2m * gamma_ns(c, nf), cc.e2delta * gamma_ns(c, nf)], ] } diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 377cb1f77..725b28c8b 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -61,12 +61,25 @@ pub fn uplike_flavors(nf: u8) -> u8 { nf / 2 } -pub fn charge_combinations(nf: u8) -> Vec { +pub struct ChargeCombinations { + pub e2avg: f64, + pub vue2m: f64, + pub vde2m: f64, + pub e2delta: f64, +} + +pub fn charge_combinations(nf: u8) -> ChargeCombinations { let nu = uplike_flavors(nf) as f64; let nd = (nf as f64) - nu; let e2avg = (nu * eu2 + nd * ed2) / (nf as f64); let vue2m = nu / (nf as f64) * (eu2 - ed2); let vde2m = nd / (nf as f64) * (eu2 - ed2); let e2delta = vde2m - vue2m + e2avg; - vec![e2avg, vue2m, vde2m, e2delta] + + ChargeCombinations { + e2avg, + vue2m, + vde2m, + e2delta, + } } From a8d72c02ea5b723e02030c09e2266ebebe517b82 Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 23 Oct 2024 15:58:53 +0200 Subject: [PATCH 05/71] some work on gamma_ns_qed --- .../unpolarized/spacelike.rs | 27 ++++++++++++++++++- crates/ekore/src/constants.rs | 9 +++++++ 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 3ac896f93..7ff05f431 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -1,6 +1,8 @@ //! The unpolarized, space-like anomalous dimensions at various couplings power. -use crate::constants::{PID_NSM, PID_NSP, PID_NSV}; +use crate::constants::{ + ed2, eu2, PID_NSM, PID_NSM_ED2, PID_NSM_EU2, PID_NSP, PID_NSP_ED2, PID_NSP_EU2, PID_NSV, +}; use crate::harmonics::cache::Cache; use num::complex::Complex; use num::Zero; @@ -56,3 +58,26 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl } gamma_S } + +/// Compute the tower of the QED non-singlet anomalous dimensions. +pub fn gamma_ns_qed( + order_qcd: usize, + order_qed: usize, + mode: u16, + c: &mut Cache, + nf: u8, +) -> Vec>> { + let row = vec![Complex::::zero(); order_qcd + 1]; + let mut gamma_ns = vec![row; order_qed + 1]; + gamma_ns[1][0] = as1::gamma_ns(c, nf); + gamma_ns[0][1] = choose_ns_as_aem1(mode, c, nf); + gamma_ns +} + +pub fn choose_ns_as_aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { + match mode { + PID_NSP_EU2 | PID_NSM_EU2 => eu2 * aem1::gamma_ns(c, nf), + PID_NSP_ED2 | PID_NSM_ED2 => ed2 * aem1::gamma_ns(c, nf), + _ => panic!("Unkown non-singlet sector element"), + } +} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 725b28c8b..7dd00d088 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -53,6 +53,15 @@ pub const PID_NSM: u16 = 10201; /// non-singlet all-valence |PID|. pub const PID_NSV: u16 = 10200; +/// QED |PID|. Need to give sensible names +pub const PID_NSP_EU2: u16 = 10102; + +pub const PID_NSP_ED2: u16 = 10103; + +pub const PID_NSM_EU2: u16 = 10202; + +pub const PID_NSM_ED2: u16 = 10203; + /// compute the number of up flavors pub fn uplike_flavors(nf: u8) -> u8 { if nf > 6 { From e8a3a8bcf8fe7af9d8115c8c9e705d373d6c83af Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 23 Oct 2024 16:22:24 +0200 Subject: [PATCH 06/71] start as1aem1.rs --- .../unpolarized/spacelike.rs | 1 + .../unpolarized/spacelike/as1aem1.rs | 37 +++++++++++++++++++ 2 files changed, 38 insertions(+) create mode 100644 crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 7ff05f431..7a5fc13e8 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -8,6 +8,7 @@ use num::complex::Complex; use num::Zero; pub mod aem1; pub mod as1; +pub mod as1aem1; pub mod as2; pub mod as3; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs new file mode 100644 index 000000000..2c9bef29e --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -0,0 +1,37 @@ +//! The $O(a_s^1a_{em}^1)$ Altarelli-Parisi splitting kernels. +use num::complex::Complex; + +use crate::constants::CF; +use crate::harmonics::cache::{Cache, K}; + +/// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. +/// +/// Implements Eq. (36) of +pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + + #[rustfmt::skip] + let tmp_const = + 2.0 + * ( + -4.0 + - 12.0 * N + - N.powu(2) + + 28.0 * N.powu(3) + + 43.0 * N.powu(4) + + 30.0 * N.powu(5) + + 12.0 * N.powu(6) + ) / ((-1.0 + N) * N.powu(3) * (1.0 + N).powu(3)); + + #[rustfmt::skip] + let tmp_S1 = -4.0 + * (10.0 + 27.0 * N + 25.0 * N.powu(2) + 13.0 * N.powu(3) + 5.0 * N.powu(4)) + / ((-1.0 + N) * N * (1.0 + N).powu(3)); + + let tmp_S12 = 4.0 * (2.0 + N + N.powu(2)) / ((-1.0 + N) * N * (1.0 + N)); + let tmp_S2 = 4.0 * (2.0 + N + N.powu(2)) / ((-1.0 + N) * N * (1.0 + N)); + + CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) +} From 136dbe8d32cf80ab222f1a3e2d3fcc5eaff8d2ed Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 23 Oct 2024 16:33:38 +0200 Subject: [PATCH 07/71] gamma_qph as1aem1 --- .../unpolarized/spacelike/as1aem1.rs | 27 ++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 2c9bef29e..429d8094e 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -1,7 +1,7 @@ //! The $O(a_s^1a_{em}^1)$ Altarelli-Parisi splitting kernels. use num::complex::Complex; -use crate::constants::CF; +use crate::constants::{CA, CF}; use crate::harmonics::cache::{Cache, K}; /// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. @@ -35,3 +35,28 @@ pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex { CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) } + +/// Compute the $O(a_s^1a_{em}^1)$ quark-photon anomalous dimension. +/// +/// Implements Eq. (26) of +pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + + let tmp_const = -2.0 + * (4.0 + + 8.0 * N + + 25.0 * N.powu(2) + + 51.0 * N.powu(3) + + 36.0 * N.powu(4) + + 15.0 * N.powu(5) + + 5.0 * N.powu(6)) + / (N.powu(3) * (1.0 + N).powu(3) * (2.0 + N)); + + let tmp_S1 = 8.0 / N.powu(2); + let tmp_S12 = -4.0 * (2.0 + N + N.powu(2)) / (N * (1.0 + N) * (2.0 + N)); + let tmp_S2 = 4.0 * (2.0 + N + N.powu(2)) / (N * (1.0 + N) * (2.0 + N)); + + 2.0 * (nf as f64) * CA * CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) +} From ce11227a32ecba72bfddbdeef0c3e91863291b35 Mon Sep 17 00:00:00 2001 From: tgiani Date: Fri, 25 Oct 2024 12:07:35 +0200 Subject: [PATCH 08/71] some work on as1aem1 --- .../unpolarized/spacelike/as1aem1.rs | 50 ++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 429d8094e..6a77690aa 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -1,7 +1,8 @@ //! The $O(a_s^1a_{em}^1)$ Altarelli-Parisi splitting kernels. +use crate::cmplx; use num::complex::Complex; -use crate::constants::{CA, CF}; +use crate::constants::{ed2, eu2, uplike_flavors, CA, CF, NC, TR}; use crate::harmonics::cache::{Cache, K}; /// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. @@ -60,3 +61,50 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { 2.0 * (nf as f64) * CA * CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) } + +/// Compute the $O(a_s^1a_{em}^1)$ gluon-photon anomalous dimension. +/// +/// Implements Eq. (27) of +pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + CF * CA + * (8.0 * (-4.0 + N * (-4.0 + N * (-5.0 + N * (-10.0 + N + 2.0 * N.powu(2) * (2.0 + N)))))) + / (N.powu(3) * (1.0 + N).powu(3) * (-2.0 + N + N.powu(2))) +} + +/// Compute the $O(a_s^1a_{em}^1)$ photon-gluon anomalous dimension. +/// +/// Implements Eq. (30) of +pub fn gamma_phg(c: &mut Cache, nf: u8) -> Complex { + TR / CF / CA * (NC as f64) * gamma_gph(c, nf) +} + +/// Compute the $O(a_s^1a_{em}^1)$ quark-gluon singlet anomalous dimension. +/// +/// Implements Eq. (29) of +pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { + TR / CF / CA * (NC as f64) * gamma_qph(c, nf) +} + +/// Compute the $O(a_s^1a_{em}^1)$ gluon-quark singlet anomalous dimension. +/// +/// Implements Eq. (35) of +pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { + gamma_phq(c, nf) +} + +/// Compute the $O(a_s^1a_{em}^1)$ photon-photon singlet anomalous dimension. +/// +/// Implements Eq. (28) of +pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { + let nu = uplike_flavors(nf); + let nd = nf - nu; + cmplx!(4.0 * CF * CA * ((nu as f64) * eu2 + (nd as f64) * ed2), 0.) +} + +/// Compute the $O(a_s^1a_{em}^1)$ gluon-gluon singlet anomalous dimension. +/// +/// Implements Eq. (31) of +pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex { + cmplx!(4.0 * TR * (NC as f64), 0.) +} From 659fe3112f4e840569b927369905a16c4dcb4961 Mon Sep 17 00:00:00 2001 From: tgiani Date: Fri, 25 Oct 2024 14:01:41 +0200 Subject: [PATCH 09/71] addingv recursive harmonics to cache --- crates/ekore/src/harmonics/cache.rs | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs index da0f4800b..f8dfd487c 100644 --- a/crates/ekore/src/harmonics/cache.rs +++ b/crates/ekore/src/harmonics/cache.rs @@ -46,6 +46,10 @@ pub enum K { Sm21e, /// $S_{-2,1}(N)$ odd moments Sm21o, + /// recursive harmonics + S1ph, + S2ph, + S3ph, } /// Hold all cached values. @@ -98,6 +102,9 @@ impl Cache { K::Sm3o => w3::Sm3o(self.get(K::S3), self.get(K::S3mh)), K::Sm21e => w3::Sm21e(self.n, self.get(K::S1), self.get(K::Sm1e)), K::Sm21o => w3::Sm21o(self.n, self.get(K::S1), self.get(K::Sm1o)), + K::S1ph => recursive_harmonic_sum(self.get(K::S1mh), (self.n - 1.) / 2., 1, 1), + K::S2ph => recursive_harmonic_sum(self.get(K::S2mh), (self.n - 1.) / 2., 1, 2), + K::S3ph => recursive_harmonic_sum(self.get(K::S3mh), (self.n - 1.) / 2., 1, 3), }; // insert self.m.insert(k, val); From 4007f0dedd5e09869d0cfd449209d6d2d0097be8 Mon Sep 17 00:00:00 2001 From: tgiani Date: Fri, 25 Oct 2024 14:24:11 +0200 Subject: [PATCH 10/71] S1p2 and G3p2. Check that this is correct --- crates/ekore/src/harmonics/cache.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs index f8dfd487c..a9dc6ad25 100644 --- a/crates/ekore/src/harmonics/cache.rs +++ b/crates/ekore/src/harmonics/cache.rs @@ -50,6 +50,8 @@ pub enum K { S1ph, S2ph, S3ph, + S1p2, + G3p2, } /// Hold all cached values. @@ -94,6 +96,7 @@ impl Cache { K::S2mh => w2::S2((self.n - 1.) / 2.), K::S3mh => w3::S3((self.n - 1.) / 2.), K::G3 => g_functions::g3(self.n, self.get(K::S1)), + K::G3p2 => g_functions::g3(self.n + 2., self.get(K::S1p2)), K::Sm1e => w1::Sm1e(self.get(K::S1), self.get(K::S1h)), K::Sm1o => w1::Sm1o(self.get(K::S1), self.get(K::S1mh)), K::Sm2e => w2::Sm2e(self.get(K::S2), self.get(K::S2h)), @@ -105,6 +108,7 @@ impl Cache { K::S1ph => recursive_harmonic_sum(self.get(K::S1mh), (self.n - 1.) / 2., 1, 1), K::S2ph => recursive_harmonic_sum(self.get(K::S2mh), (self.n - 1.) / 2., 1, 2), K::S3ph => recursive_harmonic_sum(self.get(K::S3mh), (self.n - 1.) / 2., 1, 3), + K::S1p2 => recursive_harmonic_sum(self.get(K::S1), self.n, 2, 1), }; // insert self.m.insert(k, val); From 7c160f1d9da7018a331d718b2bb68ac6a83bfa41 Mon Sep 17 00:00:00 2001 From: tgiani Date: Fri, 25 Oct 2024 14:29:58 +0200 Subject: [PATCH 11/71] gamma_nsp --- .../unpolarized/spacelike/as1aem1.rs | 34 ++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 6a77690aa..2eedf14ec 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -2,7 +2,7 @@ use crate::cmplx; use num::complex::Complex; -use crate::constants::{ed2, eu2, uplike_flavors, CA, CF, NC, TR}; +use crate::constants::{ed2, eu2, uplike_flavors, CA, CF, NC, TR, ZETA2, ZETA3}; use crate::harmonics::cache::{Cache, K}; /// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. @@ -45,6 +45,7 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { let S1 = c.get(K::S1); let S2 = c.get(K::S2); + #[rustfmt::skip] let tmp_const = -2.0 * (4.0 + 8.0 * N @@ -108,3 +109,34 @@ pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex { cmplx!(4.0 * TR * (NC as f64), 0.) } + +/// Compute the $O(a_s^1a_{em}^1)$ singlet-like non singlet anomalous dimension. +/// +/// Implements Eqs. (33-34) of +pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S3 = c.get(K::S3); + let S1h = c.get(K::S1h); + let S2h = c.get(K::S2h); + let S3h = c.get(K::S3h); + let S1p1h = c.get(K::S1ph); + let S2p1h = c.get(K::S2ph); + let S3p1h = c.get(K::S3ph); + let g3N = c.get(K::G3); + let g3Np2 = c.get(K::G3p2); + + #[rustfmt::skip] + let result = 32.0 * ZETA2 * S1h - 32.0 * ZETA2 * S1p1h + + 8.0 / (N + N.powu(2)) * S2h + - 4.0 * S3h + (24.0 + 16.0 / (N + N.powu(2))) * S2 + - 32.0 * S3 - 8.0 / (N + N.powu(2)) * S2p1h + + S1 * (16.0 * (3.0 / N.powu(2) - 3.0 / (1.0 + N).powu(2) + 2.0 * ZETA2) - 16.0 * S2h + - 32.0 * S2 + 16.0 * S2p1h ) + + (-8.0 + N * (-32.0 + N * ( -8.0 - 3.0 * N * (3.0 + N) * (3.0 + N.powu(2)) - 48.0 * (1.0 + N).powu(2) * ZETA2))) + / (N.powu(3) * (1.0 + N).powu(3)) + + 32.0 * (g3N + g3Np2) + 4.0 * S3p1h - 16.0 * ZETA3; + + CF * result +} From 37a65e8b066502729e2af48621da83caa8c6f5c3 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 11 Nov 2024 15:53:22 +0100 Subject: [PATCH 12/71] rewrite ChargeCombination struct --- .../unpolarized/spacelike/aem1.rs | 26 +++++----- .../unpolarized/spacelike/as1aem1.rs | 51 +++++++++++++++++++ crates/ekore/src/constants.rs | 41 ++++++++------- 3 files changed, 88 insertions(+), 30 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 8d8d97af1..bc70ed97f 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -2,7 +2,7 @@ use num::complex::Complex; use num::Zero; -use crate::constants::{charge_combinations, ed2, eu2, uplike_flavors, CF, NC, TR}; +use crate::constants::{ed2, eu2, uplike_flavors, ChargeCombinations, CF, NC, TR}; use crate::harmonics::cache::Cache; use crate::anomalous_dimensions::unpolarized::spacelike::as1; @@ -41,7 +41,7 @@ pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (2.5) of pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { - let cc = charge_combinations(nf); + let cc = ChargeCombinations { nf }; let gamma_ph_q = gamma_phq(c, nf); let gamma_q_ph = gamma_qph(c, nf); @@ -57,20 +57,20 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { [ Complex::::zero(), gamma_phph(c, nf), - cc.e2avg * gamma_ph_q, - cc.vue2m * gamma_ph_q, + cc.e2avg() * gamma_ph_q, + cc.vue2m() * gamma_ph_q, ], [ Complex::::zero(), - cc.e2avg * gamma_q_ph, - cc.e2avg * gamma_nonsinglet, - cc.vue2m * gamma_nonsinglet, + cc.e2avg() * gamma_q_ph, + cc.e2avg() * gamma_nonsinglet, + cc.vue2m() * gamma_nonsinglet, ], [ Complex::::zero(), - cc.vde2m * gamma_q_ph, - cc.vde2m * gamma_nonsinglet, - cc.e2delta * gamma_nonsinglet, + cc.vde2m() * gamma_q_ph, + cc.vde2m() * gamma_nonsinglet, + cc.e2delta() * gamma_nonsinglet, ], ] } @@ -79,11 +79,11 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { /// /// Implements Eq. (2.5) of pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { - let cc = charge_combinations(nf); + let cc = ChargeCombinations { nf }; [ - [cc.e2avg * gamma_ns(c, nf), cc.vue2m * gamma_ns(c, nf)], - [cc.vde2m * gamma_ns(c, nf), cc.e2delta * gamma_ns(c, nf)], + [cc.e2avg() * gamma_ns(c, nf), cc.vue2m() * gamma_ns(c, nf)], + [cc.vde2m() * gamma_ns(c, nf), cc.e2delta() * gamma_ns(c, nf)], ] } diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 2eedf14ec..917ea45a7 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -124,6 +124,7 @@ pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { let S1p1h = c.get(K::S1ph); let S2p1h = c.get(K::S2ph); let S3p1h = c.get(K::S3ph); + let g3N = c.get(K::G3); let g3Np2 = c.get(K::G3p2); @@ -140,3 +141,53 @@ pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { CF * result } + +/// Compute the $O(a_s^1a_{em}^1)$ valence-like non singlet anomalous dimension. +/// +/// Implements Eqs. (33-34) of +pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S3 = c.get(K::S3); + let S1h = c.get(K::S1h); + let S2h = c.get(K::S2h); + let S3h = c.get(K::S3h); + let S1p1h = c.get(K::S1ph); + let S2p1h = c.get(K::S2ph); + let S3p1h = c.get(K::S3ph); + let g3N = c.get(K::G3); + let g3Np2 = c.get(K::G3p2); + + #[rustfmt::skip] + let result = + -32.0 * ZETA2 * S1h + - 8.0 / (N + N.powu(2)) * S2h + + (24.0 + 16.0 / (N + N.powu(2))) * S2 + + 8.0 / (N + N.powu(2)) * S2p1h + + S1 + * ( + 16.0 * (-1.0 / N.powu(2) + 1.0 / (1.0 + N).powu(2) + 2.0 * ZETA2) + + 16.0 * S2h + - 32.0 * S2 + - 16.0 * S2p1h + ) + + ( + 72.0 + + N + * ( + 96.0 + - 3.0 * N * (8.0 + 3.0 * N * (3.0 + N) * (3.0 + N.powu(2))) + + 48.0 * N * (1.0 + N).powu(2) * ZETA2 + ) + ) + / (3.0 * N.powu(3) * (1.0 + N).powu(3)) + - 32.0 * (g3N + g3Np2) + + 32.0 * ZETA2 * S1p1h + + 4.0 * S3h + - 32.0 * S3 + - 4.0 * S3p1h + - 16.0 * ZETA3; + + CF * result +} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 7dd00d088..80c3bc807 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -71,24 +71,31 @@ pub fn uplike_flavors(nf: u8) -> u8 { } pub struct ChargeCombinations { - pub e2avg: f64, - pub vue2m: f64, - pub vde2m: f64, - pub e2delta: f64, + pub nf: u8, } -pub fn charge_combinations(nf: u8) -> ChargeCombinations { - let nu = uplike_flavors(nf) as f64; - let nd = (nf as f64) - nu; - let e2avg = (nu * eu2 + nd * ed2) / (nf as f64); - let vue2m = nu / (nf as f64) * (eu2 - ed2); - let vde2m = nd / (nf as f64) * (eu2 - ed2); - let e2delta = vde2m - vue2m + e2avg; - - ChargeCombinations { - e2avg, - vue2m, - vde2m, - e2delta, +impl ChargeCombinations { + pub fn nu(&self) -> u8 { + uplike_flavors(self.nf) + } + + pub fn nd(&self) -> u8 { + self.nf - self.nu() + } + + pub fn e2avg(&self) -> f64 { + (self.nu() as f64 * eu2 + self.nd() as f64 * ed2) / (self.nf as f64) + } + + pub fn vue2m(&self) -> f64 { + self.nu() as f64 / (self.nf as f64) * (eu2 - ed2) + } + + pub fn vde2m(&self) -> f64 { + self.nd() as f64 / (self.nf as f64) * (eu2 - ed2) + } + + pub fn e2delta(&self) -> f64 { + self.vde2m() - self.vue2m() + self.e2avg() } } From fc5de81d706fd282db40508f6e5324418c240c39 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 11 Nov 2024 16:14:21 +0100 Subject: [PATCH 13/71] missing entries in as1aem1 --- .../unpolarized/spacelike/as1aem1.rs | 53 ++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 917ea45a7..19f598466 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -2,7 +2,9 @@ use crate::cmplx; use num::complex::Complex; -use crate::constants::{ed2, eu2, uplike_flavors, CA, CF, NC, TR, ZETA2, ZETA3}; +use crate::constants::{ + ed2, eu2, uplike_flavors, ChargeCombinations, CA, CF, NC, TR, ZETA2, ZETA3, +}; use crate::harmonics::cache::{Cache, K}; /// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. @@ -191,3 +193,52 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { CF * result } + +/// Compute the $O(a_s^1a_{em}^1)$ singlet sector. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + let cc = ChargeCombinations { nf }; + // let e2avg = cc.e2avg(); + // let vue2m = cc.vue2m(); + // let vde2m = cc.vde2m(); + // let e2delta = cc.e2delta(); + let e2_tot = nf as f64 * cc.e2avg(); + + [ + [ + e2_tot * gamma_gg(c, nf), + e2_tot * gamma_gph(c, nf), + cc.e2avg() * gamma_gq(c, nf), + cc.vue2m() * gamma_gq(c, nf), + ], + [ + e2_tot * gamma_phg(c, nf), + gamma_phph(c, nf), + cc.e2avg() * gamma_phq(c, nf), + cc.vue2m() * gamma_phq(c, nf), + ], + [ + cc.e2avg() * gamma_qg(c, nf), + cc.e2avg() * gamma_qph(c, nf), + cc.e2avg() * gamma_nsp(c, nf), + cc.vue2m() * gamma_nsp(c, nf), + ], + [ + cc.vde2m() * gamma_qg(c, nf), + cc.vde2m() * gamma_qph(c, nf), + cc.vde2m() * gamma_nsp(c, nf), + cc.e2delta() * gamma_nsp(c, nf), + ], + ] +} + +/// Compute the $O(a_s^1a_{em}^1)$ valence sector. +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let cc = ChargeCombinations { nf }; + [ + [cc.e2avg() * gamma_nsm(c, nf), cc.vue2m() * gamma_nsm(c, nf)], + [ + cc.vde2m() * gamma_nsm(c, nf) * gamma_nsm(c, nf), + cc.e2delta() * gamma_nsm(c, nf), + ], + ] +} From 2a911b1b22ca45b33f698700649d4f0a38073db8 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 11 Nov 2024 16:38:29 +0100 Subject: [PATCH 14/71] adding tests --- .../unpolarized/spacelike/as1aem1.rs | 65 +++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 19f598466..b96ed3680 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -242,3 +242,68 @@ pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { ], ] } + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + const NF: u8 = 5; + + #[test] + fn number_conservation() { + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + let me = gamma_nsm(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-4); + } + + #[test] + fn gluon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let nu = uplike_flavors(nf); + let nd = nf - nu; + assert_approx_eq_cmplx!( + f64, + eu2 * gamma_qg(&mut c, nu) + + ed2 * gamma_qg(&mut c, nd) + + (nu as f64 * eu2 + nd as f64 * ed2) * gamma_phg(&mut c, nf) + + (nu as f64 * eu2 + nd as f64 * ed2) * gamma_gg(&mut c, nf), + cmplx!(0., 0.), + epsilon = 1e-14 + ); + } + } + + #[test] + fn photon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let nu = uplike_flavors(nf); + let nd = nf - nu; + assert_approx_eq_cmplx!( + f64, + eu2 * gamma_qph(&mut c, nu) + + ed2 * gamma_qph(&mut c, nd) + + gamma_phph(&mut c, nf) + + (nu as f64 * eu2 + nd as f64 * ed2) * gamma_gph(&mut c, nf), + cmplx!(0., 0.), + epsilon = 1e-14 + ); + } + } + + #[test] + fn quark_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + let me = gamma_nsp(&mut c, NF) + gamma_gq(&mut c, NF) + gamma_phq(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-4); + } +} From bad124ac40aca464729d962e05ef4bb4e6553adb Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 20 Nov 2024 11:10:25 +0100 Subject: [PATCH 15/71] some work on spacelike.rs and recasting from list to vec --- .../unpolarized/spacelike.rs | 39 +++++++++++++++++++ .../unpolarized/spacelike/aem1.rs | 4 +- .../unpolarized/spacelike/as1.rs | 32 +++++++++++++++ .../unpolarized/spacelike/as1aem1.rs | 4 +- 4 files changed, 75 insertions(+), 4 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 7a5fc13e8..f41e52266 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -72,6 +72,7 @@ pub fn gamma_ns_qed( let mut gamma_ns = vec![row; order_qed + 1]; gamma_ns[1][0] = as1::gamma_ns(c, nf); gamma_ns[0][1] = choose_ns_as_aem1(mode, c, nf); + gamma_ns[1][1] = choose_ns_as_as1aem1(mode, c, nf); gamma_ns } @@ -82,3 +83,41 @@ pub fn choose_ns_as_aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { _ => panic!("Unkown non-singlet sector element"), } } + +pub fn choose_ns_as_as1aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { + match mode { + PID_NSP_EU2 => eu2 * as1aem1::gamma_nsp(c, nf), + PID_NSP_ED2 => ed2 * as1aem1::gamma_nsp(c, nf), + PID_NSM_EU2 => eu2 * as1aem1::gamma_nsm(c, nf), + PID_NSM_ED2 => ed2 * as1aem1::gamma_nsm(c, nf), + _ => panic!("Unkown non-singlet sector element"), + } +} + +pub fn gamma_singlet_qed( + order_qcd: usize, + order_qed: usize, + mode: u16, + c: &mut Cache, + nf: u8, +) -> Vec; 4]>>> { + let mut row = vec![ + vec![ + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero() + ]; + 4 + ]; + order_qcd + 1 + ]; + + let mut gamma_s = vec![row; order_qed + 1]; + + gamma_s[1][0] = as1::gamma_singlet_qed(c, nf); + gamma_s[0][1] = aem1::gamma_singlet(c, nf); + gamma_s[1][1] = as1aem1::gamma_singlet(c, nf); + gamma_s +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index bc70ed97f..7791d9552 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -40,14 +40,14 @@ pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { /// Compute the leading-order singlet QED anomalous dimension matrix /// /// Implements Eq. (2.5) of -pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { let cc = ChargeCombinations { nf }; let gamma_ph_q = gamma_phq(c, nf); let gamma_q_ph = gamma_qph(c, nf); let gamma_nonsinglet = gamma_ns(c, nf); - [ + vec![ [ Complex::::zero(), Complex::::zero(), diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs index 6d25e1c51..dad16e819 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -1,6 +1,7 @@ //! |LO| |QCD|. use num::complex::Complex; +use num::Zero; use crate::constants::{CA, CF, TR}; use crate::harmonics::cache::{Cache, K}; @@ -52,6 +53,37 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { ] } +/// Compute the leading-order singlet anomalous dimension matrix +/// for the unified evolution basis. +pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { + vec![ + [ + gamma_gg(c, nf), + Complex::::zero(), + gamma_gq(c, nf), + Complex::::zero(), + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + gamma_qg(c, nf), + Complex::::zero(), + gamma_ns(c, nf), + Complex::::zero(), + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + gamma_ns(c, nf), + ], + ] +} + #[cfg(test)] mod tests { use super::*; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index b96ed3680..a9ae80770 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -195,7 +195,7 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { } /// Compute the $O(a_s^1a_{em}^1)$ singlet sector. -pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { let cc = ChargeCombinations { nf }; // let e2avg = cc.e2avg(); // let vue2m = cc.vue2m(); @@ -203,7 +203,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { // let e2delta = cc.e2delta(); let e2_tot = nf as f64 * cc.e2avg(); - [ + vec![ [ e2_tot * gamma_gg(c, nf), e2_tot * gamma_gph(c, nf), From d50e9c735f248189486e2990e1f8726df3a482ee Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 20 Nov 2024 15:13:56 +0100 Subject: [PATCH 16/71] valence qed --- .../unpolarized/spacelike.rs | 20 +++++++++++++++++-- .../unpolarized/spacelike/aem1.rs | 4 ++-- .../unpolarized/spacelike/as1.rs | 9 +++++++++ .../unpolarized/spacelike/as1aem1.rs | 8 ++------ 4 files changed, 31 insertions(+), 10 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index f41e52266..e7f94273f 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -94,14 +94,14 @@ pub fn choose_ns_as_as1aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { } } +/// Compute the grid of the QED singlet anomalous dimensions matrices pub fn gamma_singlet_qed( order_qcd: usize, order_qed: usize, - mode: u16, c: &mut Cache, nf: u8, ) -> Vec; 4]>>> { - let mut row = vec![ + let row = vec![ vec![ [ Complex::::zero(), @@ -121,3 +121,19 @@ pub fn gamma_singlet_qed( gamma_s[1][1] = as1aem1::gamma_singlet(c, nf); gamma_s } + +/// Compute the grid of the QED valence anomalous dimensions matrices +pub fn gamma_valence_qed( + order_qcd: usize, + order_qed: usize, + c: &mut Cache, + nf: u8, +) -> Vec; 2]>>> { + let row = vec![vec![[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; + + let mut gamma_v = vec![row; order_qed + 1]; + gamma_v[1][0] = as1::gamma_valence_qed(c, nf); + gamma_v[0][1] = aem1::gamma_valence(c, nf); + gamma_v[1][1] = as1aem1::gamma_valence(c, nf); + gamma_v +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 7791d9552..ded4825d7 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -78,10 +78,10 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { /// Compute the leading-order valence QED anomalous dimension matrix /// /// Implements Eq. (2.5) of -pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { +pub fn gamma_valence(c: &mut Cache, nf: u8) -> Vec<[Complex; 2]> { let cc = ChargeCombinations { nf }; - [ + vec![ [cc.e2avg() * gamma_ns(c, nf), cc.vue2m() * gamma_ns(c, nf)], [cc.vde2m() * gamma_ns(c, nf), cc.e2delta() * gamma_ns(c, nf)], ] diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs index dad16e819..8075804de 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -84,6 +84,15 @@ pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { ] } +/// Compute the leading-order valence anomalous dimension matrix +/// for the unified evolution basis. +pub fn gamma_valence_qed(c: &mut Cache, nf: u8) -> Vec<[Complex; 2]> { + vec![ + [gamma_ns(c, nf), Complex::::zero()], + [Complex::::zero(), gamma_ns(c, nf)], + ] +} + #[cfg(test)] mod tests { use super::*; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index a9ae80770..a7ee20db9 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -197,10 +197,6 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { /// Compute the $O(a_s^1a_{em}^1)$ singlet sector. pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { let cc = ChargeCombinations { nf }; - // let e2avg = cc.e2avg(); - // let vue2m = cc.vue2m(); - // let vde2m = cc.vde2m(); - // let e2delta = cc.e2delta(); let e2_tot = nf as f64 * cc.e2avg(); vec![ @@ -232,9 +228,9 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { } /// Compute the $O(a_s^1a_{em}^1)$ valence sector. -pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { +pub fn gamma_valence(c: &mut Cache, nf: u8) -> Vec<[Complex; 2]> { let cc = ChargeCombinations { nf }; - [ + vec![ [cc.e2avg() * gamma_nsm(c, nf), cc.vue2m() * gamma_nsm(c, nf)], [ cc.vde2m() * gamma_nsm(c, nf) * gamma_nsm(c, nf), From b0f6722f7281f2c682140ff81969a177c3996529 Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 20 Nov 2024 16:21:02 +0100 Subject: [PATCH 17/71] some unit tests --- .../unpolarized/spacelike.rs | 91 +++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index e7f94273f..5d50ca6b5 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -137,3 +137,94 @@ pub fn gamma_valence_qed( gamma_v[1][1] = as1aem1::gamma_valence(c, nf); gamma_v } + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + + #[test] + fn gamma_ns() { + const NF: u8 = 3; + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + assert_approx_eq_cmplx!( + f64, + gamma_ns_qcd(3, PID_NSP, &mut c, NF)[0], + cmplx!(0., 0.), + epsilon = 1e-14 + ); + + for i in [0, 1] { + assert_approx_eq_cmplx!( + f64, + gamma_ns_qcd(2, PID_NSM, &mut c, NF)[i], + cmplx!(0., 0.), + epsilon = 2e-6 + ); + } + + for i in 0..3 { + assert_approx_eq_cmplx!( + f64, + gamma_ns_qcd(3, PID_NSM, &mut c, NF)[i], + cmplx!(0., 0.), + epsilon = 2e-4 + ); + } + + for i in 0..3 { + assert_approx_eq_cmplx!( + f64, + gamma_ns_qcd(3, PID_NSV, &mut c, NF)[i], + cmplx!(0., 0.), + epsilon = 8e-4 + ); + } + } + + #[test] + fn test_gamma_ns_qed() { + const NF: u8 = 3; + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + + for i in [0, 1] { + for j in [0, 1] { + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, PID_NSM_EU2, &mut c, NF)[i][j], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + } + } + + for i in [0, 1] { + for j in [0, 1] { + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, PID_NSM_ED2, &mut c, NF)[i][j], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + } + } + + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, PID_NSP_EU2, &mut c, NF)[0][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, PID_NSP_ED2, &mut c, NF)[0][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + } +} From e4c621b4dbdeea576bdc77c510d76f342aa318ba Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 20 Nov 2024 17:25:44 +0100 Subject: [PATCH 18/71] use Vec<> only when dimension is not known, else use normal list --- .../unpolarized/spacelike.rs | 22 ++++++++----------- .../unpolarized/spacelike/aem1.rs | 8 +++---- .../unpolarized/spacelike/as1.rs | 8 +++---- .../unpolarized/spacelike/as1aem1.rs | 8 +++---- 4 files changed, 21 insertions(+), 25 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 5d50ca6b5..810e004a9 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -94,23 +94,19 @@ pub fn choose_ns_as_as1aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { } } -/// Compute the grid of the QED singlet anomalous dimensions matrices pub fn gamma_singlet_qed( order_qcd: usize, order_qed: usize, c: &mut Cache, nf: u8, -) -> Vec; 4]>>> { +) -> Vec; 4]; 4]>> { let row = vec![ - vec![ - [ - Complex::::zero(), - Complex::::zero(), - Complex::::zero(), - Complex::::zero() - ]; - 4 - ]; + [[ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero() + ]; 4]; order_qcd + 1 ]; @@ -128,8 +124,8 @@ pub fn gamma_valence_qed( order_qed: usize, c: &mut Cache, nf: u8, -) -> Vec; 2]>>> { - let row = vec![vec![[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; +) -> Vec; 2]; 2]>> { + let row = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; let mut gamma_v = vec![row; order_qed + 1]; gamma_v[1][0] = as1::gamma_valence_qed(c, nf); diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index ded4825d7..bc70ed97f 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -40,14 +40,14 @@ pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { /// Compute the leading-order singlet QED anomalous dimension matrix /// /// Implements Eq. (2.5) of -pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = ChargeCombinations { nf }; let gamma_ph_q = gamma_phq(c, nf); let gamma_q_ph = gamma_qph(c, nf); let gamma_nonsinglet = gamma_ns(c, nf); - vec![ + [ [ Complex::::zero(), Complex::::zero(), @@ -78,10 +78,10 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { /// Compute the leading-order valence QED anomalous dimension matrix /// /// Implements Eq. (2.5) of -pub fn gamma_valence(c: &mut Cache, nf: u8) -> Vec<[Complex; 2]> { +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; - vec![ + [ [cc.e2avg() * gamma_ns(c, nf), cc.vue2m() * gamma_ns(c, nf)], [cc.vde2m() * gamma_ns(c, nf), cc.e2delta() * gamma_ns(c, nf)], ] diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs index 8075804de..bd2c67045 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -55,8 +55,8 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { /// Compute the leading-order singlet anomalous dimension matrix /// for the unified evolution basis. -pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { - vec![ +pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + [ [ gamma_gg(c, nf), Complex::::zero(), @@ -86,8 +86,8 @@ pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { /// Compute the leading-order valence anomalous dimension matrix /// for the unified evolution basis. -pub fn gamma_valence_qed(c: &mut Cache, nf: u8) -> Vec<[Complex; 2]> { - vec![ +pub fn gamma_valence_qed(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + [ [gamma_ns(c, nf), Complex::::zero()], [Complex::::zero(), gamma_ns(c, nf)], ] diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index a7ee20db9..546acca85 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -195,11 +195,11 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { } /// Compute the $O(a_s^1a_{em}^1)$ singlet sector. -pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = ChargeCombinations { nf }; let e2_tot = nf as f64 * cc.e2avg(); - vec![ + [ [ e2_tot * gamma_gg(c, nf), e2_tot * gamma_gph(c, nf), @@ -228,9 +228,9 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> Vec<[Complex; 4]> { } /// Compute the $O(a_s^1a_{em}^1)$ valence sector. -pub fn gamma_valence(c: &mut Cache, nf: u8) -> Vec<[Complex; 2]> { +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; - vec![ + [ [cc.e2avg() * gamma_nsm(c, nf), cc.vue2m() * gamma_nsm(c, nf)], [ cc.vde2m() * gamma_nsm(c, nf) * gamma_nsm(c, nf), From 81ece22721a2d051646706da31b055e0fddd906b Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 20 Nov 2024 17:36:28 +0100 Subject: [PATCH 19/71] fix notation row/columns --- .../anomalous_dimensions/unpolarized/spacelike.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 810e004a9..38479b89f 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -68,8 +68,8 @@ pub fn gamma_ns_qed( c: &mut Cache, nf: u8, ) -> Vec>> { - let row = vec![Complex::::zero(); order_qcd + 1]; - let mut gamma_ns = vec![row; order_qed + 1]; + let col = vec![Complex::::zero(); order_qcd + 1]; + let mut gamma_ns = vec![col; order_qed + 1]; gamma_ns[1][0] = as1::gamma_ns(c, nf); gamma_ns[0][1] = choose_ns_as_aem1(mode, c, nf); gamma_ns[1][1] = choose_ns_as_as1aem1(mode, c, nf); @@ -100,7 +100,7 @@ pub fn gamma_singlet_qed( c: &mut Cache, nf: u8, ) -> Vec; 4]; 4]>> { - let row = vec![ + let col = vec![ [[ Complex::::zero(), Complex::::zero(), @@ -110,7 +110,7 @@ pub fn gamma_singlet_qed( order_qcd + 1 ]; - let mut gamma_s = vec![row; order_qed + 1]; + let mut gamma_s = vec![col; order_qed + 1]; gamma_s[1][0] = as1::gamma_singlet_qed(c, nf); gamma_s[0][1] = aem1::gamma_singlet(c, nf); @@ -125,9 +125,9 @@ pub fn gamma_valence_qed( c: &mut Cache, nf: u8, ) -> Vec; 2]; 2]>> { - let row = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; + let col = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; - let mut gamma_v = vec![row; order_qed + 1]; + let mut gamma_v = vec![col; order_qed + 1]; gamma_v[1][0] = as1::gamma_valence_qed(c, nf); gamma_v[0][1] = aem1::gamma_valence(c, nf); gamma_v[1][1] = as1aem1::gamma_valence(c, nf); From 0def93893a5c8631b287ddbf50e4a7f8ae3a5dd8 Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 20 Nov 2024 17:43:55 +0100 Subject: [PATCH 20/71] add unravel functions for qed case. To be checked --- crates/eko/src/lib.rs | 45 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index a6c5fa9ae..f8c7cb1d8 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -14,6 +14,7 @@ struct RawCmplx { } /// Map tensors to c-ordered list +/// (res is a vector with dim order_qcd filled with DIMxDIM matrices) fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: usize) -> RawCmplx { let mut target = RawCmplx { re: Vec::::new(), @@ -30,6 +31,50 @@ fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: us target } +/// Map tensors to c-ordered list in the QED singlet and valence case +/// (res is a matrix with dim order_qcd x order_qed filled with DIMxDIM matrices) +fn unravel_qed_singlet( + res: Vec; DIM]; DIM]>>, + order_qcd: usize, + order_qed: usize, +) -> RawCmplx { + let mut target = RawCmplx { + re: Vec::::new(), + im: Vec::::new(), + }; + for obj_ in res.iter().take(order_qcd) { + for obj in obj_.iter().take(order_qed) { + for col in obj.iter().take(DIM) { + for el in col.iter().take(DIM) { + target.re.push(el.re); + target.im.push(el.im); + } + } + } + } + target +} + +/// Map tensors to c-ordered list in the QED singlet and valence case +/// (res is a matrix with dim order_qcd x order_qed filled with complex numbers) +fn unravel_qed_ns( + res: Vec>>, + order_qcd: usize, + order_qed: usize, +) -> RawCmplx { + let mut target = RawCmplx { + re: Vec::::new(), + im: Vec::::new(), + }; + for col in res.iter().take(order_qcd) { + for el in col.iter().take(order_qed) { + target.re.push(el.re); + target.im.push(el.im); + } + } + target +} + /// QCD intergration kernel inside quad. /// /// # Safety From 9b58455df37c253ad15c8928a02cb5bc4d2d7d3e Mon Sep 17 00:00:00 2001 From: tgiani Date: Thu, 21 Nov 2024 11:33:44 +0100 Subject: [PATCH 21/71] modifying rust_quad_ker_qcd, PyQuadKerQCDT and QuadQCDargs to include qed case --- crates/eko/src/lib.rs | 42 +++++++++++++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 11 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index f8c7cb1d8..3e383ce92 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -57,11 +57,7 @@ fn unravel_qed_singlet( /// Map tensors to c-ordered list in the QED singlet and valence case /// (res is a matrix with dim order_qcd x order_qed filled with complex numbers) -fn unravel_qed_ns( - res: Vec>>, - order_qcd: usize, - order_qed: usize, -) -> RawCmplx { +fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usize) -> RawCmplx { let mut target = RawCmplx { re: Vec::::new(), im: Vec::::new(), @@ -115,13 +111,30 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { ); } } else if is_singlet { - let gamma_singlet_qcd = match args.is_polarized { - true => ekore::anomalous_dimensions::polarized::spacelike::gamma_singlet_qcd, - false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd, - }; - raw = unravel( - gamma_singlet_qcd(args.order_qcd, &mut c, args.nf), + if args.is_qed { + let gamma_singlet_qed = + ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qed; + raw = unravel_qed_singlet( + gamma_singlet_qed(args.order_qcd, args.order_qed, &mut c, args.nf), + args.order_qcd, + args.order_qed, + ); + } else { + let gamma_singlet_qcd = match args.is_polarized { + true => ekore::anomalous_dimensions::polarized::spacelike::gamma_singlet_qcd, + false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd, + }; + raw = unravel( + gamma_singlet_qcd(args.order_qcd, &mut c, args.nf), + args.order_qcd, + ); + } + } else if args.is_qed { + let gamma_ns_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qed; + raw = unravel_qed_ns( + gamma_ns_qed(args.order_qcd, args.order_qed, args.mode0, &mut c, args.nf), args.order_qcd, + args.order_qed, ); } else { // we can not do 1D @@ -145,6 +158,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { jac.re, jac.im, args.order_qcd, + args.order_qed, is_singlet, args.mode0, args.mode1, @@ -175,6 +189,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn( f64, f64, usize, + usize, bool, u16, u16, @@ -201,6 +216,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn( #[derive(Clone, Copy)] pub struct QuadQCDargs { pub order_qcd: usize, + pub order_qed: usize, pub mode0: u16, pub mode1: u16, pub is_polarized: bool, @@ -221,6 +237,7 @@ pub struct QuadQCDargs { pub sv_mode_num: u8, pub is_threshold: bool, pub is_ome: bool, + pub is_qed: bool, pub Lsv: f64, } @@ -236,6 +253,7 @@ pub unsafe extern "C" fn my_py( _re_jac: f64, _im_jac: f64, _order_qcd: usize, + _order_qed: usize, _is_singlet: bool, _mode0: u16, _mode1: u16, @@ -269,6 +287,7 @@ pub unsafe extern "C" fn my_py( pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { QuadQCDargs { order_qcd: 0, + order_qed: 0, mode0: 0, mode1: 0, is_polarized: false, @@ -289,6 +308,7 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { sv_mode_num: 0, is_threshold: false, is_ome: false, + is_qed: false, Lsv: 0., } } From 588496ada1b781f6a7c5a2e29ee7b7c68c64053d Mon Sep 17 00:00:00 2001 From: tgiani Date: Thu, 21 Nov 2024 12:11:03 +0100 Subject: [PATCH 22/71] add qed valence option in rust_quad_ker_qcd --- crates/eko/src/lib.rs | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 3e383ce92..c241fb759 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -33,7 +33,7 @@ fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: us /// Map tensors to c-ordered list in the QED singlet and valence case /// (res is a matrix with dim order_qcd x order_qed filled with DIMxDIM matrices) -fn unravel_qed_singlet( +fn unravel_qed( res: Vec; DIM]; DIM]>>, order_qcd: usize, order_qed: usize, @@ -55,7 +55,7 @@ fn unravel_qed_singlet( target } -/// Map tensors to c-ordered list in the QED singlet and valence case +/// Map tensors to c-ordered list in the QED non-singlet case /// (res is a matrix with dim order_qcd x order_qed filled with complex numbers) fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usize) -> RawCmplx { let mut target = RawCmplx { @@ -79,6 +79,7 @@ fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usiz 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); + let is_qed_valence = (10200 == args.mode0) || (10204 == args.mode0); // prepare Mellin stuff let path = mellin::TalbotPath::new(u, args.logx, is_singlet); let jac = path.jac() * path.prefactor(); @@ -114,7 +115,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { if args.is_qed { let gamma_singlet_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qed; - raw = unravel_qed_singlet( + raw = unravel_qed( gamma_singlet_qed(args.order_qcd, args.order_qed, &mut c, args.nf), args.order_qcd, args.order_qed, @@ -130,12 +131,22 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { ); } } else if args.is_qed { - let gamma_ns_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qed; - raw = unravel_qed_ns( - gamma_ns_qed(args.order_qcd, args.order_qed, args.mode0, &mut c, args.nf), - args.order_qcd, - args.order_qed, - ); + if is_qed_valence { + let gamma_valence_qed = + ekore::anomalous_dimensions::unpolarized::spacelike::gamma_valence_qed; + raw = unravel_qed( + gamma_valence_qed(args.order_qcd, args.order_qed, &mut c, args.nf), + args.order_qcd, + args.order_qed, + ); + } else { + let gamma_ns_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qed; + raw = unravel_qed_ns( + gamma_ns_qed(args.order_qcd, args.order_qed, args.mode0, &mut c, args.nf), + args.order_qcd, + args.order_qed, + ); + } } else { // we can not do 1D let gamma_ns_qcd = match args.is_polarized { From 3452243a1d05acdc1307f35dc84446256da40f82 Mon Sep 17 00:00:00 2001 From: tgiani Date: Thu, 21 Nov 2024 14:29:44 +0100 Subject: [PATCH 23/71] remove useless flag is_qed --- crates/eko/src/lib.rs | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index c241fb759..040a12c69 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -112,7 +112,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { ); } } else if is_singlet { - if args.is_qed { + if args.order_qed > 0 { let gamma_singlet_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qed; raw = unravel_qed( @@ -130,7 +130,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.order_qcd, ); } - } else if args.is_qed { + } else if args.order_qed > 0 { if is_qed_valence { let gamma_valence_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_valence_qed; @@ -248,7 +248,6 @@ pub struct QuadQCDargs { pub sv_mode_num: u8, pub is_threshold: bool, pub is_ome: bool, - pub is_qed: bool, pub Lsv: f64, } @@ -319,7 +318,6 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { sv_mode_num: 0, is_threshold: false, is_ome: false, - is_qed: false, Lsv: 0., } } From 9335b97ed0901020ddd7c409d3d618383888058f Mon Sep 17 00:00:00 2001 From: tgiani Date: Thu, 21 Nov 2024 15:42:56 +0100 Subject: [PATCH 24/71] extend QuadQCDargs to include arguments for c_quad_ker_qed --- crates/eko/src/lib.rs | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 040a12c69..d6f0153ec 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -188,6 +188,12 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.sv_mode_num, args.is_threshold, args.Lsv, + // additional QED params + args.as_list, + args.mu2_from, + args.mu2_to, + args.a_half, + args.alphaem_running, ) } @@ -219,6 +225,11 @@ type PyQuadKerQCDT = unsafe extern "C" fn( u8, bool, f64, + *const f64, + f64, + f64, + f64, + bool, ) -> f64; /// Additional integration parameters @@ -249,6 +260,12 @@ pub struct QuadQCDargs { pub is_threshold: bool, pub is_ome: bool, pub Lsv: f64, + // additional param required for QED + pub as_list: *const f64, + pub mu2_from: f64, + pub mu2_to: f64, + pub a_half: f64, + pub alphaem_running: bool, } /// Empty placeholder function for python callback. @@ -282,6 +299,11 @@ pub unsafe extern "C" fn my_py( _sv_mode_num: u8, _is_threshold: bool, _lsv: f64, + _as_list: *const f64, + _mu2_from: f64, + _mu2_to: f64, + _a_half: f64, + _alphaem_running: bool, ) -> f64 { 0. } @@ -319,5 +341,10 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { is_threshold: false, is_ome: false, Lsv: 0., + as_list: [].as_ptr(), + mu2_from: 0., + mu2_to: 0., + a_half: 0., + alphaem_running: false, } } From 3390acaa22e07db638c513e5159b3b8181fedfa0 Mon Sep 17 00:00:00 2001 From: tgiani Date: Thu, 28 Nov 2024 16:35:23 +0100 Subject: [PATCH 25/71] setting up input vectors for cb_quad_ker_qed --- .pre-commit-config.yaml | 2 +- crates/eko/src/lib.rs | 23 ++- src/eko/evolution_operator/quad_ker.py | 241 +++++++++++++++++++------ 3 files changed, 205 insertions(+), 61 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 11ee1e506..8b88dfbb8 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -30,7 +30,7 @@ repos: # Run the formatter. - id: ruff-format - repo: https://github.com/PyCQA/docformatter - rev: v1.7.5 + rev: "master" hooks: - id: docformatter additional_dependencies: [tomli] diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index d6f0153ec..b25520fba 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -190,9 +190,12 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.Lsv, // additional QED params args.as_list, + args.as_list_len, args.mu2_from, args.mu2_to, args.a_half, + args.a_half_x, + args.a_half_y, args.alphaem_running, ) } @@ -226,9 +229,12 @@ type PyQuadKerQCDT = unsafe extern "C" fn( bool, f64, *const f64, + u8, f64, f64, - f64, + *const f64, + u8, + u8, bool, ) -> f64; @@ -262,9 +268,12 @@ pub struct QuadQCDargs { pub Lsv: f64, // additional param required for QED pub as_list: *const f64, + pub as_list_len: u8, pub mu2_from: f64, pub mu2_to: f64, - pub a_half: f64, + pub a_half: *const f64, + pub a_half_x: u8, + pub a_half_y: u8, pub alphaem_running: bool, } @@ -300,9 +309,12 @@ pub unsafe extern "C" fn my_py( _is_threshold: bool, _lsv: f64, _as_list: *const f64, + _as_list_len: u8, _mu2_from: f64, _mu2_to: f64, - _a_half: f64, + _a_half: *const f64, + _a_half_x: u8, + _a_half_y: u8, _alphaem_running: bool, ) -> f64 { 0. @@ -342,9 +354,12 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { is_ome: false, Lsv: 0., as_list: [].as_ptr(), + as_list_len: 0, mu2_from: 0., mu2_to: 0., - a_half: 0., + a_half: [].as_ptr(), + a_half_x: 0, + a_half_y: 0, alphaem_running: false, } } diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index 087f6bbb7..0ce15d587 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -47,6 +47,7 @@ def select_singlet_element(ker, mode0, mode1): nb.types.double, # re_jac nb.types.double, # im_jac nb.types.uintc, # order_qcd + nb.types.uintc, # order_qed nb.types.bool_, # is_singlet nb.types.uintc, # mode0 nb.types.uintc, # mode1 @@ -65,6 +66,14 @@ def select_singlet_element(ker, mode0, mode1): nb.types.uintc, # sv_mode_num nb.types.bool_, # is_threshold nb.types.double, # Lsv + nb.types.CPointer(nb.types.double), # as_list + nb.types.uintc, # as_list_len + nb.types.double, # mu2_from + nb.types.double, # mu2_to + nb.types.CPointer(nb.types.double), # a_half + nb.types.uintc, # a_half_x + nb.types.uintc, # a_half_y + nb.types.bool_, # alphaem_running ) @@ -81,6 +90,7 @@ def cb_quad_ker_qcd( re_jac, im_jac, order_qcd, + _order_qed, is_singlet, mode0, mode1, @@ -99,6 +109,14 @@ def cb_quad_ker_qcd( sv_mode, is_threshold, Lsv, + _as_list, + _as_list_len, + _mu2_from, + _mu2_to, + _a_half, + _a_half_x, + _a_half_y, + _alphaem_running, ): """C Callback inside integration kernel.""" # recover complex variables @@ -225,6 +243,7 @@ def cb_quad_ker_ome( re_jac, im_jac, order_qcd, + _order_qed, is_singlet, mode0, mode1, @@ -243,6 +262,14 @@ def cb_quad_ker_ome( sv_mode, _is_threshold, Lsv, + _as_list, + _as_list_len, + _mu2_from, + _mu2_to, + _a_half, + _a_half_x, + _a_half_y, + _alphaem_running, ): """C Callback inside integration kernel.""" # recover complex variables @@ -280,68 +307,170 @@ def cb_quad_ker_ome( 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. -# @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] -# 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. -# @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.cfunc( +# CB_SIGNATURE, +# cache=True, +# nopython=True, +# ) +# def cb_quad_ker_qed( +# re_gamma_raw, +# im_gamma_raw, +# re_n, +# im_n, +# re_jac, +# im_jac, +# order_qcd, +# order_qed, +# is_singlet, +# mode0, +# mode1, +# nf, +# is_log, +# logx, +# areas_raw, +# areas_x, +# areas_y, +# _L, +# ev_method, +# as1, +# as0, +# ev_op_iterations, +# ev_op_max_order_qcd, +# sv_mode, +# is_threshold, +# Lsv, +# as_list, +# as_list_len, +# mu2_from, +# mu2_to, +# a_half, +# a_half_x, +# a_half_y, +# alphaem_running, +# ): +# """C Callback inside integration kernel.""" +# recover complex variables +# n = re_n + im_n * 1j +# jac = re_jac + im_jac * 1j +# compute basis functions +# areas = nb.carray(areas_raw, (areas_x, areas_y)) +# pj = interpolation.evaluate_grid(n, is_log, logx, areas) +# order = (order_qcd, order_qed) +# ev_op_max_order = (ev_op_max_order_qcd, order_qed) +# is_valence = (mode0 == 10200) or (mode0 == 10204) + +# if is_singlet: +# reconstruct singlet matrices +# re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, order_qed, 4, 4)) +# im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, order_qed, 4, 4)) +# gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j + +# scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_singlet = sv.exponentiated.gamma_variation_qed( +# gamma_singlet, order, nf, L, alphaem_running +# ) -# 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] +# ker = qed_s.dispatcher( +# order, +# method, +# gamma_s, +# as_list, ### +# a_half, ### +# nf, +# ev_op_iterations, +# ev_op_max_order, +# ) +# 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 is_valence: + +# 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, Lsv) +# construct eko +# ker = ns.dispatcher( +# order, +# ev_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, Lsv) * ker +# recombine everything +# res = ker * pj * jac +# return np.real(res) # @nb.njit(cache=True) # def quad_ker_qed( From 6a3af197636f2cfebeb2ddbe3374da848d2b4b44 Mon Sep 17 00:00:00 2001 From: tgiani Date: Tue, 3 Dec 2024 17:16:22 +0100 Subject: [PATCH 26/71] cb_quad_ker_qed --- src/eko/evolution_operator/quad_ker.py | 413 ++++++++++--------------- 1 file changed, 162 insertions(+), 251 deletions(-) diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index 0ce15d587..aaedb7374 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -9,7 +9,11 @@ from .. import scale_variations as sv from ..io.types import InversionMethod from ..kernels import non_singlet as ns +from ..kernels import non_singlet_qed as qed_ns from ..kernels import singlet as s +from ..kernels import singlet_qed as qed_s +from ..kernels import valence_qed as qed_v +from ..matchings import lepton_number from ..scale_variations import expanded as sv_expanded from ..scale_variations import exponentiated as sv_exponentiated @@ -365,254 +369,161 @@ def select_QEDvalence_element(ker, mode0, mode1): return ker[index1, index2] -# @nb.cfunc( -# CB_SIGNATURE, -# cache=True, -# nopython=True, -# ) -# def cb_quad_ker_qed( -# re_gamma_raw, -# im_gamma_raw, -# re_n, -# im_n, -# re_jac, -# im_jac, -# order_qcd, -# order_qed, -# is_singlet, -# mode0, -# mode1, -# nf, -# is_log, -# logx, -# areas_raw, -# areas_x, -# areas_y, -# _L, -# ev_method, -# as1, -# as0, -# ev_op_iterations, -# ev_op_max_order_qcd, -# sv_mode, -# is_threshold, -# Lsv, -# as_list, -# as_list_len, -# mu2_from, -# mu2_to, -# a_half, -# a_half_x, -# a_half_y, -# alphaem_running, -# ): -# """C Callback inside integration kernel.""" -# recover complex variables -# n = re_n + im_n * 1j -# jac = re_jac + im_jac * 1j -# compute basis functions -# areas = nb.carray(areas_raw, (areas_x, areas_y)) -# pj = interpolation.evaluate_grid(n, is_log, logx, areas) -# order = (order_qcd, order_qed) -# ev_op_max_order = (ev_op_max_order_qcd, order_qed) -# is_valence = (mode0 == 10200) or (mode0 == 10204) - -# if is_singlet: -# reconstruct singlet matrices -# re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, order_qed, 4, 4)) -# im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, order_qed, 4, 4)) -# gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j - -# scale var exponentiated is directly applied on gamma -# if sv_mode == sv.Modes.exponentiated: -# gamma_singlet = sv.exponentiated.gamma_variation_qed( -# gamma_singlet, 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, -# ) -# 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 is_valence: - -# 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, Lsv) -# construct eko -# ker = ns.dispatcher( -# order, -# ev_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, Lsv) * ker -# recombine everything -# res = ker * pj * jac -# return np.real(res) - -# @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 +@nb.cfunc( + CB_SIGNATURE, + cache=True, + nopython=True, +) +def cb_quad_ker_qed( + re_gamma_raw, + im_gamma_raw, + re_n, + im_n, + re_jac, + im_jac, + order_qcd, + order_qed, + is_singlet, + mode0, + mode1, + nf, + is_log, + logx, + areas_raw, + areas_x, + areas_y, + L, + ev_method, + as1, + as0, + ev_op_iterations, + ev_op_max_order_qcd, + sv_mode, + is_threshold, + Lsv, + as_list_raw, + as_list_len, + mu2_from, + mu2_to, + a_half_raw, + a_half_x, + a_half_y, + alphaem_running, +): + """C Callback inside integration kernel.""" + # recover complex variables + n = re_n + im_n * 1j + jac = re_jac + im_jac * 1j + # compute basis functions + areas = nb.carray(areas_raw, (areas_x, areas_y)) + pj = interpolation.evaluate_grid(n, is_log, logx, areas) + order = (order_qcd, order_qed) + ev_op_max_order = (ev_op_max_order_qcd, order_qed) + is_valence = (mode0 == 10200) or (mode0 == 10204) + + as_list = nb.carray(as_list_raw, as_list_len) + a_half = nb.carray(a_half_raw, (a_half_x, a_half_y)) + + if is_singlet: + # reconstruct singlet matrices + re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, order_qed, 4, 4)) + im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, order_qed, 4, 4)) + gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j + + # scale var exponentiated is directly applied on gamma + if sv_mode == sv.Modes.exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation_qed( + gamma_singlet, order, nf, lepton_number(mu2_to), L, alphaem_running + ) + + ker = qed_s.dispatcher( + order, + ev_method, + gamma_singlet, + as_list, + a_half, + nf, + ev_op_iterations, + ev_op_max_order, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = np.ascontiguousarray( + sv.expanded.singlet_variation_qed( + gamma_singlet, + as_list[-1], + a_half[-1][1], + alphaem_running, + order, + nf, + L, + ) + ) @ np.ascontiguousarray(ker) + ker = select_QEDsinglet_element(ker, mode0, mode1) + + elif is_valence: + # reconstruct valence matrices + re_gamma_valence = nb.carray(re_gamma_raw, (order_qcd, order_qed, 2, 2)) + im_gamma_valence = nb.carray(im_gamma_raw, (order_qcd, order_qed, 2, 2)) + gamma_valence = re_gamma_valence + im_gamma_valence * 1j + + if sv_mode == sv.Modes.exponentiated: + gamma_valence = sv.exponentiated.gamma_variation_qed( + gamma_valence, order, nf, lepton_number(mu2_to), L, alphaem_running + ) + ker = qed_v.dispatcher( + order, + ev_method, + gamma_valence, + 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_valence, + as_list[-1], + a_half[-1][1], + alphaem_running, + order, + nf, + L, + ) + ) @ np.ascontiguousarray(ker) + ker = select_QEDvalence_element(ker, mode0, mode1) + + else: + # construct non-singlet matrices + re_gamma_ns = nb.carray(re_gamma_raw, (order_qcd, order_qed)) + im_gamma_ns = nb.carray(im_gamma_raw, (order_qcd, order_qed)) + gamma_ns = re_gamma_ns + im_gamma_ns * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_ns = sv_exponentiated.gamma_variation_qed( + gamma_ns, order, nf, lepton_number(mu2_to), L, alphaem_running + ) + # construct eko + ker = qed_ns.dispatcher( + order, + ev_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 + ) + # recombine everything + res = ker * pj * jac + return np.real(res) From f0a727dd55c645cfad176ab5f63c0f92414e32b8 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 9 Dec 2024 13:52:07 +0100 Subject: [PATCH 27/71] fix modes for singlet QED and use == in if statements --- crates/eko/src/lib.rs | 2 +- src/eko/kernels/non_singlet.py | 2 +- src/eko/kernels/singlet.py | 8 ++++---- src/eko/kernels/singlet_qed.py | 2 +- src/eko/kernels/valence_qed.py | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index b25520fba..f8aac3047 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -78,7 +78,7 @@ fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usiz #[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); + let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0)|| (22 == args.mode0) || (101 == args.mode0); let is_qed_valence = (10200 == args.mode0) || (10204 == args.mode0); // prepare Mellin stuff let path = mellin::TalbotPath::new(u, args.logx, is_singlet); diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py index 1d25d5845..0e8352dda 100644 --- a/src/eko/kernels/non_singlet.py +++ b/src/eko/kernels/non_singlet.py @@ -375,7 +375,7 @@ def dispatcher(order, method, gamma_ns, a1, a0, nf, ev_op_iterations): # pylint return eko_ordered_truncated( gamma_ns, a1, a0, betalist, order, ev_op_iterations ) - if method is EvoMethods.TRUNCATED: + if method == EvoMethods.TRUNCATED: return eko_truncated(gamma_ns, a1, a0, betalist, order, ev_op_iterations) # NLO diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py index 96f02a972..2ce7bce21 100644 --- a/src/eko/kernels/singlet.py +++ b/src/eko/kernels/singlet.py @@ -613,7 +613,7 @@ def dispatcher( # pylint: disable=too-many-return-statements # Common method for NLO and NNLO if method in [EvoMethods.ITERATE_EXACT, EvoMethods.ITERATE_EXPANDED]: return eko_iterate(gamma_singlet, a1, a0, betalist, order, ev_op_iterations) - if method is EvoMethods.PERTURBATIVE_EXACT: + if method == EvoMethods.PERTURBATIVE_EXACT: return eko_perturbative( gamma_singlet, a1, @@ -624,7 +624,7 @@ def dispatcher( # pylint: disable=too-many-return-statements ev_op_max_order, True, ) - if method is EvoMethods.PERTURBATIVE_EXPANDED: + if method == EvoMethods.PERTURBATIVE_EXPANDED: return eko_perturbative( gamma_singlet, a1, @@ -638,13 +638,13 @@ def dispatcher( # pylint: disable=too-many-return-statements if method in [EvoMethods.TRUNCATED, EvoMethods.ORDERED_TRUNCATED]: return eko_truncated(gamma_singlet, a1, a0, betalist, order, ev_op_iterations) # These methods are scattered for nlo and nnlo - if method is EvoMethods.DECOMPOSE_EXACT: + if method == EvoMethods.DECOMPOSE_EXACT: if order[0] == 2: return nlo_decompose_exact(gamma_singlet, a1, a0, betalist) if order[0] == 3: return nnlo_decompose_exact(gamma_singlet, a1, a0, betalist) return n3lo_decompose_exact(gamma_singlet, a1, a0, nf) - if method is EvoMethods.DECOMPOSE_EXPANDED: + if method == EvoMethods.DECOMPOSE_EXPANDED: if order[0] == 2: return nlo_decompose_expanded(gamma_singlet, a1, a0, betalist) if order[0] == 3: diff --git a/src/eko/kernels/singlet_qed.py b/src/eko/kernels/singlet_qed.py index d63a2f1ac..0a5716df8 100644 --- a/src/eko/kernels/singlet_qed.py +++ b/src/eko/kernels/singlet_qed.py @@ -97,7 +97,7 @@ def dispatcher( e_s : numpy.ndarray singlet EKO """ - if method is EvoMethods.ITERATE_EXACT: + if method == EvoMethods.ITERATE_EXACT: return eko_iterate( gamma_singlet, as_list, a_half, nf, order, ev_op_iterations, 4 ) diff --git a/src/eko/kernels/valence_qed.py b/src/eko/kernels/valence_qed.py index 8b83e1917..b5da158b4 100644 --- a/src/eko/kernels/valence_qed.py +++ b/src/eko/kernels/valence_qed.py @@ -45,7 +45,7 @@ def dispatcher( e_v : numpy.ndarray singlet EKO """ - if method is EvoMethods.ITERATE_EXACT: + if method == EvoMethods.ITERATE_EXACT: return eko_iterate( gamma_valence, as_list, a_half, nf, order, ev_op_iterations, 2 ) From 398c89aa9db530d35bcccbf9b01a60fb3b9662e0 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 9 Dec 2024 14:21:09 +0100 Subject: [PATCH 28/71] small fix --- crates/eko/src/lib.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index f8aac3047..44fbc14fb 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -78,6 +78,7 @@ fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usiz #[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)|| (22 == args.mode0) || (101 == args.mode0); let is_qed_valence = (10200 == args.mode0) || (10204 == args.mode0); // prepare Mellin stuff From ea44cee7a0b962e492161b44d8f62ff3a43e0e75 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 9 Dec 2024 14:22:43 +0100 Subject: [PATCH 29/71] cargo fmt --- crates/eko/src/lib.rs | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 44fbc14fb..3aa6d1df8 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -78,8 +78,12 @@ fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usiz #[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)|| (22 == args.mode0) || (101 == args.mode0); + + let is_singlet = (100 == args.mode0) + || (21 == args.mode0) + || (90 == args.mode0) + || (22 == args.mode0) + || (101 == args.mode0); let is_qed_valence = (10200 == args.mode0) || (10204 == args.mode0); // prepare Mellin stuff let path = mellin::TalbotPath::new(u, args.logx, is_singlet); From 57213266a7e4fda5279454fce5d0476c7bccc9d5 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 9 Dec 2024 14:30:51 +0100 Subject: [PATCH 30/71] forgot pre-commit --- crates/eko/src/lib.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 3aa6d1df8..ea42fb1e9 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -84,6 +84,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { || (90 == args.mode0) || (22 == args.mode0) || (101 == args.mode0); + let is_qed_valence = (10200 == args.mode0) || (10204 == args.mode0); // prepare Mellin stuff let path = mellin::TalbotPath::new(u, args.logx, is_singlet); From 83af8361d52bc37832c8078fba3d52af5859f13a Mon Sep 17 00:00:00 2001 From: tgiani Date: Wed, 11 Dec 2024 11:52:26 +0100 Subject: [PATCH 31/71] updating patch file --- src/eko/evolution_operator/__init__.py.patch | 41 ++++++++++++++++---- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index a4951e687..41ddb9f9d 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,5 +1,5 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index bd1b19d6..f543f7bc 100644 +index bd1b19d6..de87651c 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -3,16 +3,16 @@ r"""Contains the central operator classes. @@ -21,11 +21,12 @@ index bd1b19d6..f543f7bc 100644 import ekore.anomalous_dimensions.polarized.space_like as ad_ps import ekore.anomalous_dimensions.unpolarized.space_like as ad_us -@@ -32,91 +32,10 @@ from ..matchings import Segment, lepton_number +@@ -32,91 +32,11 @@ from ..matchings import Segment, lepton_number from ..member import OpMember from ..scale_variations import expanded as sv_expanded from ..scale_variations import exponentiated as sv_exponentiated +from .quad_ker import cb_quad_ker_qcd ++from .quad_ker import cb_quad_ker_qed logger = logging.getLogger(__name__) @@ -114,7 +115,7 @@ index bd1b19d6..f543f7bc 100644 spec = [ ("is_singlet", nb.boolean), ("is_QEDsinglet", nb.boolean), -@@ -188,421 +107,6 @@ class QuadKerBase: +@@ -188,422 +108,6 @@ class QuadKerBase: return self.path.prefactor * pj * self.path.jac @@ -533,9 +534,10 @@ index bd1b19d6..f543f7bc 100644 - ) - return ker - - +- OpMembers = Dict[OperatorLabel, OpMember] """Map of all operators.""" + @@ -792,49 +296,6 @@ class Operator(sv.ScaleVariationModeMixin): """Return the evolution method.""" return ev_method(EvolutionMethod(self.config["method"])) @@ -586,7 +588,7 @@ index bd1b19d6..f543f7bc 100644 def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( -@@ -857,10 +318,14 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -857,10 +318,17 @@ class Operator(sv.ScaleVariationModeMixin): else: self.op_members[n] = zero.copy() @@ -598,14 +600,17 @@ index bd1b19d6..f543f7bc 100644 + """Adjust integration config.""" + cfg.as1 = self.as_list[1] + cfg.as0 = self.as_list[0] -+ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qcd.address) ++ if self.order[1] == 0: ++ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qcd.address) ++ else: ++ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qed.address) + cfg.method_num = self.ev_method + + def run_op_integration(self, log_grid): """Run the integration for each grid point. Parameters -@@ -875,18 +339,53 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -875,18 +343,75 @@ class Operator(sv.ScaleVariationModeMixin): """ column = [] k, logx = log_grid @@ -615,6 +620,7 @@ index bd1b19d6..f543f7bc 100644 + # start preparing C arguments + cfg = ekors.lib.empty_qcd_args() + cfg.order_qcd = self.order[0] ++ cfg.order_qed = self.order[1] + cfg.is_polarized = self.config["polarized"] + cfg.is_time_like = self.config["time_like"] + cfg.nf = self.nf @@ -625,6 +631,27 @@ index bd1b19d6..f543f7bc 100644 + cfg.ev_op_max_order_qcd = self.config["ev_op_max_order"][0] + cfg.sv_mode_num = self.sv_mode + cfg.is_threshold = self.is_threshold ++ cfg.mu2_from = self.q2_from ++ cfg.mu2_to = self.q2_to ++ cfg.alphaem_running=self.alphaem_running ++ ++ # prepare as_list for c ++ as_list_len = self.as_list.shape[0] ++ as_list_ffi = ekors.ffi.new(f"double[{as_list_len}]", self.as_list.tolist()) ++ cfg.as_list = as_list_ffi ++ cfg.as_list_len = as_list_len ++ ++ # prepare a_half for c ++ a_half_x = self.a_half_list.shape[0] ++ a_half_y = self.a_half_list.shape[1] ++ a_half_len = a_half_x * a_half_y ++ a_half_ffi = ekors.ffi.new( ++ f"double[{a_half_len}]", self.a_half_list.flatten().tolist() ++ ) ++ cfg.a_half = a_half_ffi ++ cfg.a_half_x = a_half_x ++ cfg.a_half_y = a_half_y ++ + self.update_cfg(cfg) + # iterate basis functions From 1af5cbf4739f3de81785a44e47fec322c3aa92d5 Mon Sep 17 00:00:00 2001 From: tgiani <33056186+tgiani@users.noreply.github.com> Date: Wed, 11 Dec 2024 11:55:01 +0100 Subject: [PATCH 32/71] Update crates/ekore/src/constants.rs Co-authored-by: Felix Hekhorn --- crates/ekore/src/constants.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 80c3bc807..61cc0d9e0 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -76,7 +76,7 @@ pub struct ChargeCombinations { impl ChargeCombinations { pub fn nu(&self) -> u8 { - uplike_flavors(self.nf) + self.nf / 2 } pub fn nd(&self) -> u8 { From 5ce8eb3b51257d272e4acd9fad580bee1b3af1c5 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 14:31:47 +0200 Subject: [PATCH 33/71] rust: Iterate ekore/constants --- .../unpolarized/spacelike.rs | 14 ++++---- .../unpolarized/spacelike/aem1.rs | 14 ++++---- .../unpolarized/spacelike/as1aem1.rs | 32 +++++++++---------- crates/ekore/src/constants.rs | 31 +++++++++--------- 4 files changed, 45 insertions(+), 46 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 38479b89f..179ed0686 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -1,7 +1,7 @@ //! The unpolarized, space-like anomalous dimensions at various couplings power. use crate::constants::{ - ed2, eu2, PID_NSM, PID_NSM_ED2, PID_NSM_EU2, PID_NSP, PID_NSP_ED2, PID_NSP_EU2, PID_NSV, + ED2, EU2, PID_NSM, PID_NSM_ED2, PID_NSM_EU2, PID_NSP, PID_NSP_ED2, PID_NSP_EU2, PID_NSV, }; use crate::harmonics::cache::Cache; use num::complex::Complex; @@ -78,18 +78,18 @@ pub fn gamma_ns_qed( pub fn choose_ns_as_aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { match mode { - PID_NSP_EU2 | PID_NSM_EU2 => eu2 * aem1::gamma_ns(c, nf), - PID_NSP_ED2 | PID_NSM_ED2 => ed2 * aem1::gamma_ns(c, nf), + PID_NSP_EU2 | PID_NSM_EU2 => EU2 * aem1::gamma_ns(c, nf), + PID_NSP_ED2 | PID_NSM_ED2 => ED2 * aem1::gamma_ns(c, nf), _ => panic!("Unkown non-singlet sector element"), } } pub fn choose_ns_as_as1aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { match mode { - PID_NSP_EU2 => eu2 * as1aem1::gamma_nsp(c, nf), - PID_NSP_ED2 => ed2 * as1aem1::gamma_nsp(c, nf), - PID_NSM_EU2 => eu2 * as1aem1::gamma_nsm(c, nf), - PID_NSM_ED2 => ed2 * as1aem1::gamma_nsm(c, nf), + PID_NSP_EU2 => EU2 * as1aem1::gamma_nsp(c, nf), + PID_NSP_ED2 => ED2 * as1aem1::gamma_nsp(c, nf), + PID_NSM_EU2 => EU2 * as1aem1::gamma_nsm(c, nf), + PID_NSM_ED2 => ED2 * as1aem1::gamma_nsm(c, nf), _ => panic!("Unkown non-singlet sector element"), } } diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index bc70ed97f..0a5feaa6b 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -2,7 +2,7 @@ use num::complex::Complex; use num::Zero; -use crate::constants::{ed2, eu2, uplike_flavors, ChargeCombinations, CF, NC, TR}; +use crate::constants::{ChargeCombinations, CF, ED2, EU2, NC, TR}; use crate::harmonics::cache::Cache; use crate::anomalous_dimensions::unpolarized::spacelike::as1; @@ -25,9 +25,8 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (2.5) of pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { - let nu = uplike_flavors(nf); - let nd = nf - nu; - (4.0 / 3.0 * (NC as f64) * ((nu as f64) * eu2 + (nd as f64) * ed2)).into() + let cc = ChargeCombinations { nf }; + (4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2)).into() } /// Compute the leading-order non-singlet QED anomalous dimension @@ -117,11 +116,12 @@ mod tests { let mut c = Cache::new(N); for nf in 2..7 { - let nu = uplike_flavors(nf); - let nd = nf - nu; + let cc = ChargeCombinations { nf }; assert_approx_eq_cmplx!( f64, - eu2 * gamma_qph(&mut c, nu) + ed2 * gamma_qph(&mut c, nd) + gamma_phph(&mut c, nf), + EU2 * gamma_qph(&mut c, cc.nu()) + + ED2 * gamma_qph(&mut c, cc.nd()) + + gamma_phph(&mut c, nf), cmplx!(0., 0.), epsilon = 2e-6 ); diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 546acca85..48220e8c3 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -2,9 +2,7 @@ use crate::cmplx; use num::complex::Complex; -use crate::constants::{ - ed2, eu2, uplike_flavors, ChargeCombinations, CA, CF, NC, TR, ZETA2, ZETA3, -}; +use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC, TR, ZETA2, ZETA3}; use crate::harmonics::cache::{Cache, K}; /// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. @@ -100,9 +98,11 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (28) of pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { - let nu = uplike_flavors(nf); - let nd = nf - nu; - cmplx!(4.0 * CF * CA * ((nu as f64) * eu2 + (nd as f64) * ed2), 0.) + let cc = ChargeCombinations { nf }; + cmplx!( + 4.0 * CF * CA * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2), + 0. + ) } /// Compute the $O(a_s^1a_{em}^1)$ gluon-gluon singlet anomalous dimension. @@ -261,14 +261,13 @@ mod tests { let mut c = Cache::new(N); for nf in 2..7 { - let nu = uplike_flavors(nf); - let nd = nf - nu; + let cc = ChargeCombinations { nf }; assert_approx_eq_cmplx!( f64, - eu2 * gamma_qg(&mut c, nu) - + ed2 * gamma_qg(&mut c, nd) - + (nu as f64 * eu2 + nd as f64 * ed2) * gamma_phg(&mut c, nf) - + (nu as f64 * eu2 + nd as f64 * ed2) * gamma_gg(&mut c, nf), + EU2 * gamma_qg(&mut c, cc.nu()) + + ED2 * gamma_qg(&mut c, cc.nd()) + + (cc.nu() as f64 * EU2 + cc.nd() as f64 * ED2) * gamma_phg(&mut c, nf) + + (cc.nu() as f64 * EU2 + cc.nd() as f64 * ED2) * gamma_gg(&mut c, nf), cmplx!(0., 0.), epsilon = 1e-14 ); @@ -281,14 +280,13 @@ mod tests { let mut c = Cache::new(N); for nf in 2..7 { - let nu = uplike_flavors(nf); - let nd = nf - nu; + let cc = ChargeCombinations { nf }; assert_approx_eq_cmplx!( f64, - eu2 * gamma_qph(&mut c, nu) - + ed2 * gamma_qph(&mut c, nd) + EU2 * gamma_qph(&mut c, cc.nu()) + + ED2 * gamma_qph(&mut c, cc.nd()) + gamma_phph(&mut c, nf) - + (nu as f64 * eu2 + nd as f64 * ed2) * gamma_gph(&mut c, nf), + + (cc.nu() as f64 * EU2 + cc.nd() as f64 * ED2) * gamma_gph(&mut c, nf), cmplx!(0., 0.), epsilon = 1e-14 ); diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 61cc0d9e0..43da87af8 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -1,5 +1,4 @@ //! Global constants. -use std::unimplemented; /// The number of colors. /// @@ -24,12 +23,12 @@ pub const CF: f64 = ((NC * NC - 1) as f64) / ((2 * NC) as f64); /// Up quark charge square. /// /// Defaults to $e_u^2 = 4./9$ -pub const eu2: f64 = 4. / 9.; +pub const EU2: f64 = 4. / 9.; /// Down quark charge square. /// /// Defaults to $e_d^2 = 1./9$ -pub const ed2: f64 = 1. / 9.; +pub const ED2: f64 = 1. / 9.; /// Riemann zeta function at z = 2. /// @@ -53,48 +52,50 @@ pub const PID_NSM: u16 = 10201; /// non-singlet all-valence |PID|. pub const PID_NSV: u16 = 10200; -/// QED |PID|. Need to give sensible names +/// singlet-like non-singlet up-sector |PID| pub const PID_NSP_EU2: u16 = 10102; +/// singlet-like non-singlet down-sector |PID| pub const PID_NSP_ED2: u16 = 10103; +/// valence-like non-singlet up-sector |PID| pub const PID_NSM_EU2: u16 = 10202; +/// valence-like non-singlet down-sector |PID| pub const PID_NSM_ED2: u16 = 10203; -/// compute the number of up flavors -pub fn uplike_flavors(nf: u8) -> u8 { - if nf > 6 { - unimplemented!("Selected nf is not implemented") - } - nf / 2 -} - +/// |QED| electric charge combinations. pub struct ChargeCombinations { pub nf: u8, } impl ChargeCombinations { + /// Number of up-like flavors. pub fn nu(&self) -> u8 { self.nf / 2 } + /// Number of down-like flavors. pub fn nd(&self) -> u8 { self.nf - self.nu() } + /// Electric charge average. pub fn e2avg(&self) -> f64 { - (self.nu() as f64 * eu2 + self.nd() as f64 * ed2) / (self.nf as f64) + (self.nu() as f64 * EU2 + self.nd() as f64 * ED2) / (self.nf as f64) } + /// Relative up contribution to charge difference. pub fn vue2m(&self) -> f64 { - self.nu() as f64 / (self.nf as f64) * (eu2 - ed2) + self.nu() as f64 / (self.nf as f64) * (EU2 - ED2) } + /// Relative down contribution to charge difference. pub fn vde2m(&self) -> f64 { - self.nd() as f64 / (self.nf as f64) * (eu2 - ed2) + self.nd() as f64 / (self.nf as f64) * (EU2 - ED2) } + /// Asymmetric charge combination. pub fn e2delta(&self) -> f64 { self.vde2m() - self.vue2m() + self.e2avg() } From 88d7ed9822c6544542dc32ee907129cc3ed71e90 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 14:33:11 +0200 Subject: [PATCH 34/71] rust: Fix clippy warning --- crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 179ed0686..b310289c3 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -139,7 +139,6 @@ mod tests { use super::*; use crate::{assert_approx_eq_cmplx, cmplx}; use num::complex::Complex; - use num::Zero; #[test] fn gamma_ns() { From 979719066f570aff96306ccf9e3cc7e269b5326e Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 14:40:00 +0200 Subject: [PATCH 35/71] rust: Add missing ref for polarized ad --- crates/ekore/refs.bib | 12 ++++++++++++ crates/ekore/src/bib.rs | 14 +++++++++++++- crates/ekore/src/constants.rs | 16 ++++++++-------- 3 files changed, 33 insertions(+), 9 deletions(-) diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib index 0035b7a39..8d155ed48 100644 --- a/crates/ekore/refs.bib +++ b/crates/ekore/refs.bib @@ -107,3 +107,15 @@ @article{Bierenbaum2009zt pages = "401--406", year = "2009" } +@article{Gluck1995yr, + author = "Gluck, M. and Reya, E. and Stratmann, M. and Vogelsang, W.", + title = "{Next-to-leading order radiative parton model analysis of polarized deep inelastic lepton - nucleon scattering}", + eprint = "hep-ph/9508347", + archivePrefix = "arXiv", + reportNumber = "DO-TH-95-13, RAL-TR-95-042", + doi = "10.1103/PhysRevD.53.4775", + journal = "Phys. Rev. D", + volume = "53", + pages = "4775--4786", + year = "1996" +} diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs index 3e8a1eb8f..175efe194 100644 --- a/crates/ekore/src/bib.rs +++ b/crates/ekore/src/bib.rs @@ -1,4 +1,4 @@ -//! List of References (autogenerated on 2024-09-10T11:22:16.645053). +//! List of References (autogenerated on 2025-01-13T14:38:04.013494). #[allow(non_snake_case)] /// The Three loop splitting functions in QCD: The Nonsinglet case @@ -107,3 +107,15 @@ pub fn Buza1996wv() {} /// /// DOI: [10.1016/j.physletb.2009.01.057](https:dx.doi.org/10.1016/j.physletb.2009.01.057) pub fn Bierenbaum2009zt() {} + +#[allow(non_snake_case)] +/// Next-to-leading order radiative parton model analysis of polarized deep inelastic lepton - nucleon scattering +/// +/// Gluck, M. and Reya, E. and Stratmann, M. and Vogelsang, W. +/// +/// Published in: Phys. Rev. D 53 (1996), 4775--4786 +/// +/// e-Print: [hep-ph/9508347](https://arxiv.org/abs/hep-ph/9508347) +/// +/// DOI: [10.1103/PhysRevD.53.4775](https:dx.doi.org/10.1103/PhysRevD.53.4775) +pub fn Gluck1995yr() {} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 43da87af8..10c8cc2d6 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -20,14 +20,14 @@ pub const CA: f64 = NC as f64; /// 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); -/// Up quark charge square. +/// Up quark electric charge square. /// -/// Defaults to $e_u^2 = 4./9$ +/// Defaults to $e_u^2 = 4./9$. pub const EU2: f64 = 4. / 9.; -/// Down quark charge square. +/// Down quark electric charge square. /// -/// Defaults to $e_d^2 = 1./9$ +/// Defaults to $e_d^2 = 1./9$. pub const ED2: f64 = 1. / 9.; /// Riemann zeta function at z = 2. @@ -52,16 +52,16 @@ pub const PID_NSM: u16 = 10201; /// non-singlet all-valence |PID|. pub const PID_NSV: u16 = 10200; -/// singlet-like non-singlet up-sector |PID| +/// singlet-like non-singlet up-sector |PID|. pub const PID_NSP_EU2: u16 = 10102; -/// singlet-like non-singlet down-sector |PID| +/// singlet-like non-singlet down-sector |PID|. pub const PID_NSP_ED2: u16 = 10103; -/// valence-like non-singlet up-sector |PID| +/// valence-like non-singlet up-sector |PID|. pub const PID_NSM_EU2: u16 = 10202; -/// valence-like non-singlet down-sector |PID| +/// valence-like non-singlet down-sector |PID|. pub const PID_NSM_ED2: u16 = 10203; /// |QED| electric charge combinations. From 5013a8fd65b9000c3ef1644b84f5cc0f502c9ecc Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 14:48:22 +0200 Subject: [PATCH 36/71] rust: Iterate aem1 --- crates/ekore/refs.bib | 10 ++++++++++ .../unpolarized/spacelike/aem1.rs | 20 +++++++++---------- crates/ekore/src/bib.rs | 14 ++++++++++++- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib index 8d155ed48..480619d12 100644 --- a/crates/ekore/refs.bib +++ b/crates/ekore/refs.bib @@ -119,3 +119,13 @@ @article{Gluck1995yr pages = "4775--4786", year = "1996" } +@phdthesis{Carrazza2015dea, + author = "Carrazza, Stefano", + title = "{Parton distribution functions with QED corrections}", + eprint = "1509.00209", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + doi = "10.13130/carrazza-stefano_phd2015-07-06", + school = "Milan U.", + year = "2015" +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 0a5feaa6b..1a9e317ba 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -9,36 +9,36 @@ use crate::anomalous_dimensions::unpolarized::spacelike::as1; /// Compute the leading-order photon-quark anomalous dimension. /// -/// Implements Eq. (2.5) of +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex { as1::gamma_gq(c, nf) / CF } /// Compute the leading-order quark-photon anomalous dimension. /// -/// Implements Eq. (2.5) of +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. +/// However, we are adding the $N_C$ and the $2n_f$ factors from $\theta$ inside the +/// definition of $\gamma_{q \gamma}^{(0)}(N)$. pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { as1::gamma_qg(c, nf) / TR * (NC as f64) } /// Compute the leading-order photon-photon anomalous dimension. /// -/// Implements Eq. (2.5) of +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { let cc = ChargeCombinations { nf }; (4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2)).into() } -/// Compute the leading-order non-singlet QED anomalous dimension +/// Compute the leading-order non-singlet |QED| anomalous dimension. /// -/// Implements Eq. (2.5) of +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { as1::gamma_ns(c, nf) / CF } -/// Compute the leading-order singlet QED anomalous dimension matrix -/// -/// Implements Eq. (2.5) of +/// Compute the leading-order singlet |QED| anomalous dimension matrix. pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = ChargeCombinations { nf }; @@ -74,9 +74,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { ] } -/// Compute the leading-order valence QED anomalous dimension matrix -/// -/// Implements Eq. (2.5) of +/// Compute the leading-order valence |QED| anomalous dimension matrix. pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs index 175efe194..06533f6e2 100644 --- a/crates/ekore/src/bib.rs +++ b/crates/ekore/src/bib.rs @@ -1,4 +1,4 @@ -//! List of References (autogenerated on 2025-01-13T14:38:04.013494). +//! List of References (autogenerated on 2025-01-13T14:46:34.886936). #[allow(non_snake_case)] /// The Three loop splitting functions in QCD: The Nonsinglet case @@ -119,3 +119,15 @@ pub fn Bierenbaum2009zt() {} /// /// DOI: [10.1103/PhysRevD.53.4775](https:dx.doi.org/10.1103/PhysRevD.53.4775) pub fn Gluck1995yr() {} + +#[allow(non_snake_case)] +/// Parton distribution functions with QED corrections +/// +/// Carrazza, Stefano +/// +/// Published as PhD thesis at Milan U. (2015) +/// +/// e-Print: [1509.00209](https://arxiv.org/abs/1509.00209) +/// +/// DOI: [10.13130/carrazza-stefano_phd2015-07-06](https:dx.doi.org/10.13130/carrazza-stefano_phd2015-07-06) +pub fn Carrazza2015dea() {} From 064c31cf6f386a9a6fabca32fe4c294b5f547a81 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 15:05:55 +0200 Subject: [PATCH 37/71] rust: Iterate as1aem1 --- crates/ekore/refs.bib | 14 ++++++ .../unpolarized/spacelike/aem1.rs | 12 ++--- .../unpolarized/spacelike/as1.rs | 4 +- .../unpolarized/spacelike/as1aem1.rs | 44 +++++++++---------- crates/ekore/src/bib.rs | 14 +++++- crates/make_bib.py | 16 ++++--- 6 files changed, 68 insertions(+), 36 deletions(-) diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib index 480619d12..dc225226e 100644 --- a/crates/ekore/refs.bib +++ b/crates/ekore/refs.bib @@ -129,3 +129,17 @@ @phdthesis{Carrazza2015dea school = "Milan U.", year = "2015" } +@article{deFlorian2015ujt, + author = "de Florian, Daniel and Sborlini, Germ\'an F. R. and Rodrigo, Germ\'an", + title = "{QED corrections to the Altarelli\textendash{}Parisi splitting functions}", + eprint = "1512.00612", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "ICAS-03-15, IFIC-15-81", + doi = "10.1140/epjc/s10052-016-4131-8", + journal = "Eur. Phys. J. C", + volume = "76", + number = "5", + pages = "282", + year = "2016" +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 1a9e317ba..872674927 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -7,14 +7,14 @@ use crate::harmonics::cache::Cache; use crate::anomalous_dimensions::unpolarized::spacelike::as1; -/// Compute the leading-order photon-quark anomalous dimension. +/// Compute the photon-quark anomalous dimension. /// /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex { as1::gamma_gq(c, nf) / CF } -/// Compute the leading-order quark-photon anomalous dimension. +/// Compute the quark-photon anomalous dimension. /// /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. /// However, we are adding the $N_C$ and the $2n_f$ factors from $\theta$ inside the @@ -23,7 +23,7 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { as1::gamma_qg(c, nf) / TR * (NC as f64) } -/// Compute the leading-order photon-photon anomalous dimension. +/// Compute the photon-photon anomalous dimension. /// /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { @@ -31,14 +31,14 @@ pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { (4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2)).into() } -/// Compute the leading-order non-singlet |QED| anomalous dimension. +/// Compute the non-singlet anomalous dimension. /// /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { as1::gamma_ns(c, nf) / CF } -/// Compute the leading-order singlet |QED| anomalous dimension matrix. +/// Compute the singlet anomalous dimension matrix. pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = ChargeCombinations { nf }; @@ -74,7 +74,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { ] } -/// Compute the leading-order valence |QED| anomalous dimension matrix. +/// Compute the valence anomalous dimension matrix. pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs index bd2c67045..f0c41c3a7 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -53,7 +53,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { ] } -/// Compute the leading-order singlet anomalous dimension matrix +/// Compute the singlet anomalous dimension matrix /// for the unified evolution basis. pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { [ @@ -84,7 +84,7 @@ pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { ] } -/// Compute the leading-order valence anomalous dimension matrix +/// Compute the valence anomalous dimension matrix /// for the unified evolution basis. pub fn gamma_valence_qed(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { [ diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 48220e8c3..e8e9cc4bb 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -1,13 +1,13 @@ -//! The $O(a_s^1a_{em}^1)$ Altarelli-Parisi splitting kernels. +//! |LO| |QED| x |LO| |QCD|. use crate::cmplx; use num::complex::Complex; use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC, TR, ZETA2, ZETA3}; use crate::harmonics::cache::{Cache, K}; -/// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension. +/// Compute the photon-quark anomalous dimension. /// -/// Implements Eq. (36) of +/// Implements Eq. (36) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex { let N = c.n(); let S1 = c.get(K::S1); @@ -37,7 +37,7 @@ pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex { CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) } -/// Compute the $O(a_s^1a_{em}^1)$ quark-photon anomalous dimension. +/// Compute the quark-photon anomalous dimension. /// /// Implements Eq. (26) of pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { @@ -63,9 +63,9 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { 2.0 * (nf as f64) * CA * CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) } -/// Compute the $O(a_s^1a_{em}^1)$ gluon-photon anomalous dimension. +/// Compute the gluon-photon anomalous dimension. /// -/// Implements Eq. (27) of +/// Implements Eq. (27) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex { let N = c.n(); CF * CA @@ -73,30 +73,30 @@ pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex { / (N.powu(3) * (1.0 + N).powu(3) * (-2.0 + N + N.powu(2))) } -/// Compute the $O(a_s^1a_{em}^1)$ photon-gluon anomalous dimension. +/// Compute the photon-gluon anomalous dimension. /// -/// Implements Eq. (30) of +/// Implements Eq. (30) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_phg(c: &mut Cache, nf: u8) -> Complex { TR / CF / CA * (NC as f64) * gamma_gph(c, nf) } -/// Compute the $O(a_s^1a_{em}^1)$ quark-gluon singlet anomalous dimension. +/// Compute the quark-gluon singlet anomalous dimension. /// -/// Implements Eq. (29) of +/// Implements Eq. (29) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { TR / CF / CA * (NC as f64) * gamma_qph(c, nf) } -/// Compute the $O(a_s^1a_{em}^1)$ gluon-quark singlet anomalous dimension. +/// Compute the gluon-quark singlet anomalous dimension. /// -/// Implements Eq. (35) of +/// Implements Eq. (35) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { gamma_phq(c, nf) } -/// Compute the $O(a_s^1a_{em}^1)$ photon-photon singlet anomalous dimension. +/// Compute the photon-photon singlet anomalous dimension. /// -/// Implements Eq. (28) of +/// Implements Eq. (28) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { let cc = ChargeCombinations { nf }; cmplx!( @@ -105,16 +105,16 @@ pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { ) } -/// Compute the $O(a_s^1a_{em}^1)$ gluon-gluon singlet anomalous dimension. +/// Compute the gluon-gluon singlet anomalous dimension. /// -/// Implements Eq. (31) of +/// Implements Eq. (31) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex { cmplx!(4.0 * TR * (NC as f64), 0.) } -/// Compute the $O(a_s^1a_{em}^1)$ singlet-like non singlet anomalous dimension. +/// Compute the singlet-like non-singlet anomalous dimension. /// -/// Implements Eqs. (33-34) of +/// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { let N = c.n(); let S1 = c.get(K::S1); @@ -144,9 +144,9 @@ pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { CF * result } -/// Compute the $O(a_s^1a_{em}^1)$ valence-like non singlet anomalous dimension. +/// Compute the valence-like non-singlet anomalous dimension. /// -/// Implements Eqs. (33-34) of +/// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { let N = c.n(); let S1 = c.get(K::S1); @@ -194,7 +194,7 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { CF * result } -/// Compute the $O(a_s^1a_{em}^1)$ singlet sector. +/// Compute the singlet anomalous dimension matrix. pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = ChargeCombinations { nf }; let e2_tot = nf as f64 * cc.e2avg(); @@ -227,7 +227,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { ] } -/// Compute the $O(a_s^1a_{em}^1)$ valence sector. +/// Compute the valence anomalous dimension matrix. pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; [ diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs index 06533f6e2..b7c65a3d4 100644 --- a/crates/ekore/src/bib.rs +++ b/crates/ekore/src/bib.rs @@ -1,4 +1,4 @@ -//! List of References (autogenerated on 2025-01-13T14:46:34.886936). +//! List of References (autogenerated on 2025-01-13T14:59:31.519662). #[allow(non_snake_case)] /// The Three loop splitting functions in QCD: The Nonsinglet case @@ -131,3 +131,15 @@ pub fn Gluck1995yr() {} /// /// DOI: [10.13130/carrazza-stefano_phd2015-07-06](https:dx.doi.org/10.13130/carrazza-stefano_phd2015-07-06) pub fn Carrazza2015dea() {} + +#[allow(non_snake_case)] +/// QED corrections to the Altarelli-Parisi splitting functions +/// +/// de Florian, Daniel and Sborlini, Germán F. R. and Rodrigo, Germán +/// +/// Published in: Eur. Phys. J. C 76 (2016), 282 +/// +/// e-Print: [1512.00612](https://arxiv.org/abs/1512.00612) +/// +/// DOI: [10.1140/epjc/s10052-016-4131-8](https:dx.doi.org/10.1140/epjc/s10052-016-4131-8) +pub fn deFlorian2015ujt() {} diff --git a/crates/make_bib.py b/crates/make_bib.py index e758c41ad..ef1f06472 100644 --- a/crates/make_bib.py +++ b/crates/make_bib.py @@ -25,9 +25,15 @@ PHD = """Published as PhD thesis at {school} ({year})""" -def clean_nl(t: str) -> str: - """Shrink whitespace to one.""" - return re.sub("\n\\s+", " ", t) +def sanitize(t: str) -> str: + """Reformat for HTML.""" + # Shrink whitespace to one. + t = re.sub("\n\\s+", " ", t) + # Recover some TeX + t = t.replace(r"\textendash{}", "-").replace(r"\'a", "á") + # Remove curly brackets + t = re.sub(r"^\{(.+)\}$", r"\1", t) + return t # collect output @@ -39,8 +45,8 @@ def clean_nl(t: str) -> str: 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 + title = sanitize(el.fields_dict["title"].value) + author = sanitize(el.fields_dict["author"].value) if el.entry_type == "phdthesis": publication = PHD.format( school=el.fields_dict["school"].value, From 1dc1ffda2604ae93b327ff8cb6da3754a36ca095 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 15:38:02 +0200 Subject: [PATCH 38/71] rust: Init ad/u/s/aem2 --- crates/ekore/refs.bib | 13 ++ .../unpolarized/spacelike.rs | 1 + .../unpolarized/spacelike/aem2.rs | 158 ++++++++++++++++++ crates/ekore/src/bib.rs | 14 +- 4 files changed, 185 insertions(+), 1 deletion(-) create mode 100644 crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib index dc225226e..13d85e8b5 100644 --- a/crates/ekore/refs.bib +++ b/crates/ekore/refs.bib @@ -143,3 +143,16 @@ @article{deFlorian2015ujt pages = "282", year = "2016" } +@article{deFlorian2016gvk, + author = "de Florian, Daniel and Sborlini, Germ\'an F. R. and Rodrigo, Germ\'an", + title = "{Two-loop QED corrections to the Altarelli-Parisi splitting functions}", + eprint = "1606.02887", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "ICAS-09-16, IFIC-15-88", + doi = "10.1007/JHEP10(2016)056", + journal = "JHEP", + volume = "10", + pages = "056", + year = "2016" +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index b310289c3..f9a3ff149 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -7,6 +7,7 @@ use crate::harmonics::cache::Cache; use num::complex::Complex; use num::Zero; pub mod aem1; +pub mod aem2; pub mod as1; pub mod as1aem1; pub mod as2; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs new file mode 100644 index 000000000..0e5fdda5e --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs @@ -0,0 +1,158 @@ +//! |NLO| |QED|. +use num::complex::Complex; +use num::Zero; + +use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC}; +use crate::harmonics::cache::{Cache, K}; + +use crate::anomalous_dimensions::unpolarized::spacelike::as1aem1; + +/// Compute the photon-photon anomalous dimension. +/// +/// Implements Eq. (68) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_phph(c: &mut Cache, nf: u8) -> Complex { + let cc = ChargeCombinations { nf }; + (NC as f64) + * ((cc.nu() as f64) * EU2.powi(2) + (cc.nd() as f64) * ED2.powi(2)) + * (as1aem1::gamma_gph(c, nf) / CF / CA + 4.) +} + +/// Compute the up-photon anomalous dimension. +/// +/// Implements Eq. (55) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u. +pub fn gamma_uph(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_qph(c, nf) / CF +} + +/// Compute the down-photon anomalous dimension. +/// +/// Implements Eq. (55) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d. +pub fn gamma_dph(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_qph(c, nf) / CF +} + +/// Compute the photon-quark anomalous dimension. +/// +/// Implements singlet part of Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +fn gamma_phq(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let gamma = (-16.0 * (-16.0 - 27.0 * N - 13.0 * N.powu(2) - 8.0 * N.powu(3))) + / (9.0 * (-1.0 + N) * N * (1.0 + N).powu(2)) + - 16.0 * (2.0 + 3.0 * N + 2.0 * N.powu(2) + N.powu(3)) + / (3.0 * (-1.0 + N) * N * (1.0 + N).powu(2)) + * S1; + let cc = ChargeCombinations { nf }; + let eSigma2 = (NC as f64) * cc.e2avg() * (nf as f64); + eSigma2 * gamma +} + +/// Compute the photon-up anomalous dimension. +/// +/// Implements Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_phu(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf) +} + +/// Compute the photon-down anomalous dimension. +/// +/// Implements Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_phd(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf) +} + +// /// Compute the non-singlet anomalous dimension. +// /// +// /// Implements Eq. (2.5) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +// pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { +// as1::gamma_ns(c, nf) / CF +// } + +// /// Compute the singlet anomalous dimension matrix. +// pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { +// let cc = ChargeCombinations { nf }; + +// let gamma_ph_q = gamma_phq(c, nf); +// let gamma_q_ph = gamma_qph(c, nf); +// let gamma_nonsinglet = gamma_ns(c, nf); + +// [ +// [ +// Complex::::zero(), +// Complex::::zero(), +// Complex::::zero(), +// Complex::::zero(), +// ], +// [ +// Complex::::zero(), +// gamma_phph(c, nf), +// cc.e2avg() * gamma_ph_q, +// cc.vue2m() * gamma_ph_q, +// ], +// [ +// Complex::::zero(), +// cc.e2avg() * gamma_q_ph, +// cc.e2avg() * gamma_nonsinglet, +// cc.vue2m() * gamma_nonsinglet, +// ], +// [ +// Complex::::zero(), +// cc.vde2m() * gamma_q_ph, +// cc.vde2m() * gamma_nonsinglet, +// cc.e2delta() * gamma_nonsinglet, +// ], +// ] +// } + +// /// Compute the valence anomalous dimension matrix. +// pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { +// let cc = ChargeCombinations { nf }; + +// [ +// [cc.e2avg() * gamma_ns(c, nf), cc.vue2m() * gamma_ns(c, nf)], +// [cc.vde2m() * gamma_ns(c, nf), cc.e2delta() * gamma_ns(c, nf)], +// ] +// } + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + 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_cmplx!(f64, me, Complex::zero(), 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_phq(&mut c, NF); + // assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); + // } + + #[test] + fn photon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + assert_approx_eq_cmplx!( + f64, + EU2 * gamma_uph(&mut c, cc.nu()) + + ED2 * gamma_dph(&mut c, cc.nd()) + + gamma_phph(&mut c, nf), + Complex::zero(), + epsilon = 1e-12 + ); + } + } +} diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs index b7c65a3d4..89ca13ccf 100644 --- a/crates/ekore/src/bib.rs +++ b/crates/ekore/src/bib.rs @@ -1,4 +1,4 @@ -//! List of References (autogenerated on 2025-01-13T14:59:31.519662). +//! List of References (autogenerated on 2025-01-13T15:14:10.528903). #[allow(non_snake_case)] /// The Three loop splitting functions in QCD: The Nonsinglet case @@ -143,3 +143,15 @@ pub fn Carrazza2015dea() {} /// /// DOI: [10.1140/epjc/s10052-016-4131-8](https:dx.doi.org/10.1140/epjc/s10052-016-4131-8) pub fn deFlorian2015ujt() {} + +#[allow(non_snake_case)] +/// Two-loop QED corrections to the Altarelli-Parisi splitting functions +/// +/// de Florian, Daniel and Sborlini, Germán F. R. and Rodrigo, Germán +/// +/// Published in: JHEP 10 (2016), 056 +/// +/// e-Print: [1606.02887](https://arxiv.org/abs/1606.02887) +/// +/// DOI: [10.1007/JHEP10(2016)056](https:dx.doi.org/10.1007/JHEP10(2016)056) +pub fn deFlorian2016gvk() {} From b7b314fa6ba708d33d567f770562e21e5571d1bb Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 15:53:44 +0200 Subject: [PATCH 39/71] rust: Expand ad/u/s/aem2 --- .../unpolarized/spacelike/aem2.rs | 65 +++++++++++++++---- 1 file changed, 52 insertions(+), 13 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs index 0e5fdda5e..378cacf94 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs @@ -61,12 +61,49 @@ pub fn gamma_phd(c: &mut Cache, nf: u8) -> Complex { ED2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf) } -// /// Compute the non-singlet anomalous dimension. -// /// -// /// Implements Eq. (2.5) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. -// pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { -// as1::gamma_ns(c, nf) / CF -// } +/// Compute the non-singlet anomalous dimension. +/// +/// Implements singlet part of Eq. (57) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +fn gamma_nsq(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let gamma = 2.0 * (-12.0 + 20.0 * N + 47.0 * N.powu(2) + 6.0 * N.powu(3) + 3.0 * N.powu(4)) + / (9.0 * N.powu(2) * (1.0 + N).powu(2)) + - 80.0 / 9.0 * S1 + + 16.0 / 3.0 * S2; + let cc = ChargeCombinations { nf }; + let eSigma2 = (NC as f64) * cc.e2avg() * (nf as f64); + eSigma2 * gamma +} + +/// Compute the singlet-like non-singlet up anomalous dimension. +/// +/// Implements sum of Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u. +pub fn gamma_nspu(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_nsp(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the singlet-like non-singlet down anomalous dimension. +/// +/// Implements sum of Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d. +pub fn gamma_nspd(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_nsp(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the valence-like non-singlet up anomalous dimension. +/// +/// Implements difference between Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u. +pub fn gamma_nsmu(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the valence-like non-singlet down anomalous dimension. +/// +/// Implements difference between Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d. +pub fn gamma_nsmd(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf) +} // /// Compute the singlet anomalous dimension matrix. // pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { @@ -122,13 +159,15 @@ mod tests { use num::Zero; 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_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); - // } + #[test] + fn number_conservation() { + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + for nf in 2..7 { + assert_approx_eq_cmplx!(f64, gamma_nsmu(&mut c, nf), Complex::zero(), epsilon = 1e-5); + assert_approx_eq_cmplx!(f64, gamma_nsmd(&mut c, nf), Complex::zero(), epsilon = 1e-5); + } + } // #[test] // fn quark_momentum_conservation() { From 180cc7e5727dd07595da7973f043ee49eec92154 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 16:39:58 +0200 Subject: [PATCH 40/71] rust: Complete ad/u/s/aem2 --- .../unpolarized/spacelike/aem2.rs | 135 +++++++++++------- crates/ekore/src/constants.rs | 16 ++- 2 files changed, 95 insertions(+), 56 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs index 378cacf94..1a39ef6a6 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs @@ -105,51 +105,75 @@ pub fn gamma_nsmd(c: &mut Cache, nf: u8) -> Complex { ED2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf) } -// /// Compute the singlet anomalous dimension matrix. -// pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { -// let cc = ChargeCombinations { nf }; - -// let gamma_ph_q = gamma_phq(c, nf); -// let gamma_q_ph = gamma_qph(c, nf); -// let gamma_nonsinglet = gamma_ns(c, nf); - -// [ -// [ -// Complex::::zero(), -// Complex::::zero(), -// Complex::::zero(), -// Complex::::zero(), -// ], -// [ -// Complex::::zero(), -// gamma_phph(c, nf), -// cc.e2avg() * gamma_ph_q, -// cc.vue2m() * gamma_ph_q, -// ], -// [ -// Complex::::zero(), -// cc.e2avg() * gamma_q_ph, -// cc.e2avg() * gamma_nonsinglet, -// cc.vue2m() * gamma_nonsinglet, -// ], -// [ -// Complex::::zero(), -// cc.vde2m() * gamma_q_ph, -// cc.vde2m() * gamma_nonsinglet, -// cc.e2delta() * gamma_nonsinglet, -// ], -// ] -// } - -// /// Compute the valence anomalous dimension matrix. -// pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { -// let cc = ChargeCombinations { nf }; - -// [ -// [cc.e2avg() * gamma_ns(c, nf), cc.vue2m() * gamma_ns(c, nf)], -// [cc.vde2m() * gamma_ns(c, nf), cc.e2delta() * gamma_ns(c, nf)], -// ] -// } +/// Compute the pure-singlet anomalous dimension. +/// +/// Implements Eq. (59) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let result = -4.0 * (2.0 + N * (5.0 + N)) * (4.0 + N * (4.0 + N * (7.0 + 5.0 * N))) + / ((-1.0 + N) * N.powu(3) * (1.0 + N).powu(3) * (2.0 + N).powu(2)); + 2. * (nf as f64) * CA * result +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + let cc = ChargeCombinations { nf }; + const E2M: f64 = EU2 - ED2; + + let gamma_phu = gamma_phu(c, nf); + let gamma_phd = gamma_phd(c, nf); + let gamma_uph = gamma_uph(c, nf); + let gamma_dph = gamma_dph(c, nf); + let gamma_nspu = gamma_nspu(c, nf); + let gamma_nspd = gamma_nspd(c, nf); + let gamma_ps = gamma_ps(c, nf); + + [ + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + Complex::::zero(), + gamma_phph(c, nf), + cc.vu() * EU2 * gamma_phu + cc.vd() * ED2 * gamma_phd, + cc.vu() * (EU2 * gamma_phu - ED2 * gamma_phd), + ], + [ + Complex::::zero(), + cc.vu() * EU2 * gamma_uph + cc.vd() * ED2 * gamma_dph, + cc.vu() * EU2 * gamma_nspu + cc.vd() * ED2 * gamma_nspd + cc.e2avg().powi(2) * gamma_ps, + cc.vu() * (EU2 * gamma_nspu - ED2 * gamma_nspd + E2M * cc.e2avg() * gamma_ps), + ], + [ + Complex::::zero(), + cc.vd() * (EU2 * gamma_uph - ED2 * gamma_dph), + cc.vd() * (EU2 * gamma_nspu - ED2 * gamma_nspd + E2M * cc.e2avg() * gamma_ps), + cc.vd() * EU2 * gamma_nspu + + cc.vu() * ED2 * gamma_nspd + + cc.vu() * cc.vd() * E2M.powi(2) * gamma_ps, + ], + ] +} + +/// Compute the valence anomalous dimension matrix. +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let cc = ChargeCombinations { nf }; + let gamma_nsmu = gamma_nsmu(c, nf); + let gamma_nsmd = gamma_nsmd(c, nf); + [ + [ + cc.vu() * EU2 * gamma_nsmu + cc.vd() * ED2 * gamma_nsmd, + cc.vu() * (EU2 * gamma_nsmu - ED2 * gamma_nsmd), + ], + [ + cc.vd() * (EU2 * gamma_nsmu - ED2 * gamma_nsmd), + cc.vd() * EU2 * gamma_nsmu + cc.vu() * ED2 * gamma_nsmd, + ], + ] +} #[cfg(test)] mod tests { @@ -157,7 +181,6 @@ mod tests { use crate::{assert_approx_eq_cmplx, cmplx}; use num::complex::Complex; use num::Zero; - const NF: u8 = 5; #[test] fn number_conservation() { @@ -169,13 +192,19 @@ mod tests { } } - // #[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_phq(&mut c, NF); - // assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); - // } + #[test] + fn quark_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + let ps = EU2 * gamma_ps(&mut c, cc.nu()) + ED2 * gamma_ps(&mut c, cc.nd()); + let me_u = gamma_nspu(&mut c, nf) + ps + gamma_phu(&mut c, nf); + assert_approx_eq_cmplx!(f64, me_u, Complex::zero(), epsilon = 1e-5); + let me_d = gamma_nspd(&mut c, nf) + ps + gamma_phd(&mut c, nf); + assert_approx_eq_cmplx!(f64, me_d, Complex::zero(), epsilon = 1e-5); + } + } #[test] fn photon_momentum_conservation() { diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 10c8cc2d6..b860b96c5 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -80,19 +80,29 @@ impl ChargeCombinations { self.nf - self.nu() } - /// Electric charge average. + /// Charge average. pub fn e2avg(&self) -> f64 { (self.nu() as f64 * EU2 + self.nd() as f64 * ED2) / (self.nf as f64) } + /// Relative up contribution. + pub fn vu(&self) -> f64 { + self.nu() as f64 / (self.nf as f64) + } + + /// Relative down contribution. + pub fn vd(&self) -> f64 { + self.nd() as f64 / (self.nf as f64) + } + /// Relative up contribution to charge difference. pub fn vue2m(&self) -> f64 { - self.nu() as f64 / (self.nf as f64) * (EU2 - ED2) + self.vu() * (EU2 - ED2) } /// Relative down contribution to charge difference. pub fn vde2m(&self) -> f64 { - self.nd() as f64 / (self.nf as f64) * (EU2 - ED2) + self.vd() * (EU2 - ED2) } /// Asymmetric charge combination. From 7bed428b6461368604c0c2881b3b06b9573515d8 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 13 Jan 2025 17:46:37 +0200 Subject: [PATCH 41/71] rust: Start generalizing ad/u/spacelike --- .../unpolarized/spacelike.rs | 67 +++++++++++-------- crates/ekore/src/constants.rs | 8 +-- 2 files changed, 44 insertions(+), 31 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index f9a3ff149..a4313755d 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -1,7 +1,7 @@ //! The unpolarized, space-like anomalous dimensions at various couplings power. use crate::constants::{ - ED2, EU2, PID_NSM, PID_NSM_ED2, PID_NSM_EU2, PID_NSP, PID_NSP_ED2, PID_NSP_EU2, PID_NSV, + ED2, EU2, PID_NSM, PID_NSM_D, PID_NSM_U, PID_NSP, PID_NSP_D, PID_NSP_U, PID_NSV, }; use crate::harmonics::cache::Cache; use num::complex::Complex; @@ -61,7 +61,7 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl gamma_S } -/// Compute the tower of the QED non-singlet anomalous dimensions. +/// Compute the tower of the |QCD| x |QED| non-singlet anomalous dimensions. pub fn gamma_ns_qed( order_qcd: usize, order_qed: usize, @@ -71,30 +71,44 @@ pub fn gamma_ns_qed( ) -> Vec>> { let col = vec![Complex::::zero(); order_qcd + 1]; let mut gamma_ns = vec![col; order_qed + 1]; - gamma_ns[1][0] = as1::gamma_ns(c, nf); - gamma_ns[0][1] = choose_ns_as_aem1(mode, c, nf); - gamma_ns[1][1] = choose_ns_as_as1aem1(mode, c, nf); - gamma_ns -} -pub fn choose_ns_as_aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { - match mode { - PID_NSP_EU2 | PID_NSM_EU2 => EU2 * aem1::gamma_ns(c, nf), - PID_NSP_ED2 | PID_NSM_ED2 => ED2 * aem1::gamma_ns(c, nf), - _ => panic!("Unkown non-singlet sector element"), + // QCD corrections + let qcd_mode = match mode { + PID_NSP_U | PID_NSP_D => PID_NSP, + PID_NSM_U | PID_NSM_D => PID_NSM, + _ => mode, + }; + let gamma_qcd = gamma_ns_qcd(order_qcd, qcd_mode, c, nf); + for j in 0..order_qcd { + gamma_ns[1 + j][0] = gamma_qcd[j]; } -} - -pub fn choose_ns_as_as1aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex { - match mode { - PID_NSP_EU2 => EU2 * as1aem1::gamma_nsp(c, nf), - PID_NSP_ED2 => ED2 * as1aem1::gamma_nsp(c, nf), - PID_NSM_EU2 => EU2 * as1aem1::gamma_nsm(c, nf), - PID_NSM_ED2 => ED2 * as1aem1::gamma_nsm(c, nf), + // QED corrections + gamma_ns[0][1] = match mode { + PID_NSP_U | PID_NSM_U => EU2 * aem1::gamma_ns(c, nf), + PID_NSP_D | PID_NSM_D => ED2 * aem1::gamma_ns(c, nf), _ => panic!("Unkown non-singlet sector element"), + }; + if order_qed >= 2 { + gamma_ns[0][2] = match mode { + PID_NSP_U => EU2 * aem2::gamma_nspu(c, nf), + PID_NSP_D => ED2 * aem2::gamma_nspd(c, nf), + PID_NSM_U => EU2 * aem2::gamma_nsmu(c, nf), + PID_NSM_D => ED2 * aem2::gamma_nsmd(c, nf), + _ => panic!("Unkown non-singlet sector element"), + }; } + // QCDxQED corrections + gamma_ns[1][1] = match mode { + PID_NSP_U => EU2 * as1aem1::gamma_nsp(c, nf), + PID_NSP_D => ED2 * as1aem1::gamma_nsp(c, nf), + PID_NSM_U => EU2 * as1aem1::gamma_nsm(c, nf), + PID_NSM_D => ED2 * as1aem1::gamma_nsm(c, nf), + _ => panic!("Unkown non-singlet sector element"), + }; + gamma_ns } +/// Compute the tower of the |QCD| x |QED| singlet anomalous dimensions matrices. pub fn gamma_singlet_qed( order_qcd: usize, order_qed: usize, @@ -110,7 +124,6 @@ pub fn gamma_singlet_qed( ]; 4]; order_qcd + 1 ]; - let mut gamma_s = vec![col; order_qed + 1]; gamma_s[1][0] = as1::gamma_singlet_qed(c, nf); @@ -119,7 +132,7 @@ pub fn gamma_singlet_qed( gamma_s } -/// Compute the grid of the QED valence anomalous dimensions matrices +/// Compute the tower of the |QCD| x |QED| valence anomalous dimensions matrices. pub fn gamma_valence_qed( order_qcd: usize, order_qed: usize, @@ -127,8 +140,8 @@ pub fn gamma_valence_qed( nf: u8, ) -> Vec; 2]; 2]>> { let col = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; - let mut gamma_v = vec![col; order_qed + 1]; + gamma_v[1][0] = as1::gamma_valence_qed(c, nf); gamma_v[0][1] = aem1::gamma_valence(c, nf); gamma_v[1][1] = as1aem1::gamma_valence(c, nf); @@ -191,7 +204,7 @@ mod tests { for j in [0, 1] { assert_approx_eq_cmplx!( f64, - gamma_ns_qed(1, 1, PID_NSM_EU2, &mut c, NF)[i][j], + gamma_ns_qed(1, 1, PID_NSM_U, &mut c, NF)[i][j], cmplx!(0., 0.), epsilon = 1e-5 ); @@ -202,7 +215,7 @@ mod tests { for j in [0, 1] { assert_approx_eq_cmplx!( f64, - gamma_ns_qed(1, 1, PID_NSM_ED2, &mut c, NF)[i][j], + gamma_ns_qed(1, 1, PID_NSM_D, &mut c, NF)[i][j], cmplx!(0., 0.), epsilon = 1e-5 ); @@ -211,14 +224,14 @@ mod tests { assert_approx_eq_cmplx!( f64, - gamma_ns_qed(1, 1, PID_NSP_EU2, &mut c, NF)[0][1], + gamma_ns_qed(1, 1, PID_NSP_U, &mut c, NF)[0][1], cmplx!(0., 0.), epsilon = 1e-5 ); assert_approx_eq_cmplx!( f64, - gamma_ns_qed(1, 1, PID_NSP_ED2, &mut c, NF)[0][1], + gamma_ns_qed(1, 1, PID_NSP_D, &mut c, NF)[0][1], cmplx!(0., 0.), epsilon = 1e-5 ); diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index b860b96c5..98023e016 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -53,16 +53,16 @@ pub const PID_NSM: u16 = 10201; pub const PID_NSV: u16 = 10200; /// singlet-like non-singlet up-sector |PID|. -pub const PID_NSP_EU2: u16 = 10102; +pub const PID_NSP_U: u16 = 10102; /// singlet-like non-singlet down-sector |PID|. -pub const PID_NSP_ED2: u16 = 10103; +pub const PID_NSP_D: u16 = 10103; /// valence-like non-singlet up-sector |PID|. -pub const PID_NSM_EU2: u16 = 10202; +pub const PID_NSM_U: u16 = 10202; /// valence-like non-singlet down-sector |PID|. -pub const PID_NSM_ED2: u16 = 10203; +pub const PID_NSM_D: u16 = 10203; /// |QED| electric charge combinations. pub struct ChargeCombinations { From 680dae8878134656bb488ff796365b24ea734c7e Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 14 Jan 2025 15:47:38 +0200 Subject: [PATCH 42/71] rust: Use QCD ads also in QED --- .../unpolarized/spacelike.rs | 53 ++++++++++++++++++- .../unpolarized/spacelike/as1.rs | 41 -------------- 2 files changed, 51 insertions(+), 43 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index a4313755d..2bcf99d52 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -126,8 +126,44 @@ pub fn gamma_singlet_qed( ]; let mut gamma_s = vec![col; order_qed + 1]; - gamma_s[1][0] = as1::gamma_singlet_qed(c, nf); + // QCD corrections + let gamma_qcd_s = gamma_singlet_qcd(order_qcd, c, nf); + let gamma_qcd_nsp = gamma_ns_qcd(order_qcd, PID_NSP, c, nf); + for j in 0..order_qcd { + let s = gamma_qcd_s[j]; + gamma_s[1 + j][0] = [ + [ + s[1][1], + Complex::::zero(), + s[1][0], + Complex::::zero(), + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + s[0][1], + Complex::::zero(), + s[0][0], + Complex::::zero(), + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + gamma_qcd_nsp[j], + ], + ]; + } + // QED corrections gamma_s[0][1] = aem1::gamma_singlet(c, nf); + if order_qed >= 2 { + gamma_s[0][2] = aem2::gamma_singlet(c, nf); + } + // QCDxQED corrections gamma_s[1][1] = as1aem1::gamma_singlet(c, nf); gamma_s } @@ -142,8 +178,21 @@ pub fn gamma_valence_qed( let col = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; let mut gamma_v = vec![col; order_qed + 1]; - gamma_v[1][0] = as1::gamma_valence_qed(c, nf); + // QCD corrections + let gamma_qcd_nsv = gamma_ns_qcd(order_qcd, PID_NSV, c, nf); + let gamma_qcd_nsm = gamma_ns_qcd(order_qcd, PID_NSM, c, nf); + for j in 0..order_qcd { + gamma_v[1 + j][0] = [ + [gamma_qcd_nsv[j], Complex::::zero()], + [Complex::::zero(), gamma_qcd_nsm[j]], + ]; + } + // QED corrections gamma_v[0][1] = aem1::gamma_valence(c, nf); + if order_qed >= 2 { + gamma_v[0][2] = aem2::gamma_valence(c, nf); + } + // QCDxQED corrections gamma_v[1][1] = as1aem1::gamma_valence(c, nf); gamma_v } diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs index f0c41c3a7..6d25e1c51 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -1,7 +1,6 @@ //! |LO| |QCD|. use num::complex::Complex; -use num::Zero; use crate::constants::{CA, CF, TR}; use crate::harmonics::cache::{Cache, K}; @@ -53,46 +52,6 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { ] } -/// Compute the singlet anomalous dimension matrix -/// for the unified evolution basis. -pub fn gamma_singlet_qed(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { - [ - [ - gamma_gg(c, nf), - Complex::::zero(), - gamma_gq(c, nf), - Complex::::zero(), - ], - [ - Complex::::zero(), - Complex::::zero(), - Complex::::zero(), - Complex::::zero(), - ], - [ - gamma_qg(c, nf), - Complex::::zero(), - gamma_ns(c, nf), - Complex::::zero(), - ], - [ - Complex::::zero(), - Complex::::zero(), - Complex::::zero(), - gamma_ns(c, nf), - ], - ] -} - -/// Compute the valence anomalous dimension matrix -/// for the unified evolution basis. -pub fn gamma_valence_qed(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { - [ - [gamma_ns(c, nf), Complex::::zero()], - [Complex::::zero(), gamma_ns(c, nf)], - ] -} - #[cfg(test)] mod tests { use super::*; From 5046f5015c925c23e1533f92b9ca52e12f8f444b Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 14 Jan 2025 16:19:07 +0200 Subject: [PATCH 43/71] rust: Iterate eko/lib.rs --- crates/eko/src/lib.rs | 31 +++++++++++--------- src/eko/evolution_operator/__init__.py.patch | 4 +-- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index ea42fb1e9..6f042ff86 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -7,14 +7,15 @@ use std::ffi::c_void; pub mod bib; pub mod mellin; -/// Wrapper to pass arguments back to Python +/// Wrapper to pass arguments back to Python. struct RawCmplx { re: Vec, im: Vec, } -/// Map tensors to c-ordered list -/// (res is a vector with dim order_qcd filled with DIMxDIM matrices) +/// Map tensor with shape (o,d,d) to c-ordered list. +/// +/// This is needed for the QCD singlet. fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: usize) -> RawCmplx { let mut target = RawCmplx { re: Vec::::new(), @@ -31,8 +32,9 @@ fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: us target } -/// Map tensors to c-ordered list in the QED singlet and valence case -/// (res is a matrix with dim order_qcd x order_qed filled with DIMxDIM matrices) +/// Map tensor with shape (o,o',d,d) to c-ordered list. +/// +/// This is needed for the QED singlet and valence. fn unravel_qed( res: Vec; DIM]; DIM]>>, order_qcd: usize, @@ -55,8 +57,9 @@ fn unravel_qed( target } -/// Map tensors to c-ordered list in the QED non-singlet case -/// (res is a matrix with dim order_qcd x order_qed filled with complex numbers) +/// Map tensor with shape (o,o',d) to c-ordered list. +/// +/// This is needed for the QED non-singlet. fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usize) -> RawCmplx { let mut target = RawCmplx { re: Vec::::new(), @@ -71,13 +74,13 @@ fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usiz target } -/// QCD intergration kernel inside quad. +/// Intergration kernel inside quad. /// /// # Safety /// This is the connection from Python, so we don't know what is on the other side. #[no_mangle] -pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { - let args = *(rargs as *mut QuadQCDargs); +pub unsafe extern "C" fn rust_quad_ker(u: f64, rargs: *mut c_void) -> f64 { + let args = *(rargs as *mut QuadArgs); let is_singlet = (100 == args.mode0) || (21 == args.mode0) @@ -248,7 +251,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn( #[allow(non_snake_case)] #[repr(C)] #[derive(Clone, Copy)] -pub struct QuadQCDargs { +pub struct QuadArgs { pub order_qcd: usize, pub order_qed: usize, pub mode0: u16, @@ -329,13 +332,13 @@ pub unsafe extern "C" fn my_py( /// 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 `rust_quad_ker_qcd`). +/// package (since it does not appear in the signature of `rust_quad_ker`). /// /// # Safety /// This is the connection from and back to Python, so we don't know what is on the other side. #[no_mangle] -pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { - QuadQCDargs { +pub unsafe extern "C" fn empty_args() -> QuadArgs { + QuadArgs { order_qcd: 0, order_qed: 0, mode0: 0, diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index 41ddb9f9d..d97f5c768 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -618,7 +618,7 @@ index bd1b19d6..de87651c 100644 + labels = self.labels start_time = time.perf_counter() + # start preparing C arguments -+ cfg = ekors.lib.empty_qcd_args() ++ cfg = ekors.lib.empty_args() + cfg.order_qcd = self.order[0] + cfg.order_qed = self.order[1] + cfg.is_polarized = self.config["polarized"] @@ -680,7 +680,7 @@ index bd1b19d6..de87651c 100644 + cfg.mode1 = label[1] + # construct the low level object + func = LowLevelCallable( -+ ekors.lib.rust_quad_ker_qcd, ekors.ffi.addressof(cfg) ++ ekors.lib.rust_quad_ker, ekors.ffi.addressof(cfg) + ) res = integrate.quad( - self.quad_ker( From d56b122c58dbac6b2141dd54d640625451c8be86 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 14:07:43 +0200 Subject: [PATCH 44/71] Fix LO FFNS apfel bench --- benchmarks/apfel_bench.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index b9194d68c..8558e1435 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -37,11 +37,7 @@ class BenchmarkVFNS(ApfelBenchmark): vfns_theory = { "FNS": "ZM-VFNS", - "ModEv": [ - "EXA", - "EXP", - "TRN", - ], + "ModEv": "EXA", "kcThr": 1.0, "kbThr": 1.0, "ktThr": 1.0, @@ -132,12 +128,10 @@ class BenchmarkFFNS(ApfelBenchmark): ffns_theory = { "FNS": "FFNS", - "ModEv": [ - "EXA", - "EXP", - "TRN", - ], + "ModEv": "EXA", "NfFF": 4, + "nf0": 4, + "nfref": 4, "kcThr": 0.0, "kbThr": np.inf, "ktThr": np.inf, @@ -316,14 +310,14 @@ def benchmark_plain(self, pto, qed): if __name__ == "__main__": # obj = BenchmarkVFNS() - # obj = BenchmarkFFNS() + obj = BenchmarkFFNS() # obj.benchmark_plain_pol(2) - # obj.benchmark_plain(2) + obj.benchmark_plain(0) # obj.benchmark_sv(2, "exponentiated") # obj.benchmark_kthr(2) # obj.benchmark_msbar(2) # obj = BenchmarkFFNS_qed() - obj = BenchmarkFFNS_qed() - obj.benchmark_plain(2, 2) + # obj = BenchmarkFFNS_qed() + # obj.benchmark_plain(2, 2) From f3f027ae884925dc9545e032fcc147ffcef940fd Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 15:24:21 +0200 Subject: [PATCH 45/71] rust: Fix range error --- .../ekore/src/anomalous_dimensions/unpolarized/spacelike.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 2bcf99d52..08d8aeedb 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -79,7 +79,7 @@ pub fn gamma_ns_qed( _ => mode, }; let gamma_qcd = gamma_ns_qcd(order_qcd, qcd_mode, c, nf); - for j in 0..order_qcd { + for j in 0..order_qcd - 1 { gamma_ns[1 + j][0] = gamma_qcd[j]; } // QED corrections @@ -129,7 +129,7 @@ pub fn gamma_singlet_qed( // QCD corrections let gamma_qcd_s = gamma_singlet_qcd(order_qcd, c, nf); let gamma_qcd_nsp = gamma_ns_qcd(order_qcd, PID_NSP, c, nf); - for j in 0..order_qcd { + for j in 0..order_qcd - 1 { let s = gamma_qcd_s[j]; gamma_s[1 + j][0] = [ [ @@ -181,7 +181,7 @@ pub fn gamma_valence_qed( // QCD corrections let gamma_qcd_nsv = gamma_ns_qcd(order_qcd, PID_NSV, c, nf); let gamma_qcd_nsm = gamma_ns_qcd(order_qcd, PID_NSM, c, nf); - for j in 0..order_qcd { + for j in 0..order_qcd - 1 { gamma_v[1 + j][0] = [ [gamma_qcd_nsv[j], Complex::::zero()], [Complex::::zero(), gamma_qcd_nsm[j]], From 6fea38acca09e28567bdd7f06c271a50a92cafda Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 15:24:47 +0200 Subject: [PATCH 46/71] Fix QED FFNS apfel bench --- benchmarks/apfel_bench.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index 8558e1435..2ef8f21a0 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -153,7 +153,6 @@ def benchmark_plain_pol(self, pto): th = self.ffns_theory.copy() th.update({"PTO": [pto]}) - th["ModEv"] = ["EXA"] # TODO for the time one is sufficient op = operators.apfel_config.copy() op["polarized"] = [True] self.run(cartesian_product(th), operators.build(op), ["ToyLH"]) @@ -203,6 +202,8 @@ class BenchmarkFFNS_qed(ApfelBenchmark): # "TRN", ], "NfFF": 5, + "nf0": 5, + "nfref": 5, "kcThr": 0.0, "kbThr": 0.0, "ktThr": np.inf, @@ -310,14 +311,14 @@ def benchmark_plain(self, pto, qed): if __name__ == "__main__": # obj = BenchmarkVFNS() - obj = BenchmarkFFNS() + # obj = BenchmarkFFNS() # obj.benchmark_plain_pol(2) - obj.benchmark_plain(0) + # obj.benchmark_plain(0) # obj.benchmark_sv(2, "exponentiated") # obj.benchmark_kthr(2) # obj.benchmark_msbar(2) # obj = BenchmarkFFNS_qed() - # obj = BenchmarkFFNS_qed() - # obj.benchmark_plain(2, 2) + obj = BenchmarkFFNS_qed() + obj.benchmark_plain(2, 2) From fa88d4fc886d85458ca50613565c74f0dfcc6a12 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 15:44:48 +0200 Subject: [PATCH 47/71] rust: Swap dimensions for ad --- .../anomalous_dimensions/unpolarized/spacelike.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 08d8aeedb..cada98f15 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -69,8 +69,8 @@ pub fn gamma_ns_qed( c: &mut Cache, nf: u8, ) -> Vec>> { - let col = vec![Complex::::zero(); order_qcd + 1]; - let mut gamma_ns = vec![col; order_qed + 1]; + let col = vec![Complex::::zero(); order_qed + 1]; + let mut gamma_ns = vec![col; order_qcd + 1]; // QCD corrections let qcd_mode = match mode { @@ -122,9 +122,9 @@ pub fn gamma_singlet_qed( Complex::::zero(), Complex::::zero() ]; 4]; - order_qcd + 1 + order_qed + 1 ]; - let mut gamma_s = vec![col; order_qed + 1]; + let mut gamma_s = vec![col; order_qcd + 1]; // QCD corrections let gamma_qcd_s = gamma_singlet_qcd(order_qcd, c, nf); @@ -175,8 +175,8 @@ pub fn gamma_valence_qed( c: &mut Cache, nf: u8, ) -> Vec; 2]; 2]>> { - let col = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qcd + 1]; - let mut gamma_v = vec![col; order_qed + 1]; + let col = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qed + 1]; + let mut gamma_v = vec![col; order_qcd + 1]; // QCD corrections let gamma_qcd_nsv = gamma_ns_qcd(order_qcd, PID_NSV, c, nf); From d53ca75376dcd1e4f12d4f831dc2bf50b91d1c0d Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 16:53:40 +0200 Subject: [PATCH 48/71] Use Lsv for SV in QED --- src/eko/evolution_operator/quad_ker.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index aaedb7374..731bbaa1e 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -257,7 +257,7 @@ def cb_quad_ker_ome( areas_raw, areas_x, areas_y, - L, + _L, backward_method, as1, _as0, @@ -392,7 +392,7 @@ def cb_quad_ker_qed( areas_raw, areas_x, areas_y, - L, + _L, ev_method, as1, as0, @@ -433,7 +433,7 @@ def cb_quad_ker_qed( # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_singlet = sv.exponentiated.gamma_variation_qed( - gamma_singlet, order, nf, lepton_number(mu2_to), L, alphaem_running + gamma_singlet, order, nf, lepton_number(mu2_to), Lsv, alphaem_running ) ker = qed_s.dispatcher( @@ -455,7 +455,7 @@ def cb_quad_ker_qed( alphaem_running, order, nf, - L, + Lsv, ) ) @ np.ascontiguousarray(ker) ker = select_QEDsinglet_element(ker, mode0, mode1) @@ -468,7 +468,7 @@ def cb_quad_ker_qed( if sv_mode == sv.Modes.exponentiated: gamma_valence = sv.exponentiated.gamma_variation_qed( - gamma_valence, order, nf, lepton_number(mu2_to), L, alphaem_running + gamma_valence, order, nf, lepton_number(mu2_to), Lsv, alphaem_running ) ker = qed_v.dispatcher( order, @@ -490,7 +490,7 @@ def cb_quad_ker_qed( alphaem_running, order, nf, - L, + Lsv, ) ) @ np.ascontiguousarray(ker) ker = select_QEDvalence_element(ker, mode0, mode1) @@ -502,7 +502,7 @@ def cb_quad_ker_qed( gamma_ns = re_gamma_ns + im_gamma_ns * 1j if sv_mode == sv.Modes.exponentiated: gamma_ns = sv_exponentiated.gamma_variation_qed( - gamma_ns, order, nf, lepton_number(mu2_to), L, alphaem_running + gamma_ns, order, nf, lepton_number(mu2_to), Lsv, alphaem_running ) # construct eko ker = qed_ns.dispatcher( @@ -520,7 +520,13 @@ def cb_quad_ker_qed( 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 + gamma_ns, + as_list[-1], + a_half[-1][1], + alphaem_running, + order, + nf, + Lsv, ) * ker ) From 8dc34521d2781fc5b7996c94a2d882c80d77bf07 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 17:14:52 +0200 Subject: [PATCH 49/71] Don't pass L back from Rust to Python --- crates/eko/src/lib.rs | 73 ++++++++++++-------------- src/eko/evolution_operator/quad_ker.py | 4 -- 2 files changed, 35 insertions(+), 42 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 6f042ff86..13ecf7913 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -188,7 +188,6 @@ pub unsafe extern "C" fn rust_quad_ker(u: f64, rargs: *mut c_void) -> f64 { args.areas, args.areas_x, args.areas_y, - args.L, args.method_num, args.as1, args.as0, @@ -210,41 +209,40 @@ pub unsafe extern "C" fn rust_quad_ker(u: f64, rargs: *mut c_void) -> f64 { } /// Python callback signature -type PyQuadKerQCDT = unsafe extern "C" fn( - *const f64, - *const f64, - f64, - f64, - f64, - f64, - usize, - usize, - bool, - u16, - u16, - u8, - bool, - f64, - *const f64, - u8, - u8, - f64, - u8, - f64, - f64, - u8, - u8, - u8, - bool, - f64, - *const f64, - u8, - f64, - f64, - *const f64, - u8, - u8, - bool, +type PyQuadKerT = unsafe extern "C" fn( + *const f64, // re_gamma + *const f64, // im_gamma + f64, // re_n + f64, // im_n + f64, // re_jac + f64, // im_jac + usize, // order_qcd + usize, // order_qed + bool, // is_singlet + u16, // mode0 + u16, // mode1 + u8, // nf + bool, // is_log + f64, // logx + *const f64, // areas + u8, // areas_x + u8, // areas_y + u8, // method_num + f64, // as1 + f64, // as0 + u8, // ev_op_iterations + u8, // ev_op_max_order_qcd + u8, // sv_mode_num + bool, // is_threshold + f64, // lsv + *const f64, // as_list + u8, // as_list_len + f64, // mu2_from + f64, // mu2_to + *const f64, // a_half + u8, // a_half_x + u8, // a_half_y + bool, // alphaem_running ) -> f64; /// Additional integration parameters @@ -259,7 +257,7 @@ pub struct QuadArgs { pub is_polarized: bool, pub is_time_like: bool, pub nf: u8, - pub py: PyQuadKerQCDT, + pub py: PyQuadKerT, pub is_log: bool, pub logx: f64, pub areas: *const f64, @@ -308,7 +306,6 @@ pub unsafe extern "C" fn my_py( _areas: *const f64, _areas_x: u8, _areas_y: u8, - _l: f64, _method_num: u8, _as1: f64, _as0: f64, diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index 731bbaa1e..f45d82e26 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -61,7 +61,6 @@ def select_singlet_element(ker, mode0, mode1): nb.types.CPointer(nb.types.double), # areas_raw nb.types.uintc, # areas_x nb.types.uintc, # areas_y - nb.types.double, # L nb.types.uintc, # method_num nb.types.double, # as1 nb.types.double, # as0 @@ -104,7 +103,6 @@ def cb_quad_ker_qcd( areas_raw, areas_x, areas_y, - _L, ev_method, as1, as0, @@ -257,7 +255,6 @@ def cb_quad_ker_ome( areas_raw, areas_x, areas_y, - _L, backward_method, as1, _as0, @@ -392,7 +389,6 @@ def cb_quad_ker_qed( areas_raw, areas_x, areas_y, - _L, ev_method, as1, as0, From 714efa9a57eb50361aa06ae0ac8c0b374497e0aa Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 15 Jan 2025 17:35:32 +0200 Subject: [PATCH 50/71] rust: Fix typo in aem1/v --- .../anomalous_dimensions/unpolarized/spacelike/aem1.rs | 6 +++--- .../anomalous_dimensions/unpolarized/spacelike/as1aem1.rs | 8 +++----- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index 872674927..fd0163a21 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -77,10 +77,10 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { /// Compute the valence anomalous dimension matrix. pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; - + let g = gamma_ns(c, nf); [ - [cc.e2avg() * gamma_ns(c, nf), cc.vue2m() * gamma_ns(c, nf)], - [cc.vde2m() * gamma_ns(c, nf), cc.e2delta() * gamma_ns(c, nf)], + [cc.e2avg() * g, cc.vue2m() * g], + [cc.vde2m() * g, cc.e2delta() * g], ] } diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index e8e9cc4bb..da2ccae35 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -230,12 +230,10 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { /// Compute the valence anomalous dimension matrix. pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { let cc = ChargeCombinations { nf }; + let g = gamma_nsm(c, nf); [ - [cc.e2avg() * gamma_nsm(c, nf), cc.vue2m() * gamma_nsm(c, nf)], - [ - cc.vde2m() * gamma_nsm(c, nf) * gamma_nsm(c, nf), - cc.e2delta() * gamma_nsm(c, nf), - ], + [cc.e2avg() * g, cc.vue2m() * g], + [cc.vde2m() * g, cc.e2delta() * g], ] } From b598f4d83f0bc6c69d06ccef6acae4616e41fb43 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 16 Jan 2025 16:45:44 +0200 Subject: [PATCH 51/71] Fix size of transferred arrays --- crates/eko/src/lib.rs | 8 ++++---- src/eko/evolution_operator/quad_ker.py | 12 ++++++------ 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 13ecf7913..b79cc48b2 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -44,8 +44,8 @@ fn unravel_qed( re: Vec::::new(), im: Vec::::new(), }; - for obj_ in res.iter().take(order_qcd) { - for obj in obj_.iter().take(order_qed) { + for obj_ in res.iter().take(order_qcd + 1) { + for obj in obj_.iter().take(order_qed + 1) { for col in obj.iter().take(DIM) { for el in col.iter().take(DIM) { target.re.push(el.re); @@ -65,8 +65,8 @@ fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usiz re: Vec::::new(), im: Vec::::new(), }; - for col in res.iter().take(order_qcd) { - for el in col.iter().take(order_qed) { + for col in res.iter().take(order_qcd + 1) { + for el in col.iter().take(order_qed + 1) { target.re.push(el.re); target.im.push(el.im); } diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index f45d82e26..50f7b8100 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -422,8 +422,8 @@ def cb_quad_ker_qed( if is_singlet: # reconstruct singlet matrices - re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, order_qed, 4, 4)) - im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, order_qed, 4, 4)) + re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd + 1, order_qed + 1, 4, 4)) + im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd + 1, order_qed + 1, 4, 4)) gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j # scale var exponentiated is directly applied on gamma @@ -458,8 +458,8 @@ def cb_quad_ker_qed( elif is_valence: # reconstruct valence matrices - re_gamma_valence = nb.carray(re_gamma_raw, (order_qcd, order_qed, 2, 2)) - im_gamma_valence = nb.carray(im_gamma_raw, (order_qcd, order_qed, 2, 2)) + re_gamma_valence = nb.carray(re_gamma_raw, (order_qcd + 1, order_qed + 1, 2, 2)) + im_gamma_valence = nb.carray(im_gamma_raw, (order_qcd + 1, order_qed + 1, 2, 2)) gamma_valence = re_gamma_valence + im_gamma_valence * 1j if sv_mode == sv.Modes.exponentiated: @@ -493,8 +493,8 @@ def cb_quad_ker_qed( else: # construct non-singlet matrices - re_gamma_ns = nb.carray(re_gamma_raw, (order_qcd, order_qed)) - im_gamma_ns = nb.carray(im_gamma_raw, (order_qcd, order_qed)) + re_gamma_ns = nb.carray(re_gamma_raw, (order_qcd + 1, order_qed + 1)) + im_gamma_ns = nb.carray(im_gamma_raw, (order_qcd + 1, order_qed + 1)) gamma_ns = re_gamma_ns + im_gamma_ns * 1j if sv_mode == sv.Modes.exponentiated: gamma_ns = sv_exponentiated.gamma_variation_qed( From 4633ca91f9d9ead238ed070a2da4c0adbd5749a7 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 17 Jan 2025 11:37:07 +0200 Subject: [PATCH 52/71] rust: Add minor fixes to ad/u/s/{as1}aem1 --- .../unpolarized/spacelike/aem1.rs | 6 +++- .../unpolarized/spacelike/as1aem1.rs | 32 +++++++++++-------- 2 files changed, 23 insertions(+), 15 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs index fd0163a21..2cbb04b24 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -2,6 +2,7 @@ use num::complex::Complex; use num::Zero; +use crate::cmplx; use crate::constants::{ChargeCombinations, CF, ED2, EU2, NC, TR}; use crate::harmonics::cache::Cache; @@ -28,7 +29,10 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { let cc = ChargeCombinations { nf }; - (4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2)).into() + cmplx!( + 4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2), + 0. + ) } /// Compute the non-singlet anomalous dimension. diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index da2ccae35..99e469151 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -39,7 +39,7 @@ pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex { /// Compute the quark-photon anomalous dimension. /// -/// Implements Eq. (26) of +/// Implements Eq. (26) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { let N = c.n(); let S1 = c.get(K::S1); @@ -198,31 +198,35 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { let cc = ChargeCombinations { nf }; let e2_tot = nf as f64 * cc.e2avg(); - + let gq = gamma_gq(c, nf); + let phq = gamma_phq(c, nf); + let qg = gamma_qg(c, nf); + let qph = gamma_qph(c, nf); + let nsp = gamma_nsp(c, nf); [ [ e2_tot * gamma_gg(c, nf), e2_tot * gamma_gph(c, nf), - cc.e2avg() * gamma_gq(c, nf), - cc.vue2m() * gamma_gq(c, nf), + cc.e2avg() * gq, + cc.vue2m() * gq, ], [ e2_tot * gamma_phg(c, nf), gamma_phph(c, nf), - cc.e2avg() * gamma_phq(c, nf), - cc.vue2m() * gamma_phq(c, nf), + cc.e2avg() * phq, + cc.vue2m() * phq, ], [ - cc.e2avg() * gamma_qg(c, nf), - cc.e2avg() * gamma_qph(c, nf), - cc.e2avg() * gamma_nsp(c, nf), - cc.vue2m() * gamma_nsp(c, nf), + cc.e2avg() * qg, + cc.e2avg() * qph, + cc.e2avg() * nsp, + cc.vue2m() * nsp, ], [ - cc.vde2m() * gamma_qg(c, nf), - cc.vde2m() * gamma_qph(c, nf), - cc.vde2m() * gamma_nsp(c, nf), - cc.e2delta() * gamma_nsp(c, nf), + cc.vde2m() * qg, + cc.vde2m() * qph, + cc.vde2m() * nsp, + cc.e2delta() * nsp, ], ] } From ec743caee77af00f04625e0ce35cd191d82c9bd9 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 17 Jan 2025 12:33:30 +0200 Subject: [PATCH 53/71] py: Make repetition in ad/u/s/aem2 more explicit --- src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py index db0f01c38..54be58b16 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py @@ -254,7 +254,7 @@ def gamma_nsmu(N, nf, cache): nd = nf - nu eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) tmp = ( - 2 + 2.0 * (-12.0 + 20.0 * N + 47.0 * N**2 + 6.0 * N**3 + 3.0 * N**4) / (9.0 * N**2 * (1.0 + N) ** 2) - 80.0 / 9.0 * S1 From 1a77d689981c54f54550cc206a18fc9a71346b18 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 22 Jan 2025 18:01:13 +0200 Subject: [PATCH 54/71] Add void external --- benchmarks/sandbox.py | 18 +++++++++--------- src/ekomark/benchmark/runner.py | 15 +++++++++++++++ 2 files changed, 24 insertions(+), 9 deletions(-) diff --git a/benchmarks/sandbox.py b/benchmarks/sandbox.py index 0a9a04e48..9fd576950 100644 --- a/benchmarks/sandbox.py +++ b/benchmarks/sandbox.py @@ -44,20 +44,20 @@ class Sandbox(Runner): sandbox = True # select here the external program between LHA, LHAPDF, apfel, pegasus - external = "apfel" + external = "void" # external = "pegasus" # select to plot operators - plot_operator = True + # plot_operator = True - rotate_to_evolution_basis = True + # rotate_to_evolution_basis = True @staticmethod def generate_operators(): ops = { "ev_op_iterations": [1], # "ev_op_max_order": [20], - "mugrid": [[10]], + "mugrid": [[10.0]], # "debug_skip_singlet": [True], } return ops @@ -118,12 +118,12 @@ def lha(self): "FNS": "FFNS", "NfFF": 4, "ModEv": "EXA", - "XIF": np.sqrt(2), - "Q0": np.sqrt(2), + "XIF": float(np.sqrt(2)), + "Q0": float(np.sqrt(2)), "kcThr": 0.0, "kbThr": np.inf, "ktThr": np.inf, - "Qref": np.sqrt(2.0), + "Qref": float(np.sqrt(2.0)), "alphas": 0.35, } self.skip_pdfs = lambda _theory: [ @@ -139,11 +139,11 @@ def lha(self): ] self.run( [theory_updates], - [{"mugrid": [100], "debug_skip_singlet": True}], + [{"mugrid": [100.0], "debug_skip_singlet": True}], ["ToyLH"], ) if __name__ == "__main__": sand = Sandbox() - sand.doit() + sand.lha() diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index e2ea83091..6a75f1437 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -147,6 +147,21 @@ def run_external(self, theory, ocard, pdf): DGLAP result """ # pylint:disable=import-error,import-outside-toplevel + if self.external.lower() == "void": + xgrid = ocard["interpolation_xgrid"] + mugrid = ocard["mugrid"] + labels = ( + br.evol_basis + if self.rotate_to_evolution_basis + else br.flavor_basis_pids + ) + void = { + "target_xgrid": xgrid, + "values": { + mu**2: {pid: [0.0] * len(xgrid) for pid in labels} for mu in mugrid + }, + } + return void if self.external.lower() == "lha": from .external import LHA_utils From efe3c9310014153cd515f4d6ad39ffd07514683d Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 22 Jan 2025 18:01:35 +0200 Subject: [PATCH 55/71] Fix bug in benchmark runner --- src/ekomark/benchmark/runner.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index 6a75f1437..449b96335 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -226,7 +226,7 @@ def log(self, theory, _, pdf, me, ext): # LHA NNLO VFNS needs a special treatment # Valence contains only u and d rotate_to_evolution = None - labels = None + labels = br.flavor_basis_pids qed = theory["QED"] > 0 if self.rotate_to_evolution_basis: if not qed: From ab59da8c520e9b4c9a53f5f68e34d8d8932798ff Mon Sep 17 00:00:00 2001 From: tgiani Date: Tue, 28 Jan 2025 16:37:18 +0100 Subject: [PATCH 56/71] fix bug in qed anomalous dimensions --- .../ekore/src/anomalous_dimensions/unpolarized/spacelike.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index cada98f15..fa841f346 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -129,7 +129,7 @@ pub fn gamma_singlet_qed( // QCD corrections let gamma_qcd_s = gamma_singlet_qcd(order_qcd, c, nf); let gamma_qcd_nsp = gamma_ns_qcd(order_qcd, PID_NSP, c, nf); - for j in 0..order_qcd - 1 { + for j in 0..order_qcd { let s = gamma_qcd_s[j]; gamma_s[1 + j][0] = [ [ @@ -181,7 +181,7 @@ pub fn gamma_valence_qed( // QCD corrections let gamma_qcd_nsv = gamma_ns_qcd(order_qcd, PID_NSV, c, nf); let gamma_qcd_nsm = gamma_ns_qcd(order_qcd, PID_NSM, c, nf); - for j in 0..order_qcd - 1 { + for j in 0..order_qcd { gamma_v[1 + j][0] = [ [gamma_qcd_nsv[j], Complex::::zero()], [Complex::::zero(), gamma_qcd_nsm[j]], From ffc6f701d694cc5b3ae6a49ea26fcdbc4ba12b76 Mon Sep 17 00:00:00 2001 From: tgiani Date: Tue, 28 Jan 2025 16:54:42 +0100 Subject: [PATCH 57/71] same bug in ns sector --- crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index fa841f346..b07488b88 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -79,7 +79,7 @@ pub fn gamma_ns_qed( _ => mode, }; let gamma_qcd = gamma_ns_qcd(order_qcd, qcd_mode, c, nf); - for j in 0..order_qcd - 1 { + for j in 0..order_qcd { gamma_ns[1 + j][0] = gamma_qcd[j]; } // QED corrections From 3a1bf31cc7594e3b1c3847e9d99f88d89e30fada Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 11:16:10 +0200 Subject: [PATCH 58/71] rust: Remove G3p2 from cache --- .../unpolarized/spacelike/as1aem1.rs | 29 +++++++++++++++++-- crates/ekore/src/harmonics/cache.rs | 2 -- 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index 99e469151..db1c0d3dd 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -4,6 +4,7 @@ use num::complex::Complex; use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC, TR, ZETA2, ZETA3}; use crate::harmonics::cache::{Cache, K}; +use std::f64::consts::PI; /// Compute the photon-quark anomalous dimension. /// @@ -112,6 +113,16 @@ pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex { cmplx!(4.0 * TR * (NC as f64), 0.) } +/// Shift for $g_3(N)$ by 2. +/// +/// This is $g_3(N+2) - g_3(N)$. +fn g3_shift(c: &mut Cache) -> Complex { + let n = c.n(); + let S1 = c.get(K::S1); + (6. * (n + 1.) * (2. * n + 1.) * S1 + n * (-PI.powi(2) * (n + 1.).powu(2) - 6. * n)) + / (6. * n.powu(2) * (n + 1.).powu(3)) +} + /// Compute the singlet-like non-singlet anomalous dimension. /// /// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. @@ -128,7 +139,7 @@ pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { let S3p1h = c.get(K::S3ph); let g3N = c.get(K::G3); - let g3Np2 = c.get(K::G3p2); + let g3Np2 = g3N + g3_shift(c); #[rustfmt::skip] let result = 32.0 * ZETA2 * S1h - 32.0 * ZETA2 * S1p1h @@ -159,7 +170,7 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { let S2p1h = c.get(K::S2ph); let S3p1h = c.get(K::S3ph); let g3N = c.get(K::G3); - let g3Np2 = c.get(K::G3p2); + let g3Np2 = g3N + g3_shift(c); #[rustfmt::skip] let result = @@ -302,4 +313,18 @@ mod tests { let me = gamma_nsp(&mut c, NF) + gamma_gq(&mut c, NF) + gamma_phq(&mut c, NF); assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-4); } + + #[test] + fn test_g3_shift() { + for N in [cmplx!(1.23, 0.), cmplx!(5., 0.), cmplx!(2., 1.)] { + let mut c = Cache::new(N); + let mut c2 = Cache::new(N + 2.); + assert_approx_eq_cmplx!( + f64, + g3_shift(&mut c), + c2.get(K::G3) - c.get(K::G3), + epsilon = 1e-6 + ); + } + } } diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs index a9dc6ad25..1e15a0a30 100644 --- a/crates/ekore/src/harmonics/cache.rs +++ b/crates/ekore/src/harmonics/cache.rs @@ -51,7 +51,6 @@ pub enum K { S2ph, S3ph, S1p2, - G3p2, } /// Hold all cached values. @@ -96,7 +95,6 @@ impl Cache { K::S2mh => w2::S2((self.n - 1.) / 2.), K::S3mh => w3::S3((self.n - 1.) / 2.), K::G3 => g_functions::g3(self.n, self.get(K::S1)), - K::G3p2 => g_functions::g3(self.n + 2., self.get(K::S1p2)), K::Sm1e => w1::Sm1e(self.get(K::S1), self.get(K::S1h)), K::Sm1o => w1::Sm1o(self.get(K::S1), self.get(K::S1mh)), K::Sm2e => w2::Sm2e(self.get(K::S2), self.get(K::S2h)), From 044432b7b790a8fda0a66bde2fbde00efef0ccf1 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 11:20:33 +0200 Subject: [PATCH 59/71] rust: Drop recursive definitions from cache --- .../unpolarized/spacelike/as1aem1.rs | 14 +++++++------- crates/ekore/src/harmonics/cache.rs | 9 --------- 2 files changed, 7 insertions(+), 16 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs index db1c0d3dd..ac4596edf 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -3,7 +3,7 @@ use crate::cmplx; use num::complex::Complex; use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC, TR, ZETA2, ZETA3}; -use crate::harmonics::cache::{Cache, K}; +use crate::harmonics::cache::{recursive_harmonic_sum, Cache, K}; use std::f64::consts::PI; /// Compute the photon-quark anomalous dimension. @@ -134,9 +134,9 @@ pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { let S1h = c.get(K::S1h); let S2h = c.get(K::S2h); let S3h = c.get(K::S3h); - let S1p1h = c.get(K::S1ph); - let S2p1h = c.get(K::S2ph); - let S3p1h = c.get(K::S3ph); + let S1p1h = recursive_harmonic_sum(c.get(K::S1mh), (N - 1.) / 2., 1, 1); + let S2p1h = recursive_harmonic_sum(c.get(K::S2mh), (N - 1.) / 2., 1, 2); + let S3p1h = recursive_harmonic_sum(c.get(K::S3mh), (N - 1.) / 2., 1, 3); let g3N = c.get(K::G3); let g3Np2 = g3N + g3_shift(c); @@ -166,9 +166,9 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { let S1h = c.get(K::S1h); let S2h = c.get(K::S2h); let S3h = c.get(K::S3h); - let S1p1h = c.get(K::S1ph); - let S2p1h = c.get(K::S2ph); - let S3p1h = c.get(K::S3ph); + let S1p1h = recursive_harmonic_sum(c.get(K::S1mh), (N - 1.) / 2., 1, 1); + let S2p1h = recursive_harmonic_sum(c.get(K::S2mh), (N - 1.) / 2., 1, 2); + let S3p1h = recursive_harmonic_sum(c.get(K::S3mh), (N - 1.) / 2., 1, 3); let g3N = c.get(K::G3); let g3Np2 = g3N + g3_shift(c); diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs index 1e15a0a30..da0f4800b 100644 --- a/crates/ekore/src/harmonics/cache.rs +++ b/crates/ekore/src/harmonics/cache.rs @@ -46,11 +46,6 @@ pub enum K { Sm21e, /// $S_{-2,1}(N)$ odd moments Sm21o, - /// recursive harmonics - S1ph, - S2ph, - S3ph, - S1p2, } /// Hold all cached values. @@ -103,10 +98,6 @@ impl Cache { K::Sm3o => w3::Sm3o(self.get(K::S3), self.get(K::S3mh)), K::Sm21e => w3::Sm21e(self.n, self.get(K::S1), self.get(K::Sm1e)), K::Sm21o => w3::Sm21o(self.n, self.get(K::S1), self.get(K::Sm1o)), - K::S1ph => recursive_harmonic_sum(self.get(K::S1mh), (self.n - 1.) / 2., 1, 1), - K::S2ph => recursive_harmonic_sum(self.get(K::S2mh), (self.n - 1.) / 2., 1, 2), - K::S3ph => recursive_harmonic_sum(self.get(K::S3mh), (self.n - 1.) / 2., 1, 3), - K::S1p2 => recursive_harmonic_sum(self.get(K::S1), self.n, 2, 1), }; // insert self.m.insert(k, val); From 4507b753af4e7c755d425447308fec7b46dc9d17 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 11:30:54 +0200 Subject: [PATCH 60/71] Avoid setting more couplings in OME --- src/eko/evolution_operator/__init__.py.patch | 72 +++++++++++--------- 1 file changed, 41 insertions(+), 31 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index d97f5c768..a1c71f73a 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,8 +1,8 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index bd1b19d6..de87651c 100644 +index bd1b19d6..f4e2f5ea 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py -@@ -3,16 +3,16 @@ r"""Contains the central operator classes. +@@ -3,120 +3,28 @@ r"""Contains the central operator classes. See :doc:`Operator overview `. """ @@ -17,16 +17,28 @@ index bd1b19d6..de87651c 100644 import numba as nb import numpy as np -from scipy import integrate +- +-import ekore.anomalous_dimensions.polarized.space_like as ad_ps +-import ekore.anomalous_dimensions.unpolarized.space_like as ad_us +-import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut +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 -@@ -32,91 +32,11 @@ from ..matchings import Segment, lepton_number + from .. import basis_rotation as br + from .. import interpolation, mellin + from .. import scale_variations as sv + from ..io.types import EvolutionMethod, OperatorLabel + from ..kernels import ev_method +-from ..kernels import non_singlet as ns +-from ..kernels import non_singlet_qed as qed_ns +-from ..kernels import singlet as s +-from ..kernels import singlet_qed as qed_s +-from ..kernels import valence_qed as qed_v +-from ..matchings import Segment, lepton_number ++from ..matchings import Segment from ..member import OpMember - from ..scale_variations import expanded as sv_expanded - from ..scale_variations import exponentiated as sv_exponentiated -+from .quad_ker import cb_quad_ker_qcd -+from .quad_ker import cb_quad_ker_qed +-from ..scale_variations import expanded as sv_expanded +-from ..scale_variations import exponentiated as sv_exponentiated ++from .quad_ker import cb_quad_ker_qcd, cb_quad_ker_qed logger = logging.getLogger(__name__) @@ -115,7 +127,7 @@ index bd1b19d6..de87651c 100644 spec = [ ("is_singlet", nb.boolean), ("is_QEDsinglet", nb.boolean), -@@ -188,422 +108,6 @@ class QuadKerBase: +@@ -188,422 +96,6 @@ class QuadKerBase: return self.path.prefactor * pj * self.path.jac @@ -538,7 +550,7 @@ index bd1b19d6..de87651c 100644 OpMembers = Dict[OperatorLabel, OpMember] """Map of all operators.""" -@@ -792,49 +296,6 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -792,49 +284,6 @@ class Operator(sv.ScaleVariationModeMixin): """Return the evolution method.""" return ev_method(EvolutionMethod(self.config["method"])) @@ -588,7 +600,7 @@ index bd1b19d6..de87651c 100644 def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( -@@ -857,10 +318,17 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -857,10 +306,32 @@ class Operator(sv.ScaleVariationModeMixin): else: self.op_members[n] = zero.copy() @@ -598,8 +610,23 @@ index bd1b19d6..de87651c 100644 - ): + def update_cfg(self, cfg): + """Adjust integration config.""" ++ # prepare couplings for c + cfg.as1 = self.as_list[1] + cfg.as0 = self.as_list[0] ++ as_list_len = self.as_list.shape[0] ++ as_list_ffi = ekors.ffi.new(f"double[{as_list_len}]", self.as_list.tolist()) ++ cfg.as_list = as_list_ffi ++ cfg.as_list_len = as_list_len ++ a_half_x = self.a_half_list.shape[0] ++ a_half_y = self.a_half_list.shape[1] ++ a_half_len = a_half_x * a_half_y ++ a_half_ffi = ekors.ffi.new( ++ f"double[{a_half_len}]", self.a_half_list.flatten().tolist() ++ ) ++ cfg.a_half = a_half_ffi ++ cfg.a_half_x = a_half_x ++ cfg.a_half_y = a_half_y ++ + if self.order[1] == 0: + cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qcd.address) + else: @@ -610,7 +637,7 @@ index bd1b19d6..de87651c 100644 """Run the integration for each grid point. Parameters -@@ -875,18 +343,75 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -875,18 +346,58 @@ class Operator(sv.ScaleVariationModeMixin): """ column = [] k, logx = log_grid @@ -633,24 +660,7 @@ index bd1b19d6..de87651c 100644 + cfg.is_threshold = self.is_threshold + cfg.mu2_from = self.q2_from + cfg.mu2_to = self.q2_to -+ cfg.alphaem_running=self.alphaem_running -+ -+ # prepare as_list for c -+ as_list_len = self.as_list.shape[0] -+ as_list_ffi = ekors.ffi.new(f"double[{as_list_len}]", self.as_list.tolist()) -+ cfg.as_list = as_list_ffi -+ cfg.as_list_len = as_list_len -+ -+ # prepare a_half for c -+ a_half_x = self.a_half_list.shape[0] -+ a_half_y = self.a_half_list.shape[1] -+ a_half_len = a_half_x * a_half_y -+ a_half_ffi = ekors.ffi.new( -+ f"double[{a_half_len}]", self.a_half_list.flatten().tolist() -+ ) -+ cfg.a_half = a_half_ffi -+ cfg.a_half_x = a_half_x -+ cfg.a_half_y = a_half_y ++ cfg.alphaem_running = self.alphaem_running + + self.update_cfg(cfg) + From 0a4717e31743d6b282854739b857d876c8834ea4 Mon Sep 17 00:00:00 2001 From: tgiani Date: Mon, 3 Feb 2025 11:13:20 +0100 Subject: [PATCH 61/71] unit tests for rust qed anomalous dimensions --- .../unpolarized/spacelike.rs | 136 ++++++++++++++++++ 1 file changed, 136 insertions(+) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index b07488b88..ced5d40a4 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -285,4 +285,140 @@ mod tests { epsilon = 1e-5 ); } + + #[test] + fn test_gamma_valence_qed() { + const NF: u8 = 3; + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + assert_approx_eq_cmplx!( + f64, + gamma_valence_qed(3, 2, &mut c, NF)[3][0][0][0], + cmplx!(459.646893789751, 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_valence_qed(3, 2, &mut c, NF)[3][0][1][1], + cmplx!(437.60340375, 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_valence_qed(3, 2, &mut c, NF)[3][0][1][0], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_valence_qed(3, 2, &mut c, NF)[3][0][0][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + } + + #[test] + fn test_gamma_singlet_qed() { + const NF: u8 = 3; + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][0], + cmplx!(3.857918949669738, 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][2], + cmplx!(-290.42193908689745, 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][3], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][0], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][2], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][3], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][0], + cmplx!(-3.859554320251334, 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][2], + cmplx!(290.4252052962147, 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][3], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][0], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][2], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + assert_approx_eq_cmplx!( + f64, + gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][3], + cmplx!(448.0752570151872, 0.), + epsilon = 1e-5 + ); + } } From 938378a9ca774eaf535c4ccdaffcf7db08ddaa37 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 12:48:44 +0200 Subject: [PATCH 62/71] rust: Add more test macros for complex --- .../unpolarized/spacelike.rs | 30 +++++-------------- crates/ekore/src/util.rs | 13 ++++++++ 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index ced5d40a4..499c66e00 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -200,7 +200,7 @@ pub fn gamma_valence_qed( #[cfg(test)] mod tests { use super::*; - use crate::{assert_approx_eq_cmplx, cmplx}; + use crate::{assert_approx_eq_cmplx, assert_approx_eq_cmplx_2d, cmplx}; use num::complex::Complex; #[test] @@ -292,28 +292,14 @@ mod tests { const N: Complex = cmplx!(2., 0.); let mut c = Cache::new(N); - assert_approx_eq_cmplx!( - f64, - gamma_valence_qed(3, 2, &mut c, NF)[3][0][0][0], - cmplx!(459.646893789751, 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_valence_qed(3, 2, &mut c, NF)[3][0][1][1], - cmplx!(437.60340375, 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( + assert_approx_eq_cmplx_2d!( f64, - gamma_valence_qed(3, 2, &mut c, NF)[3][0][1][0], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_valence_qed(3, 2, &mut c, NF)[3][0][0][1], - cmplx!(0., 0.), + gamma_valence_qed(3, 2, &mut c, NF)[3][0], + [ + [cmplx!(459.646893789751, 0.), cmplx!(0., 0.)], + [cmplx!(0., 0.), cmplx!(437.60340375, 0.)] + ], + 2, epsilon = 1e-5 ); } diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs index b5c812165..c564bcdcc 100644 --- a/crates/ekore/src/util.rs +++ b/crates/ekore/src/util.rs @@ -22,3 +22,16 @@ macro_rules! assert_approx_eq_cmplx { float_cmp::assert_approx_eq!($size, $ref.im, $target.im $(, $set = $val)*); }; } + +/// Shorthand complex matrix comparators. +#[cfg(test)] +#[macro_export] +macro_rules! assert_approx_eq_cmplx_2d { + ($size:ty, $ref:expr, $target:expr, $d:expr $(, $set:ident = $val:expr)*) => { + for j in 0..$d { + for k in 0..$d { + assert_approx_eq_cmplx!($size, $ref[j][k], $target[j][k] $(, $set = $val)*); + } + } + } +} From 79d5437fe0b7918d21aa3f858ab9836ae291ce78 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 13:52:25 +0200 Subject: [PATCH 63/71] ekomark: Fix void+QED --- src/ekomark/benchmark/runner.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index 449b96335..0451d136c 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -150,11 +150,10 @@ def run_external(self, theory, ocard, pdf): if self.external.lower() == "void": xgrid = ocard["interpolation_xgrid"] mugrid = ocard["mugrid"] - labels = ( - br.evol_basis - if self.rotate_to_evolution_basis - else br.flavor_basis_pids - ) + labels = br.flavor_basis_pids + if self.rotate_to_evolution_basis: + labels = br.unified_evol_basis if theory["QED"] > 0 else br.evol_basis + void = { "target_xgrid": xgrid, "values": { From bf5cb424255de79888664ae727ce244d445f3e0c Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 14:11:06 +0200 Subject: [PATCH 64/71] rust: Simplify test_gamma_singlet_qed --- .../unpolarized/spacelike.rs | 122 +++++------------- 1 file changed, 29 insertions(+), 93 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 499c66e00..7458357dc 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -310,100 +310,36 @@ mod tests { const N: Complex = cmplx!(2., 0.); let mut c = Cache::new(N); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][0], - cmplx!(3.857918949669738, 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][1], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][2], - cmplx!(-290.42193908689745, 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][0][3], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][0], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][1], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][2], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][1][3], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][0], - cmplx!(-3.859554320251334, 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][1], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][2], - cmplx!(290.4252052962147, 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][2][3], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][0], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][1], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( - f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][2], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - assert_approx_eq_cmplx!( + assert_approx_eq_cmplx_2d!( f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0][3][3], - cmplx!(448.0752570151872, 0.), + gamma_singlet_qed(3, 2, &mut c, NF)[3][0], + [ + [ + cmplx!(3.857918949669738, 0.), + cmplx!(0., 0.), + cmplx!(-290.42193908689745, 0.), + cmplx!(0., 0.) + ], + [ + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(0., 0.) + ], + [ + cmplx!(-3.859554320251334, 0.), + cmplx!(0., 0.), + cmplx!(290.4252052962147, 0.), + cmplx!(0., 0.) + ], + [ + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(448.0752570151872, 0.) + ] + ], + 4, epsilon = 1e-5 ); } From e8126f4cc8cc03e514171ff9b1bf227e47eca5fe Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 14:12:55 +0200 Subject: [PATCH 65/71] Explicitly hide asX in cb_quad_ker_qed --- src/eko/evolution_operator/quad_ker.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index 50f7b8100..ba111e3bc 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -390,8 +390,8 @@ def cb_quad_ker_qed( areas_x, areas_y, ev_method, - as1, - as0, + _as1, + _as0, ev_op_iterations, ev_op_max_order_qcd, sv_mode, From 01985ec593e5902ad21a7b4358c157a4d082625e Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 16:11:32 +0200 Subject: [PATCH 66/71] rust: Improve ad/u/s unit tests --- .../unpolarized/spacelike.rs | 118 +++++++++--------- crates/ekore/src/util.rs | 11 ++ 2 files changed, 69 insertions(+), 60 deletions(-) diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index 7458357dc..08de553fe 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -200,11 +200,13 @@ pub fn gamma_valence_qed( #[cfg(test)] mod tests { use super::*; - use crate::{assert_approx_eq_cmplx, assert_approx_eq_cmplx_2d, cmplx}; + use crate::{ + assert_approx_eq_cmplx, assert_approx_eq_cmplx_1d, assert_approx_eq_cmplx_2d, cmplx, + }; use num::complex::Complex; #[test] - fn gamma_ns() { + fn test_gamma_ns_qcd() { const NF: u8 = 3; const N: Complex = cmplx!(1., 0.); let mut c = Cache::new(N); @@ -215,32 +217,29 @@ mod tests { epsilon = 1e-14 ); - for i in [0, 1] { - assert_approx_eq_cmplx!( - f64, - gamma_ns_qcd(2, PID_NSM, &mut c, NF)[i], - cmplx!(0., 0.), - epsilon = 2e-6 - ); - } + assert_approx_eq_cmplx_1d!( + f64, + gamma_ns_qcd(2, PID_NSM, &mut c, NF), + [cmplx!(0., 0.); 2], + 2, + epsilon = 2e-6 + ); - for i in 0..3 { - assert_approx_eq_cmplx!( - f64, - gamma_ns_qcd(3, PID_NSM, &mut c, NF)[i], - cmplx!(0., 0.), - epsilon = 2e-4 - ); - } + assert_approx_eq_cmplx_1d!( + f64, + gamma_ns_qcd(3, PID_NSM, &mut c, NF), + [cmplx!(0., 0.); 3], + 3, + epsilon = 2e-4 + ); - for i in 0..3 { - assert_approx_eq_cmplx!( - f64, - gamma_ns_qcd(3, PID_NSV, &mut c, NF)[i], - cmplx!(0., 0.), - epsilon = 8e-4 - ); - } + assert_approx_eq_cmplx_1d!( + f64, + gamma_ns_qcd(3, PID_NSV, &mut c, NF), + [cmplx!(0., 0.); 3], + 3, + epsilon = 8e-4 + ); } #[test] @@ -249,41 +248,32 @@ mod tests { const N: Complex = cmplx!(1., 0.); let mut c = Cache::new(N); - for i in [0, 1] { - for j in [0, 1] { - assert_approx_eq_cmplx!( - f64, - gamma_ns_qed(1, 1, PID_NSM_U, &mut c, NF)[i][j], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - } + // ns- + for pid in [PID_NSM_U, PID_NSM_D] { + assert_approx_eq_cmplx_2d!( + f64, + gamma_ns_qed(1, 1, pid, &mut c, NF), + [[cmplx!(0., 0.); 2]; 2], + 2, + epsilon = 1e-5 + ); } - for i in [0, 1] { - for j in [0, 1] { - assert_approx_eq_cmplx!( - f64, - gamma_ns_qed(1, 1, PID_NSM_D, &mut c, NF)[i][j], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - } + // ns+ + for pid in [PID_NSP_U, PID_NSP_U] { + // as^0 a^0 must be trivial + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, pid, &mut c, NF)[0][0], + cmplx!(0., 0.) + ); + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, pid, &mut c, NF)[0][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); } - - assert_approx_eq_cmplx!( - f64, - gamma_ns_qed(1, 1, PID_NSP_U, &mut c, NF)[0][1], - cmplx!(0., 0.), - epsilon = 1e-5 - ); - - assert_approx_eq_cmplx!( - f64, - gamma_ns_qed(1, 1, PID_NSP_D, &mut c, NF)[0][1], - cmplx!(0., 0.), - epsilon = 1e-5 - ); } #[test] @@ -292,9 +282,13 @@ mod tests { const N: Complex = cmplx!(2., 0.); let mut c = Cache::new(N); + let g = gamma_valence_qed(3, 2, &mut c, NF); + // as^0 a^0 must be trivial + assert_approx_eq_cmplx_2d!(f64, g[0][0], [[cmplx!(0., 0.); 2]; 2], 2); + // reference from Python side assert_approx_eq_cmplx_2d!( f64, - gamma_valence_qed(3, 2, &mut c, NF)[3][0], + g[3][0], [ [cmplx!(459.646893789751, 0.), cmplx!(0., 0.)], [cmplx!(0., 0.), cmplx!(437.60340375, 0.)] @@ -309,10 +303,14 @@ mod tests { const NF: u8 = 3; const N: Complex = cmplx!(2., 0.); let mut c = Cache::new(N); + let g = gamma_singlet_qed(3, 2, &mut c, NF); + // as^0 a^0 must be trivial + assert_approx_eq_cmplx_2d!(f64, g[0][0], [[cmplx!(0., 0.); 4]; 4], 4); + // reference from Python side assert_approx_eq_cmplx_2d!( f64, - gamma_singlet_qed(3, 2, &mut c, NF)[3][0], + g[3][0], [ [ cmplx!(3.857918949669738, 0.), diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs index c564bcdcc..a9a8a818d 100644 --- a/crates/ekore/src/util.rs +++ b/crates/ekore/src/util.rs @@ -23,6 +23,17 @@ macro_rules! assert_approx_eq_cmplx { }; } +/// Shorthand complex list comparators. +#[cfg(test)] +#[macro_export] +macro_rules! assert_approx_eq_cmplx_1d { + ($size:ty, $ref:expr, $target:expr, $d:expr $(, $set:ident = $val:expr)*) => { + for j in 0..$d { + assert_approx_eq_cmplx!($size, $ref[j], $target[j] $(, $set = $val)*); + } + } +} + /// Shorthand complex matrix comparators. #[cfg(test)] #[macro_export] From da131193c1108b1990860b8b01f6285953f0be47 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 3 Feb 2025 18:54:53 +0200 Subject: [PATCH 67/71] Initialize pointers with max size --- src/eko/evolution_operator/__init__.py.patch | 33 ++++++++++---------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index a1c71f73a..aead1680b 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,5 +1,5 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index bd1b19d6..f4e2f5ea 100644 +index bd1b19d6..bbc9d682 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -3,120 +3,28 @@ r"""Contains the central operator classes. @@ -600,7 +600,7 @@ index bd1b19d6..f4e2f5ea 100644 def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( -@@ -857,10 +306,32 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -857,10 +306,29 @@ class Operator(sv.ScaleVariationModeMixin): else: self.op_members[n] = zero.copy() @@ -614,16 +614,13 @@ index bd1b19d6..f4e2f5ea 100644 + cfg.as1 = self.as_list[1] + cfg.as0 = self.as_list[0] + as_list_len = self.as_list.shape[0] -+ as_list_ffi = ekors.ffi.new(f"double[{as_list_len}]", self.as_list.tolist()) -+ cfg.as_list = as_list_ffi ++ for j, a in enumerate(self.as_list.tolist()): ++ cfg.as_list[j] = a + cfg.as_list_len = as_list_len + a_half_x = self.a_half_list.shape[0] + a_half_y = self.a_half_list.shape[1] -+ a_half_len = a_half_x * a_half_y -+ a_half_ffi = ekors.ffi.new( -+ f"double[{a_half_len}]", self.a_half_list.flatten().tolist() -+ ) -+ cfg.a_half = a_half_ffi ++ for j, a in enumerate(self.a_half_list.flatten().tolist()): ++ cfg.a_half[j] = a + cfg.a_half_x = a_half_x + cfg.a_half_y = a_half_y + @@ -637,7 +634,7 @@ index bd1b19d6..f4e2f5ea 100644 """Run the integration for each grid point. Parameters -@@ -875,18 +346,58 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -875,18 +343,62 @@ class Operator(sv.ScaleVariationModeMixin): """ column = [] k, logx = log_grid @@ -661,6 +658,14 @@ index bd1b19d6..f4e2f5ea 100644 + cfg.mu2_from = self.q2_from + cfg.mu2_to = self.q2_to + cfg.alphaem_running = self.alphaem_running ++ a_half_ffi = ekors.ffi.new("double[]", 2 * self.config["ev_op_iterations"]) ++ cfg.a_half = a_half_ffi ++ as_list_ffi = ekors.ffi.new("double[]", self.config["ev_op_iterations"] + 1) ++ cfg.as_list = as_list_ffi ++ pd = self.int_disp.polynomial_degree + 1 ++ max_area_len = (2 * pd - 1) * (pd + 2) ++ areas_ffi = ekors.ffi.new("double[]", max_area_len) ++ cfg.areas = areas_ffi + + self.update_cfg(cfg) + @@ -675,12 +680,8 @@ index bd1b19d6..f4e2f5ea 100644 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 ++ for j, x in enumerate(curareas.flatten().tolist()): ++ cfg.areas[j] = x + cfg.areas_x = curareas.shape[0] + cfg.areas_y = curareas.shape[1] # iterate sectors From 2ef7a2979ee8cf8eeacaa2ac34c7304f9b615ff4 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 4 Feb 2025 13:32:20 +0200 Subject: [PATCH 68/71] Add max_areas_shape to interpolation --- src/eko/interpolation.py | 5 +++++ tests/eko/test_interpolation.py | 19 ++++++++++++++++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/eko/interpolation.py b/src/eko/interpolation.py index 2948f028d..526bd5c76 100644 --- a/src/eko/interpolation.py +++ b/src/eko/interpolation.py @@ -569,6 +569,11 @@ def __iter__(self): # return iter(self.basis) yield from self.basis + def max_areas_shape(self) -> tuple[int, int]: + """Maximum dimensions of the polynomial areas.""" + # 2. dim: xmin, xmax, 1+degree coefficients + return (max(len(bf.areas) for bf in self), 2 + 1 + self.polynomial_degree) + def __getitem__(self, item): return self.basis[item] diff --git a/tests/eko/test_interpolation.py b/tests/eko/test_interpolation.py index b3d032415..9bb5d3de7 100644 --- a/tests/eko/test_interpolation.py +++ b/tests/eko/test_interpolation.py @@ -133,7 +133,7 @@ def test_get_interpolation(self): inter_x = interpolation.InterpolatorDispatcher(xg, 1, False) i = inter_x.get_interpolation(xg.raw) np.testing.assert_array_almost_equal(i, np.eye(len(xg))) - # .75 is exactly inbetween + # .75 is exactly in between i = inter_x.get_interpolation([0.75]) np.testing.assert_array_almost_equal(i, [[0.5, 0.5]]) @@ -163,6 +163,23 @@ def test_math(self): inter_N = interpolation.InterpolatorDispatcher(xgrid, poly_deg, mode_N=True) check_correspondence_interpolators(inter_x, inter_N) + def test_max_areas_shape(self): + for ng in [3, 6, 7, 8, 10, 11, 14]: + for pd in range(1, ng // 2): + xgrid = interpolation.XGrid(np.linspace(0.1, 1, ng), False) + interp = interpolation.InterpolatorDispatcher(xgrid, pd) + max_shape = interp.max_areas_shape() + # all areas must respect the maximum + for bf in interp: + s = bf.areas_to_const().shape + assert s[0] <= max_shape[0], f"{ng=},{pd=},k={bf.poly_number}" + assert s[1] <= max_shape[1], f"{ng=},{pd=},k={bf.poly_number}" + # the end must be the maximum, because we are leaning upwards + end = ng - pd - 1 + ms = interp[end].areas_to_const().shape + assert ms[0] == max_shape[0], f"{ng=},{pd=},{end=}" + assert ms[1] == max_shape[1], f"{ng=},{pd=},{end=}" + class TestBasisFunction: def test_init(self): From e956e0a916a6db17dced7095611dae6fdcc58548 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 4 Feb 2025 13:36:14 +0200 Subject: [PATCH 69/71] Use max_areas_shape for pointer init --- src/eko/evolution_operator/__init__.py.patch | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index aead1680b..da1547d28 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,5 +1,5 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index bd1b19d6..bbc9d682 100644 +index bd1b19d6..3fa63a70 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -3,120 +3,28 @@ r"""Contains the central operator classes. @@ -662,8 +662,8 @@ index bd1b19d6..bbc9d682 100644 + cfg.a_half = a_half_ffi + as_list_ffi = ekors.ffi.new("double[]", self.config["ev_op_iterations"] + 1) + cfg.as_list = as_list_ffi -+ pd = self.int_disp.polynomial_degree + 1 -+ max_area_len = (2 * pd - 1) * (pd + 2) ++ max_areas_shape = self.int_disp.max_areas_shape ++ max_area_len = max_areas_shape[0] * max_areas_shape[1] + areas_ffi = ekors.ffi.new("double[]", max_area_len) + cfg.areas = areas_ffi + From 74fd52236f5d001e1ac38f61c1f63a087bfc3ecb Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 4 Feb 2025 14:24:50 +0200 Subject: [PATCH 70/71] Call max_areas_shape for pointer init --- src/eko/evolution_operator/__init__.py.patch | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index da1547d28..576c34ae6 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,5 +1,5 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index bd1b19d6..3fa63a70 100644 +index bd1b19d6..24f99075 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -3,120 +3,28 @@ r"""Contains the central operator classes. @@ -662,7 +662,7 @@ index bd1b19d6..3fa63a70 100644 + cfg.a_half = a_half_ffi + as_list_ffi = ekors.ffi.new("double[]", self.config["ev_op_iterations"] + 1) + cfg.as_list = as_list_ffi -+ max_areas_shape = self.int_disp.max_areas_shape ++ max_areas_shape = self.int_disp.max_areas_shape() + max_area_len = max_areas_shape[0] * max_areas_shape[1] + areas_ffi = ekors.ffi.new("double[]", max_area_len) + cfg.areas = areas_ffi From 2bdad40d2473f294e77203973a2496d75bf72d0f Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 6 Feb 2025 13:44:06 +0200 Subject: [PATCH 71/71] Spell nf in apfel VFNS bench --- benchmarks/apfel_bench.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index 2ef8f21a0..c51cae83b 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -286,6 +286,8 @@ class BenchmarkVFNS_qed(ApfelBenchmark): "kbThr": 1.0, "ktThr": 1.0, "Q0": 1.25, + "nf0": 3, + "nfref": 5, "alphas": 0.118000, "alphaqed": 0.007496, }