Approximate massive N3LO for F2 and FL

extras/n3lo_grids/.gitignore

@@ -0,0 +1,2 @@
+*.npy
+

extras/n3lo_grids/yad_grids.py

@@ -0,0 +1,105 @@
+import pathlib
+import sys
+import time
+from multiprocessing import Pool
+
+import adani
+import numpy as np
+
+here = pathlib.Path(__file__).parent / "yadism_grids"
+here.mkdir(exist_ok=True)
+
+nf = int(sys.argv[1])
+n_threads = int(sys.argv[2])
+channel = sys.argv[3]
+order = int(sys.argv[4])
+try:
+    variation = int(sys.argv[5])
+except IndexError:
+    variation = 0
+mufrac = 1.0
+verbose = True
+
+
+def function_to_exe_in_parallel(pair):
+    x, xi = pair
+    m2Q2 = 1 / xi
+    m2mu2 = 1 / xi
+    if order == 3:
+        if channel == "2g":
+            return adani.C2_g3_approximation(
+                x, m2Q2, m2mu2, nf, v=variation, method_flag=1
+            )
+        elif channel == "2q":
+            return adani.C2_ps3_approximation(x, m2Q2, m2mu2, nf, v=variation)
+        elif channel == "Lg":
+            return adani.CL_g3_approximation(
+                x, m2Q2, m2mu2, nf, v=variation, method_flag=1
+            )
+        elif channel == "Lq":
+            return adani.CL_ps3_approximation(x, m2Q2, m2mu2, nf, v=variation)
+        else:
+            raise ValueError("Set channel to one of these: 2g 2q Lg Lq")
+    ## NNLO approximated
+    elif order == 2:
+        if channel == "2g":
+            return adani.C2_g2_approximation(x, m2Q2, m2mu2, v=variation)
+        elif channel == "2q":
+            return adani.C2_ps2_approximation(x, m2Q2, m2mu2, v=variation)
+        elif channel == "Lg":
+            return adani.CL_g2_approximation(x, m2Q2, m2mu2, v=variation)
+        elif channel == "Lq":
+            return adani.CL_ps2_approximation(x, m2Q2, m2mu2, v=variation)
+        else:
+            raise ValueError("Set channel to one of these: 2g 2q Lg Lq")
+    ## NNLO exact
+    elif order == 0:
+        if channel == "2g":
+            return adani.C2_g2(x, m2Q2, m2mu2)
+        elif channel == "2q":
+            return adani.C2_ps2(x, m2Q2, m2mu2)
+        elif channel == "Lg":
+            return adani.CL_g2(x, m2Q2, m2mu2)
+        elif channel == "Lq":
+            return adani.CL_ps2(x, m2Q2, m2mu2)
+        else:
+            raise ValueError("Set channel to one of these: 2g 2q Lg Lq")
+
+
+def run(n_threads, x_grid, xi_grid):
+    grid = []
+    for xi in xi_grid:
+        for x in x_grid:
+            grid.append((x, xi))
+    args = (function_to_exe_in_parallel, grid)
+    with Pool(n_threads) as pool:
+        result = pool.map(*args)
+    return result
+
+
+if __name__ == "__main__":
+    output_file = f"C{channel}_nf{nf}_var{variation}.npy"
+    xfname = here / "x.npy"
+    x_grid = np.load(xfname)
+    xifname = here / "xi.npy"
+    xi_grid = np.load(xifname)
+
+    if verbose:
+        print(
+            f"Computation of the grid for the coefficient function C{channel} for nf = {nf}, and µ/Q = {mufrac}, variation = {variation}"
+        )
+        print(f"Size of the grid (x,xi) = ({len(x_grid)},{len(xi_grid)})")
+        print(
+            "This may take a while (depending on the number of threads you choose). In order to spend this time, I would suggest you this interesting view:"
+        )
+        print("https://www.youtube.com/watch?v=53pG68KCUMI")
+
+    start = time.perf_counter()
+    res_vec = np.array(run(n_threads, x_grid, xi_grid))
+    if verbose:
+        print("total running time: ", time.perf_counter() - start)
+
+    res_mat = res_vec.reshape(len(xi_grid), len(x_grid))
+    if verbose:
+        print("Saving grid in ", output_file)
+    np.save(here / output_file, res_mat)

pyproject.toml

@@ -43,11 +43,12 @@ numba = "^0.57.0"
 eko = "^0.13.4"
 pandas = "^1.3.0"
 rich = "^12.4.4"
+adani = "^0.13"
 LeProHQ = "^0.2.5"
 # pineappl dependency
 pineappl = { version = "^0.6.0", optional = true }
 # yadmark
-banana-hep = { version = "^0.6.6", optional = true }
+banana-hep = { version = "^0.6.11", optional = true }
 sqlalchemy = { version = "^1.4.21", optional = true }
 
 [tool.poetry.group.dev.dependencies]
@@ -161,6 +162,7 @@ extension-pkg-whitelist = [
   "scipy.special",
   "pineappl",
   "apfelpy",
+  "adani"
 ]
 ignore-paths = ["benchmarks/", "docs/", "tests/"]
 jobs = 1 # has to be 1 as pylint is NOT threadsafe

src/yadism/coefficient_functions/asy/f2_nc.py

@@ -1,3 +1,5 @@
+import adani
+
 from ..light import f2_nc as light
 from ..partonic_channel import RSL
 from . import partonic_channel as pc
@@ -19,6 +21,14 @@ def cg_LL_NNLO(z, args):
 
         return RSL(cg_LL_NNLO, args=[self.L])
 
+    def N3LO(self):
+        def cg_LL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.C2_g3_highscale_LL(z, nf) * L**3
+
+        return RSL(cg_LL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNLLGluon(pc.NeutralCurrentBaseAsy):
     def NLO(self):
@@ -34,6 +44,14 @@ def cg_NLL_NNLO(z, args):
 
         return RSL(cg_NLL_NNLO, args=[self.L])
 
+    def N3LO(self):
+        def cg_NLL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.C2_g3_highscale_NLL(z, nf) * L**2
+
+        return RSL(cg_NLL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNNLLGluon(pc.NeutralCurrentBaseAsy):
     def NNLO(self):
@@ -42,9 +60,23 @@ def cg_NNLL_NNLO(z, _args):
 
         return RSL(cg_NNLL_NNLO)
 
+    def N3LO(self):
+        def cg_NNLL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.C2_g3_highscale_N2LL(z, nf) * L
+
+        return RSL(cg_NNLL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNNNLLGluon(pc.NeutralCurrentBaseAsy):
-    pass
+    def N3LO(self):
+        def cg_NNNLL_N3LO(z, args):
+            nf = int(args[0])
+            variation = int(args[1])
+            return adani.C2_g3_highscale_N3LL(z, nf, variation)
+
+        return RSL(cg_NNNLL_N3LO, args=[self.nf, self.n3lo_cf_variation])
 
 
 class AsyLLSinglet(pc.NeutralCurrentBaseAsy):
@@ -55,6 +87,14 @@ def cps_LL_NNLO(z, args):
 
         return RSL(cps_LL_NNLO, args=[self.L])
 
+    def N3LO(self):
+        def cps_LL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.C2_ps3_highscale_LL(z, nf) * L**3
+
+        return RSL(cps_LL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNLLSinglet(pc.NeutralCurrentBaseAsy):
     def NNLO(self):
@@ -64,6 +104,14 @@ def cps_NLL_NNLO(z, args):
 
         return RSL(cps_NLL_NNLO, args=[self.L])
 
+    def N3LO(self):
+        def cps_NLL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.C2_ps3_highscale_NLL(z, nf) * L**2
+
+        return RSL(cps_NLL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNNLLSinglet(pc.NeutralCurrentBaseAsy):
     def NNLO(self):
@@ -72,9 +120,22 @@ def cps_NNLL_NNLO(z, _args):
 
         return RSL(cps_NNLL_NNLO)
 
+    def N3LO(self):
+        def cps_NNLL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.C2_ps3_highscale_N2LL(z, nf) * L
+
+        return RSL(cps_NNLL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNNNLLSinglet(pc.NeutralCurrentBaseAsy):
-    pass
+    def N3LO(self):
+        def cps_NNNLL_N3LO(z, args):
+            nf = int(args[0])
+            return adani.C2_ps3_highscale_N3LL(z, nf)
+
+        return RSL(cps_NNNLL_N3LO, args=[self.nf])
 
 
 class AsyLLNonSinglet(pc.NeutralCurrentBaseAsy):

src/yadism/coefficient_functions/asy/fl_nc.py

@@ -1,3 +1,5 @@
+import adani
+
 from ..light import fl_nc as light
 from ..partonic_channel import RSL, EmptyPartonicChannel
 from . import partonic_channel as pc
@@ -22,6 +24,13 @@ def cg_NLL_NNLO(z, args):
 
         return RSL(cg_NLL_NNLO, args=[self.L])
 
+    def N3LO(self):
+        def cg_NLL_N3LO(z, args):
+            L = -args[0]
+            return adani.CL_g3_highscale_NLL(z) * L**2
+
+        return RSL(cg_NLL_N3LO, args=[self.L])
+
 
 class AsyNNLLGluon(pc.NeutralCurrentBaseAsy):
     def NNLO(self):
@@ -30,9 +39,22 @@ def cg_NNLL_NNLO(z, _args):
 
         return RSL(cg_NNLL_NNLO)
 
+    def N3LO(self):
+        def cg_NNLL_N3LO(z, args):
+            L = -args[0]
+            nf = int(args[1])
+            return adani.CL_g3_highscale_N2LL(z, nf) * L
+
+        return RSL(cg_NNLL_N3LO, args=[self.L, self.nf])
+
 
 class AsyNNNLLGluon(pc.NeutralCurrentBaseAsy):
-    pass
+    def N3LO(self):
+        def cg_NNNLL_N3LO(z, args):
+            nf = int(args[0])
+            return adani.CL_g3_highscale_N3LL(z, nf)
+
+        return RSL(cg_NNNLL_N3LO, args=[self.nf])
 
 
 class AsyLLSinglet(EmptyPartonicChannel):
@@ -47,6 +69,13 @@ def cps_NLL_NNLO(z, args):
 
         return RSL(cps_NLL_NNLO, args=[self.L])
 
+    def N3LO(self):
+        def cps_NLL_N3LO(z, args):
+            L = -args[0]
+            return adani.CL_ps3_highscale_NLL(z) * L**2
+
+        return RSL(cps_NLL_N3LO, args=[self.L])
+
 
 class AsyNNLLSinglet(pc.NeutralCurrentBaseAsy):
     def NNLO(self):
@@ -55,9 +84,21 @@ def cps_NNLL_NNLO(z, _args):
 
         return RSL(cps_NNLL_NNLO)
 
+    def N3LO(self):
+        def cps_NNLL_N3LO(z, args):
+            L = -args[0]
+            return adani.CL_ps3_highscale_N2LL(z) * L
+
+        return RSL(cps_NNLL_N3LO, args=[self.L])
+
 
 class AsyNNNLLSinglet(pc.NeutralCurrentBaseAsy):
-    pass
+    def N3LO(self):
+        def cps_NNNLL_N3LO(z, args):
+            nf = int(args[0])
+            return adani.CL_ps3_highscale_N3LL(z, nf)
+
+        return RSL(cps_NNNLL_N3LO, args=[self.nf])
 
 
 class AsyLLNonSinglet(EmptyPartonicChannel):

src/yadism/coefficient_functions/asy/kernels.py

@@ -106,14 +106,20 @@ def generate_heavy_asy(esf, nf, pto_evol, ihq):
             is_pv,
         )
         if not is_pv:
+            n3lo_cf_variation = esf.info.theory["n3lo_cf_variation"]
             for c, channel in (("g", "Gluon"), ("s", "Singlet")):
                 for res in range(pto_evol + 1):
                     name = "Asy" + ("N" * res) + "LL" + channel
                     for av in ("AA", "VV"):
                         asys.append(
                             kernels.Kernel(
                                 asy_weights[f"{c}{av}"],
-                                asy_cfs.__getattribute__(name)(esf, nf, m2hq=m2hq),
+                                asy_cfs.__getattribute__(name)(
+                                    esf,
+                                    nf,
+                                    m2hq=m2hq,
+                                    n3lo_cf_variation=n3lo_cf_variation,
+                                ),
                             )
                         )
     return asys

src/yadism/coefficient_functions/asy/partonic_channel.py

@@ -11,9 +11,10 @@
 class PartonicChannelAsy(pc.PartonicChannel):
     """Massless limit of a coeficient function."""
 
-    def __init__(self, *args, m2hq):
+    def __init__(self, *args, m2hq, n3lo_cf_variation=0):
         super().__init__(*args)
         self.L = np.log(self.ESF.Q2 / m2hq)
+        self.n3lo_cf_variation = n3lo_cf_variation
 
 
 class PartonicChannelAsyLLIntrinsic(PartonicChannelAsy):