Rustify QED

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,
@@ -159,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"])
@@ -209,6 +202,8 @@ class BenchmarkFFNS_qed(ApfelBenchmark):
             # "TRN",
         ],
         "NfFF": 5,
+        "nf0": 5,
+        "nfref": 5,
         "kcThr": 0.0,
         "kbThr": 0.0,
         "ktThr": np.inf,
@@ -319,7 +314,7 @@ def benchmark_plain(self, pto, qed):
     # 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)

crates/eko/src/lib.rs

@@ -7,13 +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<f64>,
     im: Vec<f64>,
 }
 
-/// Map tensors to c-ordered list
+/// Map tensor with shape (o,d,d) to c-ordered list.
+///
+/// This is needed for the QCD singlet.
 fn unravel<const DIM: usize>(res: Vec<[[Complex<f64>; DIM]; DIM]>, order_qcd: usize) -> RawCmplx {
     let mut target = RawCmplx {
         re: Vec::<f64>::new(),
@@ -30,14 +32,63 @@ fn unravel<const DIM: usize>(res: Vec<[[Complex<f64>; DIM]; DIM]>, order_qcd: us
     target
 }
 
-/// QCD intergration kernel inside quad.
+/// Map tensor with shape (o,o',d,d) to c-ordered list.
+///
+/// This is needed for the QED singlet and valence.
+fn unravel_qed<const DIM: usize>(
+    res: Vec<Vec<[[Complex<f64>; DIM]; DIM]>>,
+    order_qcd: usize,
+    order_qed: usize,
+) -> RawCmplx {
+    let mut target = RawCmplx {
+        re: Vec::<f64>::new(),
+        im: Vec::<f64>::new(),
+    };
+    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);
+                    target.im.push(el.im);
+                }
+            }
+        }
+    }
+    target
+}
+
+/// 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<Vec<Complex<f64>>>, order_qcd: usize, order_qed: usize) -> RawCmplx {
+    let mut target = RawCmplx {
+        re: Vec::<f64>::new(),
+        im: Vec::<f64>::new(),
+    };
+    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);
+        }
+    }
+    target
+}
+
+/// 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);
-    let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0);
+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)
+        || (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);
     let jac = path.jac() * path.prefactor();
@@ -70,14 +121,41 @@ 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),
-            args.order_qcd,
-        );
+        if args.order_qed > 0 {
+            let gamma_singlet_qed =
+                ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qed;
+            raw = unravel_qed(
+                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.order_qed > 0 {
+        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 {
@@ -100,6 +178,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,
@@ -109,7 +188,6 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 {
         args.areas,
         args.areas_x,
         args.areas_y,
-        args.L,
         args.method_num,
         args.as1,
         args.as0,
@@ -118,50 +196,68 @@ 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.as_list_len,
+        args.mu2_from,
+        args.mu2_to,
+        args.a_half,
+        args.a_half_x,
+        args.a_half_y,
+        args.alphaem_running,
     )
 }
 
 /// Python callback signature
-type PyQuadKerQCDT = unsafe extern "C" fn(
-    *const f64,
-    *const f64,
-    f64,
-    f64,
-    f64,
-    f64,
-    usize,
-    bool,
-    u16,
-    u16,
-    u8,
-    bool,
-    f64,
-    *const f64,
-    u8,
-    u8,
-    f64,
-    u8,
-    f64,
-    f64,
-    u8,
-    u8,
-    u8,
-    bool,
-    f64,
+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
 #[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,
     pub mode1: u16,
     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,
@@ -177,6 +273,15 @@ 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 as_list_len: u8,
+    pub mu2_from: f64,
+    pub mu2_to: f64,
+    pub a_half: *const f64,
+    pub a_half_x: u8,
+    pub a_half_y: u8,
+    pub alphaem_running: bool,
 }
 
 /// Empty placeholder function for python callback.
@@ -191,6 +296,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,
@@ -200,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,
@@ -209,21 +314,30 @@ pub unsafe extern "C" fn my_py(
     _sv_mode_num: u8,
     _is_threshold: bool,
     _lsv: f64,
+    _as_list: *const f64,
+    _as_list_len: u8,
+    _mu2_from: f64,
+    _mu2_to: f64,
+    _a_half: *const f64,
+    _a_half_x: u8,
+    _a_half_y: u8,
+    _alphaem_running: bool,
 ) -> f64 {
     0.
 }
 
 /// Return empty additional arguments.
 ///
 /// This is required to make the arguments part of the API, otherwise it won't be added to the compiled
-/// package (since it does not appear in the signature of `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,
         mode1: 0,
         is_polarized: false,
@@ -245,5 +359,13 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs {
         is_threshold: false,
         is_ome: false,
         Lsv: 0.,
+        as_list: [].as_ptr(),
+        as_list_len: 0,
+        mu2_from: 0.,
+        mu2_to: 0.,
+        a_half: [].as_ptr(),
+        a_half_x: 0,
+        a_half_y: 0,
+        alphaem_running: false,
     }
 }