Rustify ad.p.s.{as1,as2}

.github/workflows/lha_bot_rust.yml

@@ -17,11 +17,6 @@ jobs:
   lhabench:
     name: LHA paper Benchmarks
     runs-on: ubuntu-latest
-    # container:
-    #   image: ghcr.io/nnpdf/bench-evol:v2
-    #   credentials:
-    #     username: ${{ github.repository_owner }}
-    #     password: ${{ secrets.GITHUB_TOKEN }}
 
     strategy:
       matrix:
@@ -30,9 +25,6 @@ jobs:
 
     steps:
       - uses: actions/checkout@v2
-        # with:
-        #   # tags needed for dynamic versioning
-        #   fetch-depth: 0
       - name: Set up Python ${{ matrix.python-version }} 🐍
         id: setup-python
         uses: actions/setup-python@v5
@@ -53,5 +45,4 @@ jobs:
           ./rustify.sh
           poe compile
           poe lha -m "nnlo and sv"
-          # TODO wait for polarized to reactivate
-          # poe lha -m "ffns_pol and sv"
+          poe lha -m "ffns_pol and sv"

crates/eko/src/lib.rs

@@ -70,22 +70,21 @@ 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(
-            ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
-                args.order_qcd,
-                &mut c,
-                args.nf,
-            ),
+            gamma_singlet_qcd(args.order_qcd, &mut c, args.nf),
             args.order_qcd,
         );
     } else {
         // we can not do 1D
-        let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd(
-            args.order_qcd,
-            args.mode0,
-            &mut c,
-            args.nf,
-        );
+        let gamma_ns_qcd = match args.is_polarized {
+            true => ekore::anomalous_dimensions::polarized::spacelike::gamma_ns_qcd,
+            false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd,
+        };
+        let res = gamma_ns_qcd(args.order_qcd, args.mode0, &mut c, args.nf);
         for el in res.iter().take(args.order_qcd) {
             raw.re.push(el.re);
             raw.im.push(el.im);

crates/ekore/src/anomalous_dimensions.rs

@@ -1,3 +1,4 @@
 //! The anomalous dimensions for |DGLAP| evolution.
 
+pub mod polarized;
 pub mod unpolarized;

crates/ekore/src/anomalous_dimensions/polarized.rs

@@ -0,0 +1,3 @@
+//! The polarized anomalous dimensions for space-like kinematics.
+
+pub mod spacelike;

crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs

@@ -0,0 +1,57 @@
+//! The polarized, space-like anomalous dimensions at various couplings power.
+
+use crate::constants::{PID_NSM, PID_NSP, PID_NSV};
+use crate::harmonics::cache::Cache;
+use num::complex::Complex;
+use num::Zero;
+pub mod as1;
+pub mod as2;
+// pub mod as3;
+
+/// Compute the tower of the non-singlet anomalous dimensions.
+pub fn gamma_ns_qcd(order_qcd: usize, mode: u16, c: &mut Cache, nf: u8) -> Vec<Complex<f64>> {
+    let mut gamma_ns = vec![Complex::<f64>::zero(); order_qcd];
+    gamma_ns[0] = as1::gamma_ns(c, nf);
+    // NLO and beyond
+    if order_qcd >= 2 {
+        let gamma_ns_1 = match mode {
+            PID_NSP => as2::gamma_nsp(c, nf),
+            // To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
+            PID_NSM | PID_NSV => as2::gamma_nsm(c, nf),
+            _ => panic!("Unkown non-singlet sector element"),
+        };
+        gamma_ns[1] = gamma_ns_1
+    }
+    // // NNLO and beyond
+    // if order_qcd >= 3 {
+    //     let gamma_ns_2 = match mode {
+    //         PID_NSP => as3::gamma_nsp(c, nf),
+    //         PID_NSM => as3::gamma_nsm(c, nf),
+    //         PID_NSV => as3::gamma_nsv(c, nf),
+    //         _ => panic!("Unkown non-singlet sector element"),
+    //     };
+    //     gamma_ns[2] = gamma_ns_2
+    // }
+    gamma_ns
+}
+
+/// Compute the tower of the singlet anomalous dimension matrices.
+pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Complex<f64>; 2]; 2]> {
+    let mut gamma_S = vec![
+        [
+            [Complex::<f64>::zero(), Complex::<f64>::zero()],
+            [Complex::<f64>::zero(), Complex::<f64>::zero()]
+        ];
+        order_qcd
+    ];
+    gamma_S[0] = as1::gamma_singlet(c, nf);
+    // NLO and beyond
+    if order_qcd >= 2 {
+        gamma_S[1] = as2::gamma_singlet(c, nf);
+    }
+    // // NNLO and beyond
+    // if order_qcd >= 3 {
+    //     gamma_S[2] = as3::gamma_singlet(c, nf);
+    // }
+    gamma_S
+}

crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs

@@ -0,0 +1,85 @@
+//! |LO| |QCD|.
+
+use num::complex::Complex;
+
+use super::super::super::unpolarized::spacelike::as1::gamma_ns as unpol;
+use crate::constants::{CA, CF, TR};
+use crate::harmonics::cache::{Cache, K};
+
+/// Compute the non-singlet anomalous dimension.
+///
+/// Identical to the unpolarized counterpart.
+pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex<f64> {
+    unpol(c, nf)
+}
+
+/// Compute the quark-gluon anomalous dimension.
+///
+/// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
+pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n();
+    let gamma = -(N - 1.) / N / (N + 1.);
+    2.0 * TR * 2.0 * (nf as f64) * gamma
+}
+
+/// Compute the gluon-quark anomalous dimension.
+///
+/// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
+pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
+    let N = c.n();
+    let gamma = -(N + 2.) / N / (N + 1.);
+    2.0 * CF * gamma
+}
+
+/// Compute the gluon-gluon anomalous dimension.
+///
+/// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
+pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n();
+    let S1 = c.get(K::S1);
+    let gamma = -S1 + 2. / N / (N + 1.);
+    CA * (-4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * (nf as f64)
+}
+
+/// Compute the singlet anomalous dimension matrix.
+pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+    let gamma_qq = gamma_ns(c, nf);
+    [
+        [gamma_qq, gamma_qg(c, nf)],
+        [gamma_gq(c, nf), gamma_gg(c, nf)],
+    ]
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::harmonics::cache::Cache;
+    use crate::{assert_approx_eq_cmplx, cmplx};
+    use num::complex::Complex;
+    use num::Zero;
+    const NF: u8 = 5;
+
+    #[test]
+    fn quark_momentum_conservation() {
+        const N: Complex<f64> = cmplx!(2., 0.);
+        let mut c = Cache::new(N);
+        let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF);
+        assert_approx_eq_cmplx!(f64, me, cmplx!((4. * CF) / 3., 0.), epsilon = 1e-12);
+    }
+
+    #[test]
+    fn gluon_momentum_conservation() {
+        const N: Complex<f64> = cmplx!(2., 0.);
+        let mut c = Cache::new(N);
+        let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF);
+        assert_approx_eq_cmplx!(f64, me, cmplx!(3. + (NF as f64) / 3., 0.), epsilon = 1e-12);
+    }
+
+    #[test]
+    fn qg_helicity_conservation() {
+        const N: Complex<f64> = cmplx!(1., 0.);
+        let mut c = Cache::new(N);
+        let me = gamma_qg(&mut c, NF);
+        assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
+    }
+}