Add fhmruvv option to qed

src/eko/evolution_operator/__init__.py

@@ -312,6 +312,7 @@ def quad_ker(
             sv_mode,
             is_threshold,
             n3lo_ad_variation,
+            use_fhmruvv,
         )
 
     # recombine everything
@@ -461,6 +462,7 @@ def quad_ker_qed(
     sv_mode,
     is_threshold,
     n3lo_ad_variation,
+    use_fhmruvv,
 ):
     """Raw evolution kernel inside quad.
 
@@ -501,7 +503,9 @@ def quad_ker_qed(
     is_threshold : boolean
         is this an itermediate threshold operator?
     n3lo_ad_variation : tuple
-        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)``
+        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
+    use_fhmruvv : bool
+        if True use the |FHMRUVV| |N3LO| anomalous dimensions
 
     Returns
     -------
@@ -510,7 +514,9 @@ def quad_ker_qed(
     """
     # compute the actual evolution kernel for QEDxQCD
     if ker_base.is_QEDsinglet:
-        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation)
+        gamma_s = ad_us.gamma_singlet_qed(
+            order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+        )
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
@@ -538,7 +544,9 @@ def quad_ker_qed(
             ) @ 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)
+        gamma_v = ad_us.gamma_valence_qed(
+            order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+        )
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_v = sv.exponentiated.gamma_variation_qed(
@@ -563,7 +571,9 @@ def quad_ker_qed(
             ) @ np.ascontiguousarray(ker)
         ker = select_QEDvalence_element(ker, mode0, mode1)
     else:
-        gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf)
+        gamma_ns = ad_us.gamma_ns_qed(
+            order, mode0, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+        )
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation_qed(

src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py

@@ -137,20 +137,24 @@ def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
 
 
 @nb.njit(cache=True)
-def gamma_ns_qed(order, mode, n, nf):
+def gamma_ns_qed(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False):
     r"""
     Compute the grid of the QED non-singlet anomalous dimensions.
 
     Parameters
     ----------
-        order : tuple(int,int)
-            perturbative orders
-        mode : 10102 | 10103 | 10202 | 10203
-            sector identifier
-        n : complex
-            Mellin variable
-        nf : int
-            Number of active flavors
+    order : tuple(int,int)
+        perturbative orders
+    mode : 10102 | 10103 | 10202 | 10203
+        sector identifier
+    n : complex
+        Mellin variable
+    nf : int
+        Number of active flavors
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
+    use_fhmruvv: bool
+        if True use the |FHMRUVV| N3LO anomalous dimensions
 
     Returns
     -------
@@ -179,10 +183,20 @@ def gamma_ns_qed(order, mode, n, nf):
         elif mode in [10202, 10203]:
             gamma_ns[3, 0] = as3.gamma_nsm(n, nf, cache)
     if order[0] >= 4:
-        if mode in [10102, 10103]:
-            gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache)
-        elif mode in [10202, 10203]:
-            gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache)
+        if use_fhmruvv:
+            if mode in [10102, 10103]:
+                gamma_ns[4, 0] = as4.fhmruvv.gamma_nsp(
+                    n, nf, cache, n3lo_ad_variation[4]
+                )
+            elif mode in [10202, 10203]:
+                gamma_ns[4, 0] = as4.fhmruvv.gamma_nsm(
+                    n, nf, cache, n3lo_ad_variation[5]
+                )
+        else:
+            if mode in [10102, 10103]:
+                gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache)
+            elif mode in [10202, 10203]:
+                gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache)
     return gamma_ns
 
 
@@ -278,7 +292,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache):
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
+def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
     r"""
     Compute the grid of the QED singlet anomalous dimensions matrices.
 
@@ -291,7 +305,9 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
     nf : int
         Number of active flavors
     n3lo_ad_variation : tuple
-        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)``
+        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
+    use_fhmruvv: bool
+        if True use the |FHMRUVV| N3LO anomalous dimensions
 
     Returns
     -------
@@ -311,12 +327,17 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
     if order[0] >= 3:
         gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, cache)
     if order[0] >= 4:
-        gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation)
+        if use_fhmruvv:
+            gamma_s[4, 0] = as4.fhmruvv.gamma_singlet_qed(
+                n, nf, cache, n3lo_ad_variation
+            )
+        else:
+            gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation)
     return gamma_s
 
 
 @nb.njit(cache=True)
-def gamma_valence_qed(order, n, nf):
+def gamma_valence_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
     r"""
     Compute the grid of the QED valence anomalous dimensions matrices.
 
@@ -328,6 +349,11 @@ def gamma_valence_qed(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
+    use_fhmruvv: bool
+        if True use the |FHMRUVV| N3LO anomalous dimensions
+
 
     Returns
     -------
@@ -347,5 +373,10 @@ def gamma_valence_qed(order, n, nf):
     if order[0] >= 3:
         gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, cache)
     if order[0] >= 4:
-        gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache)
+        if use_fhmruvv:
+            gamma_v[4, 0] = as4.fhmruvv.gamma_valence_qed(
+                n, nf, cache, n3lo_ad_variation
+            )
+        else:
+            gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache)
     return gamma_v

src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py

@@ -56,3 +56,88 @@ def gamma_singlet(N, nf, cache, variation):
         np.complex_,
     )
     return gamma_S_0
+
+
+@nb.njit(cache=True)
+def gamma_singlet_qed(N, nf, cache, variation):
+    r"""Compute the leading-order singlet anomalous dimension matrix for the unified evolution basis.
+
+    .. math::
+        \\gamma_S^{(3,0)} = \\left(\begin{array}{cccc}
+        \\gamma_{gg}^{(3,0)} & 0 & \\gamma_{gq}^{(3,0)} & 0\\
+        0 & 0 & 0 & 0 \\
+        \\gamma_{qg}^{(3,0)} & 0 & \\gamma_{qq}^{(3,0)} & 0 \\
+        0 & 0 & 0 & \\gamma_{qq}^{(3,0)} \\
+        \\end{array}\right)
+
+    Parameters
+    ----------
+    N : complex
+        Mellin moment
+    nf : int
+        Number of active flavors
+    cache: numpy.ndarray
+        Harmonic sum cache
+    variation : tuple
+        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)``
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(3,0)}(N)`
+
+    """
+    gamma_np_p = gamma_nsp(N, nf, cache, variation[3])
+    gamma_qq = gamma_np_p + gamma_ps(N, nf, cache, variation[3])
+    gamma_S = np.array(
+        [
+            [
+                gamma_gg(N, nf, cache, variation[0]),
+                0.0 + 0.0j,
+                gamma_gq(N, nf, cache, variation[1]),
+                0.0 + 0.0j,
+            ],
+            [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+            [gamma_qg(N, nf, cache, variation[2]), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j],
+            [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_np_p],
+        ],
+        np.complex_,
+    )
+    return gamma_S
+
+
+@nb.njit(cache=True)
+def gamma_valence_qed(N, nf, cache, variation):
+    r"""Compute the leading-order valence anomalous dimension matrix for the unified evolution basis.
+
+    .. math::
+        \\gamma_V^{(3,0)} = \\left(\begin{array}{cc}
+        \\gamma_{nsV}^{(3,0)} & 0\\
+        0 & \\gamma_{ns-}^{(3,0)}
+        \\end{array}\right)
+
+    Parameters
+    ----------
+    N : complex
+        Mellin moment
+    nf : int
+        Number of active flavors
+    cache: numpy.ndarray
+        Harmonic sum cache
+    variation : tuple
+        |N3LO| anomalous dimension variation ``(nsm, nsv)``
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{V}^{(3,0)}(N)`
+
+    """
+    gamma_V = np.array(
+        [
+            [gamma_nsv(N, nf, cache, variation[-1]), 0.0],
+            [0.0, gamma_nsm(N, nf, cache, variation[-2])],
+        ],
+        np.complex_,
+    )
+    return gamma_V

tests/eko/kernels/test_kernels_QEDns.py

@@ -77,7 +77,7 @@ def test_zero_true_gamma():
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 n = np.random.rand()
-                gamma_ns = ad.gamma_ns_qed(order, mode, n, nf)
+                gamma_ns = ad.gamma_ns_qed(order, mode, n, nf, (0, 0, 0, 0, 0, 0, 0))
                 for method in methods:
                     for running in running_alpha:
                         np.testing.assert_allclose(
@@ -248,7 +248,9 @@ def test_ode_true_gamma():
                     )
                     a1 = sc.a_s(mu2_to)
                     n = 3 + np.random.rand()
-                    gamma_ns = ad.gamma_ns_qed(order, mode, n, nf)
+                    gamma_ns = ad.gamma_ns_qed(
+                        order, mode, n, nf, (0, 0, 0, 0, 0, 0, 0)
+                    )
                     gammatot = 0.0
                     for i in range(0, order[0] + 1):
                         for j in range(0, order[1] + 1):

tests/eko/kernels/test_kernels_QEDvalence.py

@@ -67,7 +67,7 @@ def test_zero_true_gamma(monkeypatch):
             order = (qcd, qed)
             n = np.random.rand()
             gamma_v = anomalous_dimensions.unpolarized.space_like.gamma_valence_qed(
-                order, n, nf
+                order, n, nf, (0, 0, 0, 0, 0, 0, 0)
             )
             for method in methods:
                 np.testing.assert_allclose(