setting up input vectors for cb_quad_ker_qed

.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]

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,
     }
 }

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(