Rustify ome.u.s.{as1, as2}

crates/dekoder/src/eko.rs

@@ -92,7 +92,7 @@ impl EKO {
     pub fn write(&self, dst: PathBuf) -> Result<()> {
         self.assert_working_dir()?;
         // create writer
-        let dst_file = File::create(&dst)?;
+        let dst_file = File::create(dst)?;
         let dst_file = BufWriter::with_capacity(TAR_WRITER_CAPACITY, dst_file);
         let mut ar = tar::Builder::new(dst_file);
         // do it!
@@ -101,7 +101,7 @@ impl EKO {
 
     /// Extract tar file from `src` to `dst`.
     pub fn extract(src: PathBuf, dst: PathBuf) -> Result<Self> {
-        let mut ar = tar::Archive::new(File::open(&src)?);
+        let mut ar = tar::Archive::new(File::open(src)?);
         ar.unpack(&dst)?;
         Self::load_opened(dst)
     }

crates/eko/src/lib.rs

@@ -1,11 +1,35 @@
 //! Interface to the eko Python package.
 
 use ekore::harmonics::cache::Cache;
+use num::Complex;
 use std::ffi::c_void;
 
 pub mod bib;
 pub mod mellin;
 
+/// Wrapper to pass arguments back to Python
+struct RawCmplx {
+    re: Vec<f64>,
+    im: Vec<f64>,
+}
+
+/// Map tensors to c-ordered list
+fn unravel<const DIM: usize>(res: Vec<[[Complex<f64>; DIM]; DIM]>, order_qcd: usize) -> RawCmplx {
+    let mut target = RawCmplx {
+        re: Vec::<f64>::new(),
+        im: Vec::<f64>::new(),
+    };
+    for obj in res.iter().take(order_qcd) {
+        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
+}
+
 /// QCD intergration kernel inside quad.
 ///
 /// # Safety
@@ -14,44 +38,66 @@ pub mod mellin;
 pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 {
     let args = *(rargs as *mut QuadQCDargs);
     let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0);
-    // prepare gamma
+    // prepare Mellin stuff
     let path = mellin::TalbotPath::new(u, args.logx, is_singlet);
     let jac = path.jac() * path.prefactor();
     let mut c = Cache::new(path.n());
-    let mut re = Vec::<f64>::new();
-    let mut im = Vec::<f64>::new();
-    if is_singlet {
-        let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
+    let mut raw = RawCmplx {
+        re: Vec::<f64>::new(),
+        im: Vec::<f64>::new(),
+    };
+
+    if args.is_ome {
+        if is_singlet {
+            raw = unravel(
+                ekore::operator_matrix_elements::unpolarized::spacelike::A_singlet(
+                    args.order_qcd,
+                    &mut c,
+                    args.nf,
+                    args.L,
+                ),
+                args.order_qcd,
+            );
+        } else {
+            raw = unravel(
+                ekore::operator_matrix_elements::unpolarized::spacelike::A_non_singlet(
+                    args.order_qcd,
+                    &mut c,
+                    args.nf,
+                    args.L,
+                ),
+                args.order_qcd,
+            );
+        }
+    } else if is_singlet {
+        raw = unravel(
+            ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
+                args.order_qcd,
+                &mut c,
+                args.nf,
+            ),
             args.order_qcd,
-            &mut c,
-            args.nf,
         );
-        for gamma_s in res.iter().take(args.order_qcd) {
-            for col in gamma_s.iter().take(2) {
-                for el in col.iter().take(2) {
-                    re.push(el.re);
-                    im.push(el.im);
-                }
-            }
-        }
     } 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,
         );
         for el in res.iter().take(args.order_qcd) {
-            re.push(el.re);
-            im.push(el.im);
+            raw.re.push(el.re);
+            raw.im.push(el.im);
         }
     }
+
     // pass on
     (args.py)(
-        re.as_ptr(),
-        im.as_ptr(),
-        c.n.re,
-        c.n.im,
+        raw.re.as_ptr(),
+        raw.im.as_ptr(),
+        c.n().re,
+        c.n().im,
         jac.re,
         jac.im,
         args.order_qcd,
@@ -72,6 +118,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 {
         args.ev_op_max_order_qcd,
         args.sv_mode_num,
         args.is_threshold,
+        args.Lsv,
     )
 }
 
@@ -101,6 +148,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn(
     u8,
     u8,
     bool,
+    f64,
 ) -> f64;
 
 /// Additional integration parameters
@@ -128,6 +176,8 @@ pub struct QuadQCDargs {
     pub ev_op_max_order_qcd: u8,
     pub sv_mode_num: u8,
     pub is_threshold: bool,
+    pub is_ome: bool,
+    pub Lsv: f64,
 }
 
 /// Empty placeholder function for python callback.
@@ -159,14 +209,15 @@ pub unsafe extern "C" fn my_py(
     _ev_op_max_order_qcd: u8,
     _sv_mode_num: u8,
     _is_threshold: bool,
+    _lsv: f64,
 ) -> f64 {
     0.
 }
 
 /// Return empty additional arguments.
 ///
 /// This is required to make the arguments part of the API, otherwise it won't be added to the compiled
-/// package (since it does not appear in the signature of `quad_ker_qcd`).
+/// package (since it does not appear in the signature of `rust_quad_ker_qcd`).
 ///
 /// # Safety
 /// This is the connection from and back to Python, so we don't know what is on the other side.
@@ -193,5 +244,7 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs {
         ev_op_max_order_qcd: 0,
         sv_mode_num: 0,
         is_threshold: false,
+        is_ome: false,
+        Lsv: 0.,
     }
 }

crates/ekore/src/anomalous_dimensions.rs

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

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

@@ -1,4 +1,4 @@
-//! LO QCD
+//! |LO| |QCD|.
 
 use num::complex::Complex;
 
@@ -9,7 +9,7 @@ use crate::harmonics::cache::{Cache, K};
 ///
 /// Implements Eq. (3.4) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
 pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0));
     CF * gamma
@@ -19,7 +19,7 @@ pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
 pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N + 2.0));
     2.0 * TR * 2.0 * (nf as f64) * gamma
 }
@@ -28,7 +28,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
 pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N - 1.0));
     2.0 * CF * gamma
 }
@@ -37,7 +37,7 @@ pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
 pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let gamma = S1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0);
     CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * (nf as f64)

crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs

@@ -1,4 +1,4 @@
-//! NLO QCD
+//! |NLO| |QCD|.
 
 use ::num::complex::Complex;
 use std::f64::consts::LN_2;
@@ -10,7 +10,7 @@ use crate::harmonics::cache::{Cache, K};
 ///
 /// Implements Eq. (3.6) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
 pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let Sp1m = c.get(K::S1mh);
@@ -50,7 +50,7 @@ pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.5) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
 pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let Sp1p = c.get(K::S1h);
@@ -90,7 +90,7 @@ pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.6) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let gqqps1_nfcf = (-4. * (2. + N * (5. + N)) * (4. + N * (4. + N * (7. + 5. * N))))
         / ((-1. + N) * N.powu(3) * (1. + N).powu(3) * (2. + N).powu(2));
     2.0 * TR * (nf as f64) * CF * gqqps1_nfcf
@@ -100,7 +100,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.7) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let Sp2p = c.get(K::S2h);
@@ -128,7 +128,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.8) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let Sp2p = c.get(K::S2h);
@@ -164,7 +164,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.9) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let Sp1p = c.get(K::S1h);
     let Sp2p = c.get(K::S2h);

crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs

@@ -1,4 +1,4 @@
-//! NNLO QCD
+//! |NNLO| |QCD|.
 
 use ::num::complex::Complex;
 use num::traits::Pow;
@@ -11,7 +11,7 @@ use crate::harmonics::cache::{Cache, K};
 ///
 /// Implements Eq. (3.8) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
 pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);
@@ -68,7 +68,7 @@ pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.7) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
 pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);
@@ -125,7 +125,7 @@ pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.9) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
 pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);
@@ -162,7 +162,7 @@ pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.10) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);
@@ -218,7 +218,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.11) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);
@@ -275,7 +275,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.12) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);
@@ -348,7 +348,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
 ///
 /// Implements Eq. (3.13) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
 pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
-    let N = c.n;
+    let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
     let S3 = c.get(K::S3);

crates/ekore/src/constants.rs

@@ -33,11 +33,11 @@ pub const ZETA3: f64 = 1.2020569031595942;
 /// $\zeta(4) = \pi^4 / 90$.
 pub const ZETA4: f64 = 1.082323233711138;
 
-/// singlet-like non-singlet PID
+/// singlet-like non-singlet |PID|.
 pub const PID_NSP: u16 = 10101;
 
-/// valence-like non-singlet PID
+/// valence-like non-singlet |PID|.
 pub const PID_NSM: u16 = 10201;
 
-/// non-singlet all-valence PID
+/// non-singlet all-valence |PID|.
 pub const PID_NSV: u16 = 10200;