From 5bb03896653b8db360f9454f19afbfecadf0e50a Mon Sep 17 00:00:00 2001
From: Felix Hekhorn <felix.hekhorn@uni-tuebingen.de>
Date: Wed, 9 Aug 2023 18:42:13 +0200
Subject: [PATCH 01/41] Init LH23 benchmark

---
 extras/lh_bench_23/.gitignore |   1 +
 extras/lh_bench_23/run.py     | 151 ++++++++++++++++++++++++++++++++++
 2 files changed, 152 insertions(+)
 create mode 100644 extras/lh_bench_23/.gitignore
 create mode 100644 extras/lh_bench_23/run.py

diff --git a/extras/lh_bench_23/.gitignore b/extras/lh_bench_23/.gitignore
new file mode 100644
index 000000000..d874ad67c
--- /dev/null
+++ b/extras/lh_bench_23/.gitignore
@@ -0,0 +1 @@
+*.tar
diff --git a/extras/lh_bench_23/run.py b/extras/lh_bench_23/run.py
new file mode 100644
index 000000000..eded3c7e5
--- /dev/null
+++ b/extras/lh_bench_23/run.py
@@ -0,0 +1,151 @@
+import argparse
+import copy
+import pathlib
+from math import inf, nan
+
+import numpy as np
+import pandas as pd
+import yaml
+from banana import toy
+
+import eko
+from eko import basis_rotation as br
+from eko.interpolation import lambertgrid
+from eko.io import runcards
+from eko.io.types import ReferenceRunning
+from eko.runner.managed import solve
+from ekobox import apply
+from ekomark.benchmark.external.LHA_utils import here as there
+
+_sqrt2 = float(np.sqrt(2))
+
+# VFNS theory settings
+_t_vfns = dict(
+    order=[3, 0],
+    couplings=dict(
+        alphas=0.35,
+        alphaem=0.007496,
+        scale=_sqrt2,
+        num_flavs_ref=3,
+        max_num_flavs=6,
+    ),
+    heavy=dict(
+        num_flavs_init=3,
+        num_flavs_max_pdf=6,
+        intrinsic_flavors=[],
+        masses=[ReferenceRunning([mq, nan]) for mq in (_sqrt2, 4.5, 175.0)],
+        masses_scheme="POLE",
+        matching_ratios=[1.0, 1.0, 1.0],
+    ),
+    xif=1.0,
+    n3lo_ad_variation=(0, 0, 0, 0),
+)
+t_vfns = runcards.TheoryCard.from_dict(_t_vfns)
+
+# FFNS theory settings
+_t_ffns = copy.deepcopy(_t_vfns)
+_t_ffns["couplings"]["num_flavs_ref"] = 4
+_t_ffns["heavy"]["num_flavs_init"] = 4
+_t_ffns["heavy"]["masses"] = [
+    ReferenceRunning([0, nan]),
+    ReferenceRunning([inf, nan]),
+    ReferenceRunning([inf, nan]),
+]
+t_ffns = runcards.TheoryCard.from_dict(_t_ffns)
+
+# VFNS operator settings
+_o_vfns = dict(
+    mu0=_sqrt2,
+    mugrid=[(100.0, 5)],
+    xgrid=lambertgrid(60).tolist(),
+    configs=dict(
+        evolution_method="iterate-exact",
+        ev_op_max_order=[10, 0],
+        ev_op_iterations=30,
+        interpolation_polynomial_degree=4,
+        interpolation_is_log=True,
+        scvar_method="exponentiated",
+        inversion_method=None,
+        n_integration_cores=1,
+        polarized=False,
+        time_like=False,
+    ),
+    debug=dict(
+        skip_singlet=False,
+        skip_non_singlet=False,
+    ),
+)
+o_vfns = runcards.OperatorCard.from_dict(_o_vfns)
+
+# FFNS operator settings
+_o_ffns = copy.deepcopy(_o_vfns)
+_o_ffns["mugrid"] = [(100.0, 4)]
+o_ffns = runcards.OperatorCard.from_dict(_o_ffns)
+
+# setup flavor rotations
+labels = ["u_v", "d_v", "L_m", "L_p", "s_v", "s_p", "c_p", "g"]
+rotate_to_LHA = np.zeros((len(labels), 14))
+# u_v = u - ubar
+rotate_to_LHA[0][br.flavor_basis_pids.index(-2)] = -1
+rotate_to_LHA[0][br.flavor_basis_pids.index(2)] = 1
+# d_v = d - dbar
+rotate_to_LHA[1][br.flavor_basis_pids.index(-1)] = -1
+rotate_to_LHA[1][br.flavor_basis_pids.index(1)] = 1
+# L_- = dbar - ubar
+rotate_to_LHA[2][br.flavor_basis_pids.index(-1)] = 1
+rotate_to_LHA[2][br.flavor_basis_pids.index(-2)] = -1
+# 2L_+ = 2dbar + 2ubar
+rotate_to_LHA[3][br.flavor_basis_pids.index(-1)] = 2
+rotate_to_LHA[3][br.flavor_basis_pids.index(-2)] = 2
+# s_v = s - sbar
+rotate_to_LHA[4][br.flavor_basis_pids.index(-3)] = -1
+rotate_to_LHA[4][br.flavor_basis_pids.index(3)] = 1
+# s_+ = s + sbar
+rotate_to_LHA[5][br.flavor_basis_pids.index(-3)] = 1
+rotate_to_LHA[5][br.flavor_basis_pids.index(3)] = 1
+# c_+ = c + cbar
+rotate_to_LHA[6][br.flavor_basis_pids.index(-4)] = 1
+rotate_to_LHA[6][br.flavor_basis_pids.index(4)] = 1
+# g = g
+rotate_to_LHA[7][br.flavor_basis_pids.index(21)] = 1
+
+# setup x rotation
+xgrid = np.array([1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.1, 0.3, 0.5, 0.7, 0.9])
+
+# eko path
+p = pathlib.Path("FFNS.tar")
+
+# reference values
+with open(there / "LHA.yaml", encoding="utf-8") as o:
+    ref_data = yaml.safe_load(o)
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--rerun", help="Rerun eko", action="store_true")
+    args = parser.parse_args()
+
+    # recompute?
+    if args.rerun:
+        print("Rerunning eko ...")
+        p.unlink(True)
+        solve(t_ffns, o_ffns, p)
+
+    # apply PDF
+    out = {}
+    with eko.EKO.read(p) as eko_:
+        pdf = apply.apply_pdf_flavor(
+            eko_, toy.mkPDF("ToyLH", 0), xgrid, rotate_to_LHA, labels
+        )
+        for lab, f in pdf[(10000.0, 4)]["pdfs"].items():
+            out[lab] = xgrid * f
+
+    # display result
+    pd.set_option("display.float_format", "{:.4e}".format)
+    me = pd.DataFrame(out)
+    print("EKO")
+    print(me)
+
+    # load reference
+    ref = pd.DataFrame(ref_data["table14"]["part1"])
+    print("rel. distance to reference")
+    print((me - ref) / ref)

From 36049b323e63e86da8b4d50b919b9e6f61932222 Mon Sep 17 00:00:00 2001
From: Felix Hekhorn <felix.hekhorn@uni-tuebingen.de>
Date: Wed, 9 Aug 2023 18:43:29 +0200
Subject: [PATCH 02/41] Pass labels in apply_pdf

---
 src/ekobox/apply.py             | 20 ++++++++++----------
 src/ekomark/benchmark/runner.py |  5 ++++-
 2 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/src/ekobox/apply.py b/src/ekobox/apply.py
index c392ffe78..354461c0a 100644
--- a/src/ekobox/apply.py
+++ b/src/ekobox/apply.py
@@ -37,10 +37,12 @@ def apply_pdf(
     if rotate_to_evolution_basis:
         if not qed:
             rotate_flavor_to_evolution = br.rotate_flavor_to_evolution
+            labels = br.evol_basis
         else:
             rotate_flavor_to_evolution = br.rotate_flavor_to_unified_evolution
+            labels = br.unified_evol_basis
         return apply_pdf_flavor(
-            eko, lhapdf_like, targetgrid, rotate_flavor_to_evolution, qed
+            eko, lhapdf_like, targetgrid, rotate_flavor_to_evolution, labels=labels
         )
     return apply_pdf_flavor(eko, lhapdf_like, targetgrid)
 
@@ -49,7 +51,7 @@ def apply_pdf(
 
 
 def apply_pdf_flavor(
-    eko: EKO, lhapdf_like, targetgrid=None, flavor_rotation=None, qed=False
+    eko: EKO, lhapdf_like, targetgrid=None, flavor_rotation=None, labels=None
 ):
     """
     Apply all available operators to the input PDFs.
@@ -65,8 +67,8 @@ def apply_pdf_flavor(
             if given, interpolates to the pdfs given at targetgrid (instead of xgrid)
         flavor_rotation : np.ndarray
             Rotation matrix in flavor space
-        qed : bool
-            activate qed
+        labels : list
+            list of labels
 
     Returns
     -------
@@ -103,16 +105,14 @@ def apply_pdf_flavor(
             pdf = flavor_rotation @ np.array(
                 [op["pdfs"][pid] for pid in br.flavor_basis_pids]
             )
-            if not qed:
-                evol_basis = br.evol_basis
-            else:
-                evol_basis = br.unified_evol_basis
-            op["pdfs"] = dict(zip(evol_basis, pdf))
+            if labels is None:
+                labels = list(range(flavor_rotation.shape[0]))
+            op["pdfs"] = dict(zip(labels, pdf))
             if op["errors"] is not None:
                 errors = flavor_rotation @ np.array(
                     [op["errors"][pid] for pid in br.flavor_basis_pids]
                 )
-                op["errors"] = dict(zip(evol_basis, errors))
+                op["errors"] = dict(zip(labels, errors))
 
     # rotate/interpolate to target grid
     if targetgrid is not None:
diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py
index 8c07e7153..180f439e6 100644
--- a/src/ekomark/benchmark/runner.py
+++ b/src/ekomark/benchmark/runner.py
@@ -224,12 +224,15 @@ def log(self, theory, _, pdf, me, ext):
         # LHA NNLO VFNS needs a special treatment
         # Valence contains only u and d
         rotate_to_evolution = None
+        labels = None
         qed = theory["QED"] > 0
         if self.rotate_to_evolution_basis:
             if not qed:
                 rotate_to_evolution = br.rotate_flavor_to_evolution.copy()
+                labels = br.evol_basis
             else:
                 rotate_to_evolution = br.rotate_flavor_to_unified_evolution.copy()
+                labels = br.unified_evol_basis
             if (
                 self.external == "LHA"
                 and theory["PTO"] == 2
@@ -243,7 +246,7 @@ def log(self, theory, _, pdf, me, ext):
                 pdf,
                 xgrid,
                 flavor_rotation=rotate_to_evolution,
-                qed=qed,
+                labels=labels,
             )
         for q2, ref_pdfs in ext["values"].items():
             log_tab = dfdict.DFdict()

From 031527e5015a4b5247d8f87b74bfdd176bd95585 Mon Sep 17 00:00:00 2001
From: Felix Hekhorn <felix.hekhorn@uni-tuebingen.de>
Date: Thu, 17 Aug 2023 12:37:55 +0200
Subject: [PATCH 03/41] Split and enhance lh bench script

---
 extras/lh_bench_23/cards.py |  87 ++++++++++++++++++++++++
 extras/lh_bench_23/run.py   | 129 +++++++++++++++---------------------
 2 files changed, 142 insertions(+), 74 deletions(-)
 create mode 100644 extras/lh_bench_23/cards.py

diff --git a/extras/lh_bench_23/cards.py b/extras/lh_bench_23/cards.py
new file mode 100644
index 000000000..a4c488148
--- /dev/null
+++ b/extras/lh_bench_23/cards.py
@@ -0,0 +1,87 @@
+import copy
+from math import inf, nan
+
+import numpy as np
+
+from eko.interpolation import lambertgrid
+from eko.io import runcards
+from eko.io.types import ReferenceRunning
+
+_sqrt2 = float(np.sqrt(2))
+
+# VFNS theory settings
+_t_vfns = dict(
+    order=[3, 0],
+    couplings=dict(
+        alphas=0.35,
+        alphaem=0.007496,
+        scale=_sqrt2,
+        num_flavs_ref=3,
+        max_num_flavs=6,
+    ),
+    heavy=dict(
+        num_flavs_init=3,
+        num_flavs_max_pdf=6,
+        intrinsic_flavors=[],
+        masses=[ReferenceRunning([mq, nan]) for mq in (_sqrt2, 4.5, 175.0)],
+        masses_scheme="POLE",
+        matching_ratios=[1.0, 1.0, 1.0],
+    ),
+    xif=1.0,
+    n3lo_ad_variation=(0, 0, 0, 0),
+)
+
+
+def vfns_theory(xif=1.0):
+    """Generate a VFNS theory card."""
+    tt = copy.deepcopy(_t_vfns)
+    tt["xif"] = xif
+    return runcards.TheoryCard.from_dict(tt)
+
+
+# FFNS theory settings
+_t_ffns = copy.deepcopy(_t_vfns)
+_t_ffns["couplings"]["num_flavs_ref"] = 4
+_t_ffns["heavy"]["num_flavs_init"] = 4
+_t_ffns["heavy"]["masses"] = [
+    ReferenceRunning([0, nan]),
+    ReferenceRunning([inf, nan]),
+    ReferenceRunning([inf, nan]),
+]
+
+
+def ffns_theory(xif=1.0):
+    """Generate a VFNS theory card."""
+    tt = copy.deepcopy(_t_ffns)
+    tt["xif"] = xif
+    return runcards.TheoryCard.from_dict(tt)
+
+
+# VFNS operator settings
+_o_vfns = dict(
+    mu0=_sqrt2,
+    mugrid=[(100.0, 5)],
+    xgrid=lambertgrid(60).tolist(),
+    configs=dict(
+        evolution_method="iterate-exact",
+        ev_op_max_order=[10, 0],
+        ev_op_iterations=30,
+        interpolation_polynomial_degree=4,
+        interpolation_is_log=True,
+        scvar_method="exponentiated",
+        inversion_method=None,
+        n_integration_cores=1,
+        polarized=False,
+        time_like=False,
+    ),
+    debug=dict(
+        skip_singlet=False,
+        skip_non_singlet=False,
+    ),
+)
+vfns_operator = runcards.OperatorCard.from_dict(_o_vfns)
+
+# FFNS operator settings
+_o_ffns = copy.deepcopy(_o_vfns)
+_o_ffns["mugrid"] = [(100.0, 4)]
+ffns_operator = runcards.OperatorCard.from_dict(_o_ffns)
diff --git a/extras/lh_bench_23/run.py b/extras/lh_bench_23/run.py
index eded3c7e5..b9146c1f1 100644
--- a/extras/lh_bench_23/run.py
+++ b/extras/lh_bench_23/run.py
@@ -1,87 +1,22 @@
 import argparse
-import copy
+import logging
 import pathlib
-from math import inf, nan
+import sys
 
 import numpy as np
 import pandas as pd
 import yaml
 from banana import toy
+from cards import ffns_operator, ffns_theory, vfns_operator, vfns_theory
 
 import eko
 from eko import basis_rotation as br
-from eko.interpolation import lambertgrid
-from eko.io import runcards
-from eko.io.types import ReferenceRunning
 from eko.runner.managed import solve
 from ekobox import apply
 from ekomark.benchmark.external.LHA_utils import here as there
 
 _sqrt2 = float(np.sqrt(2))
 
-# VFNS theory settings
-_t_vfns = dict(
-    order=[3, 0],
-    couplings=dict(
-        alphas=0.35,
-        alphaem=0.007496,
-        scale=_sqrt2,
-        num_flavs_ref=3,
-        max_num_flavs=6,
-    ),
-    heavy=dict(
-        num_flavs_init=3,
-        num_flavs_max_pdf=6,
-        intrinsic_flavors=[],
-        masses=[ReferenceRunning([mq, nan]) for mq in (_sqrt2, 4.5, 175.0)],
-        masses_scheme="POLE",
-        matching_ratios=[1.0, 1.0, 1.0],
-    ),
-    xif=1.0,
-    n3lo_ad_variation=(0, 0, 0, 0),
-)
-t_vfns = runcards.TheoryCard.from_dict(_t_vfns)
-
-# FFNS theory settings
-_t_ffns = copy.deepcopy(_t_vfns)
-_t_ffns["couplings"]["num_flavs_ref"] = 4
-_t_ffns["heavy"]["num_flavs_init"] = 4
-_t_ffns["heavy"]["masses"] = [
-    ReferenceRunning([0, nan]),
-    ReferenceRunning([inf, nan]),
-    ReferenceRunning([inf, nan]),
-]
-t_ffns = runcards.TheoryCard.from_dict(_t_ffns)
-
-# VFNS operator settings
-_o_vfns = dict(
-    mu0=_sqrt2,
-    mugrid=[(100.0, 5)],
-    xgrid=lambertgrid(60).tolist(),
-    configs=dict(
-        evolution_method="iterate-exact",
-        ev_op_max_order=[10, 0],
-        ev_op_iterations=30,
-        interpolation_polynomial_degree=4,
-        interpolation_is_log=True,
-        scvar_method="exponentiated",
-        inversion_method=None,
-        n_integration_cores=1,
-        polarized=False,
-        time_like=False,
-    ),
-    debug=dict(
-        skip_singlet=False,
-        skip_non_singlet=False,
-    ),
-)
-o_vfns = runcards.OperatorCard.from_dict(_o_vfns)
-
-# FFNS operator settings
-_o_ffns = copy.deepcopy(_o_vfns)
-_o_ffns["mugrid"] = [(100.0, 4)]
-o_ffns = runcards.OperatorCard.from_dict(_o_ffns)
-
 # setup flavor rotations
 labels = ["u_v", "d_v", "L_m", "L_p", "s_v", "s_p", "c_p", "g"]
 rotate_to_LHA = np.zeros((len(labels), 14))
@@ -112,8 +47,6 @@
 # setup x rotation
 xgrid = np.array([1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.1, 0.3, 0.5, 0.7, 0.9])
 
-# eko path
-p = pathlib.Path("FFNS.tar")
 
 # reference values
 with open(there / "LHA.yaml", encoding="utf-8") as o:
@@ -121,14 +54,61 @@
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
+    parser.add_argument("scheme", help="FFNS or VFNS?")
+    parser.add_argument("sv", help="scale variation: up, central, or down")
     parser.add_argument("--rerun", help="Rerun eko", action="store_true")
+    parser.add_argument(
+        "-v", "--verbose", help="Print eko log to screen", action="store_true"
+    )
     args = parser.parse_args()
 
+    # determine xif
+    if "central".startswith(args.sv):
+        xif = 1.0
+        sv = "central"
+        part = 1
+    elif "up".startswith(args.sv):
+        xif = _sqrt2
+        sv = "up"
+        part = 2
+    elif "down".startswith(args.sv):
+        xif = 1.0 / _sqrt2
+        sv = "down"
+        part = 3
+    else:
+        raise ValueError(
+            "sv has to be up, central, or down - or any abbreviation there of"
+        )
+
+    # determine scheme
+    if args.scheme == "FFNS":
+        scheme = "FFNS"
+        t = ffns_theory(xif)
+        o = ffns_operator
+        tab = 14
+    elif args.scheme == "VFNS":
+        scheme = "VFNS"
+        t = vfns_theory(xif)
+        o = vfns_operator
+        tab = 15
+    else:
+        raise ValueError("scheme has to be FFNS or VFNS")
+
+    # eko path
+    p = pathlib.Path(f"{scheme}-{sv}.tar")
+
     # recompute?
-    if args.rerun:
-        print("Rerunning eko ...")
+    if not p.exists() or args.rerun:
+        print("(Re)running eko ...")
         p.unlink(True)
-        solve(t_ffns, o_ffns, p)
+        if args.verbose:
+            logStdout = logging.StreamHandler(sys.stdout)
+            logStdout.setLevel(logging.INFO)
+            logStdout.setFormatter(logging.Formatter("%(message)s"))
+            logging.getLogger("eko").handlers = []
+            logging.getLogger("eko").addHandler(logStdout)
+            logging.getLogger("eko").setLevel(logging.INFO)
+        solve(t, o, p)
 
     # apply PDF
     out = {}
@@ -146,6 +126,7 @@
     print(me)
 
     # load reference
-    ref = pd.DataFrame(ref_data["table14"]["part1"])
+    ref = pd.DataFrame(ref_data[f"table{tab}"][f"part{part}"])
+    print()
     print("rel. distance to reference")
     print((me - ref) / ref)

From 403afea5d522219d56335eda6fd43e7f7726a0a5 Mon Sep 17 00:00:00 2001
From: Felix Hekhorn <felix.hekhorn@uni-tuebingen.de>
Date: Thu, 17 Aug 2023 13:13:27 +0200
Subject: [PATCH 04/41] Fix and improve LH bench script

---
 extras/lh_bench_23/.gitignore              |   1 +
 extras/lh_bench_23/cards.py                |  87 -------------
 extras/lh_bench_23/cfg.py                  | 144 +++++++++++++++++++++
 extras/lh_bench_23/{run.py => run-nnlo.py} |  50 +++----
 4 files changed, 163 insertions(+), 119 deletions(-)
 delete mode 100644 extras/lh_bench_23/cards.py
 create mode 100644 extras/lh_bench_23/cfg.py
 rename extras/lh_bench_23/{run.py => run-nnlo.py} (66%)

diff --git a/extras/lh_bench_23/.gitignore b/extras/lh_bench_23/.gitignore
index d874ad67c..11fa526db 100644
--- a/extras/lh_bench_23/.gitignore
+++ b/extras/lh_bench_23/.gitignore
@@ -1 +1,2 @@
 *.tar
+*.csv
diff --git a/extras/lh_bench_23/cards.py b/extras/lh_bench_23/cards.py
deleted file mode 100644
index a4c488148..000000000
--- a/extras/lh_bench_23/cards.py
+++ /dev/null
@@ -1,87 +0,0 @@
-import copy
-from math import inf, nan
-
-import numpy as np
-
-from eko.interpolation import lambertgrid
-from eko.io import runcards
-from eko.io.types import ReferenceRunning
-
-_sqrt2 = float(np.sqrt(2))
-
-# VFNS theory settings
-_t_vfns = dict(
-    order=[3, 0],
-    couplings=dict(
-        alphas=0.35,
-        alphaem=0.007496,
-        scale=_sqrt2,
-        num_flavs_ref=3,
-        max_num_flavs=6,
-    ),
-    heavy=dict(
-        num_flavs_init=3,
-        num_flavs_max_pdf=6,
-        intrinsic_flavors=[],
-        masses=[ReferenceRunning([mq, nan]) for mq in (_sqrt2, 4.5, 175.0)],
-        masses_scheme="POLE",
-        matching_ratios=[1.0, 1.0, 1.0],
-    ),
-    xif=1.0,
-    n3lo_ad_variation=(0, 0, 0, 0),
-)
-
-
-def vfns_theory(xif=1.0):
-    """Generate a VFNS theory card."""
-    tt = copy.deepcopy(_t_vfns)
-    tt["xif"] = xif
-    return runcards.TheoryCard.from_dict(tt)
-
-
-# FFNS theory settings
-_t_ffns = copy.deepcopy(_t_vfns)
-_t_ffns["couplings"]["num_flavs_ref"] = 4
-_t_ffns["heavy"]["num_flavs_init"] = 4
-_t_ffns["heavy"]["masses"] = [
-    ReferenceRunning([0, nan]),
-    ReferenceRunning([inf, nan]),
-    ReferenceRunning([inf, nan]),
-]
-
-
-def ffns_theory(xif=1.0):
-    """Generate a VFNS theory card."""
-    tt = copy.deepcopy(_t_ffns)
-    tt["xif"] = xif
-    return runcards.TheoryCard.from_dict(tt)
-
-
-# VFNS operator settings
-_o_vfns = dict(
-    mu0=_sqrt2,
-    mugrid=[(100.0, 5)],
-    xgrid=lambertgrid(60).tolist(),
-    configs=dict(
-        evolution_method="iterate-exact",
-        ev_op_max_order=[10, 0],
-        ev_op_iterations=30,
-        interpolation_polynomial_degree=4,
-        interpolation_is_log=True,
-        scvar_method="exponentiated",
-        inversion_method=None,
-        n_integration_cores=1,
-        polarized=False,
-        time_like=False,
-    ),
-    debug=dict(
-        skip_singlet=False,
-        skip_non_singlet=False,
-    ),
-)
-vfns_operator = runcards.OperatorCard.from_dict(_o_vfns)
-
-# FFNS operator settings
-_o_ffns = copy.deepcopy(_o_vfns)
-_o_ffns["mugrid"] = [(100.0, 4)]
-ffns_operator = runcards.OperatorCard.from_dict(_o_ffns)
diff --git a/extras/lh_bench_23/cfg.py b/extras/lh_bench_23/cfg.py
new file mode 100644
index 000000000..ebafc1a13
--- /dev/null
+++ b/extras/lh_bench_23/cfg.py
@@ -0,0 +1,144 @@
+import copy
+from math import inf, nan
+
+import numpy as np
+
+from eko import basis_rotation as br
+from eko.interpolation import lambertgrid
+from eko.io import runcards
+from eko.io.types import ReferenceRunning
+
+_sqrt2 = float(np.sqrt(2))
+
+# theory settings
+# ---------------
+_t_vfns = dict(
+    order=[3, 0],
+    couplings=dict(
+        alphas=0.35,
+        alphaem=0.007496,
+        scale=_sqrt2,
+        num_flavs_ref=3,
+        max_num_flavs=6,
+    ),
+    heavy=dict(
+        num_flavs_init=3,
+        num_flavs_max_pdf=6,
+        intrinsic_flavors=[],
+        masses=[ReferenceRunning([mq, nan]) for mq in (_sqrt2, 4.5, 175.0)],
+        masses_scheme="POLE",
+        matching_ratios=[1.0, 1.0, 1.0],
+    ),
+    xif=1.0,
+    n3lo_ad_variation=(0, 0, 0, 0),
+)
+
+
+def vfns_theory(xif=1.0):
+    """Generate a VFNS theory card."""
+    tt = copy.deepcopy(_t_vfns)
+    tt["xif"] = xif
+    return runcards.TheoryCard.from_dict(tt)
+
+
+_t_ffns = copy.deepcopy(_t_vfns)
+_t_ffns["couplings"]["num_flavs_ref"] = 4
+_t_ffns["heavy"]["num_flavs_init"] = 4
+_t_ffns["heavy"]["masses"] = [
+    ReferenceRunning([0, nan]),
+    ReferenceRunning([inf, nan]),
+    ReferenceRunning([inf, nan]),
+]
+
+
+def ffns_theory(xif=1.0):
+    """Generate a VFNS theory card."""
+    tt = copy.deepcopy(_t_ffns)
+    tt["xif"] = xif
+    return runcards.TheoryCard.from_dict(tt)
+
+
+# operator settings
+# -----------------
+_o_vfns = dict(
+    mu0=_sqrt2,
+    mugrid=[(100.0, 5)],
+    xgrid=lambertgrid(60).tolist(),
+    configs=dict(
+        evolution_method="iterate-exact",
+        ev_op_max_order=[10, 0],
+        ev_op_iterations=30,
+        interpolation_polynomial_degree=4,
+        interpolation_is_log=True,
+        scvar_method="exponentiated",
+        inversion_method=None,
+        n_integration_cores=-2,
+        polarized=False,
+        time_like=False,
+    ),
+    debug=dict(
+        skip_singlet=False,
+        skip_non_singlet=False,
+    ),
+)
+vfns_operator = runcards.OperatorCard.from_dict(_o_vfns)
+
+_o_ffns = copy.deepcopy(_o_vfns)
+_o_ffns["mugrid"] = [(100.0, 4)]
+ffns_operator = runcards.OperatorCard.from_dict(_o_ffns)
+
+
+# flavor rotations
+# ----------------
+
+ffns_labels = ["u_v", "d_v", "L_m", "L_p", "s_v", "s_p", "c_p", "g"]
+ffns_rotate_to_LHA = np.zeros((len(ffns_labels), 14))
+# u_v = u - ubar
+ffns_rotate_to_LHA[0][br.flavor_basis_pids.index(-2)] = -1
+ffns_rotate_to_LHA[0][br.flavor_basis_pids.index(2)] = 1
+# d_v = d - dbar
+ffns_rotate_to_LHA[1][br.flavor_basis_pids.index(-1)] = -1
+ffns_rotate_to_LHA[1][br.flavor_basis_pids.index(1)] = 1
+# L_- = dbar - ubar
+ffns_rotate_to_LHA[2][br.flavor_basis_pids.index(-1)] = 1
+ffns_rotate_to_LHA[2][br.flavor_basis_pids.index(-2)] = -1
+# 2L_+ = 2dbar + 2ubar
+ffns_rotate_to_LHA[3][br.flavor_basis_pids.index(-1)] = 2
+ffns_rotate_to_LHA[3][br.flavor_basis_pids.index(-2)] = 2
+# s_v = s - sbar
+ffns_rotate_to_LHA[4][br.flavor_basis_pids.index(-3)] = -1
+ffns_rotate_to_LHA[4][br.flavor_basis_pids.index(3)] = 1
+# s_+ = s + sbar
+ffns_rotate_to_LHA[5][br.flavor_basis_pids.index(-3)] = 1
+ffns_rotate_to_LHA[5][br.flavor_basis_pids.index(3)] = 1
+# c_+ = c + cbar
+ffns_rotate_to_LHA[6][br.flavor_basis_pids.index(-4)] = 1
+ffns_rotate_to_LHA[6][br.flavor_basis_pids.index(4)] = 1
+# g = g
+ffns_rotate_to_LHA[7][br.flavor_basis_pids.index(21)] = 1
+
+vfns_labels = ["u_v", "d_v", "L_m", "L_p", "s_p", "c_p", "b_p", "g"]
+vfns_rotate_to_LHA = np.zeros((len(vfns_labels), 14))
+# u_v = u - ubar
+vfns_rotate_to_LHA[0][br.flavor_basis_pids.index(-2)] = -1
+vfns_rotate_to_LHA[0][br.flavor_basis_pids.index(2)] = 1
+# d_v = d - dbar
+vfns_rotate_to_LHA[1][br.flavor_basis_pids.index(-1)] = -1
+vfns_rotate_to_LHA[1][br.flavor_basis_pids.index(1)] = 1
+# L_- = dbar - ubar
+vfns_rotate_to_LHA[2][br.flavor_basis_pids.index(-1)] = 1
+vfns_rotate_to_LHA[2][br.flavor_basis_pids.index(-2)] = -1
+# 2L_+ = 2dbar + 2ubar
+vfns_rotate_to_LHA[3][br.flavor_basis_pids.index(-1)] = 2
+vfns_rotate_to_LHA[3][br.flavor_basis_pids.index(-2)] = 2
+# s_+ = s + sbar
+vfns_rotate_to_LHA[4][br.flavor_basis_pids.index(-3)] = 1
+vfns_rotate_to_LHA[4][br.flavor_basis_pids.index(3)] = 1
+# c_+ = c + cbar
+vfns_rotate_to_LHA[5][br.flavor_basis_pids.index(-4)] = 1
+vfns_rotate_to_LHA[5][br.flavor_basis_pids.index(4)] = 1
+# b_+ = b + bbar
+vfns_rotate_to_LHA[6][br.flavor_basis_pids.index(-5)] = 1
+vfns_rotate_to_LHA[6][br.flavor_basis_pids.index(5)] = 1
+# g = g
+vfns_rotate_to_LHA[7][br.flavor_basis_pids.index(21)] = 1
diff --git a/extras/lh_bench_23/run.py b/extras/lh_bench_23/run-nnlo.py
similarity index 66%
rename from extras/lh_bench_23/run.py
rename to extras/lh_bench_23/run-nnlo.py
index b9146c1f1..1a006f7bc 100644
--- a/extras/lh_bench_23/run.py
+++ b/extras/lh_bench_23/run-nnlo.py
@@ -7,42 +7,24 @@
 import pandas as pd
 import yaml
 from banana import toy
-from cards import ffns_operator, ffns_theory, vfns_operator, vfns_theory
+from cfg import (
+    ffns_labels,
+    ffns_operator,
+    ffns_rotate_to_LHA,
+    ffns_theory,
+    vfns_labels,
+    vfns_operator,
+    vfns_rotate_to_LHA,
+    vfns_theory,
+)
 
 import eko
-from eko import basis_rotation as br
 from eko.runner.managed import solve
 from ekobox import apply
 from ekomark.benchmark.external.LHA_utils import here as there
 
 _sqrt2 = float(np.sqrt(2))
 
-# setup flavor rotations
-labels = ["u_v", "d_v", "L_m", "L_p", "s_v", "s_p", "c_p", "g"]
-rotate_to_LHA = np.zeros((len(labels), 14))
-# u_v = u - ubar
-rotate_to_LHA[0][br.flavor_basis_pids.index(-2)] = -1
-rotate_to_LHA[0][br.flavor_basis_pids.index(2)] = 1
-# d_v = d - dbar
-rotate_to_LHA[1][br.flavor_basis_pids.index(-1)] = -1
-rotate_to_LHA[1][br.flavor_basis_pids.index(1)] = 1
-# L_- = dbar - ubar
-rotate_to_LHA[2][br.flavor_basis_pids.index(-1)] = 1
-rotate_to_LHA[2][br.flavor_basis_pids.index(-2)] = -1
-# 2L_+ = 2dbar + 2ubar
-rotate_to_LHA[3][br.flavor_basis_pids.index(-1)] = 2
-rotate_to_LHA[3][br.flavor_basis_pids.index(-2)] = 2
-# s_v = s - sbar
-rotate_to_LHA[4][br.flavor_basis_pids.index(-3)] = -1
-rotate_to_LHA[4][br.flavor_basis_pids.index(3)] = 1
-# s_+ = s + sbar
-rotate_to_LHA[5][br.flavor_basis_pids.index(-3)] = 1
-rotate_to_LHA[5][br.flavor_basis_pids.index(3)] = 1
-# c_+ = c + cbar
-rotate_to_LHA[6][br.flavor_basis_pids.index(-4)] = 1
-rotate_to_LHA[6][br.flavor_basis_pids.index(4)] = 1
-# g = g
-rotate_to_LHA[7][br.flavor_basis_pids.index(21)] = 1
 
 # setup x rotation
 xgrid = np.array([1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.1, 0.3, 0.5, 0.7, 0.9])
@@ -86,11 +68,15 @@
         t = ffns_theory(xif)
         o = ffns_operator
         tab = 14
+        lab = ffns_labels
+        rot = ffns_rotate_to_LHA
     elif args.scheme == "VFNS":
         scheme = "VFNS"
         t = vfns_theory(xif)
         o = vfns_operator
         tab = 15
+        lab = vfns_labels
+        rot = vfns_rotate_to_LHA
     else:
         raise ValueError("scheme has to be FFNS or VFNS")
 
@@ -113,10 +99,8 @@
     # apply PDF
     out = {}
     with eko.EKO.read(p) as eko_:
-        pdf = apply.apply_pdf_flavor(
-            eko_, toy.mkPDF("ToyLH", 0), xgrid, rotate_to_LHA, labels
-        )
-        for lab, f in pdf[(10000.0, 4)]["pdfs"].items():
+        pdf = apply.apply_pdf_flavor(eko_, toy.mkPDF("ToyLH", 0), xgrid, rot, lab)
+        for lab, f in list(pdf.values())[0]["pdfs"].items():
             out[lab] = xgrid * f
 
     # display result
@@ -124,6 +108,8 @@
     me = pd.DataFrame(out)
     print("EKO")
     print(me)
+    # dump to file
+    me.to_csv(f"table{tab}-part{part}.csv")
 
     # load reference
     ref = pd.DataFrame(ref_data[f"table{tab}"][f"part{part}"])

From e8a007964616744acf09606261c0dfff56a95036 Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Wed, 20 Sep 2023 13:57:53 +0200
Subject: [PATCH 05/41] intorduce PTO_match

---
 extras/lh_bench_23/cfg.py                             | 1 +
 src/eko/evolution_operator/grid.py                    | 2 ++
 src/eko/evolution_operator/operator_matrix_element.py | 2 +-
 src/eko/io/runcards.py                                | 4 ++++
 src/eko/runner/legacy.py                              | 1 +
 src/eko/runner/parts.py                               | 1 +
 src/ekobox/cards.py                                   | 1 +
 7 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/extras/lh_bench_23/cfg.py b/extras/lh_bench_23/cfg.py
index ebafc1a13..51762de73 100644
--- a/extras/lh_bench_23/cfg.py
+++ b/extras/lh_bench_23/cfg.py
@@ -31,6 +31,7 @@
     ),
     xif=1.0,
     n3lo_ad_variation=(0, 0, 0, 0),
+    matching_order=[2, 0],
 )
 
 
diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py
index ca8a46d51..3e066a6de 100644
--- a/src/eko/evolution_operator/grid.py
+++ b/src/eko/evolution_operator/grid.py
@@ -54,6 +54,7 @@ def __init__(
         intrinsic_flavors: list,
         xif: float,
         n3lo_ad_variation: tuple,
+        matching_order: Order,
         configs: Configs,
         debug: Debug,
         atlas: Atlas,
@@ -81,6 +82,7 @@ def __init__(
         config["debug_skip_non_singlet"] = debug.skip_non_singlet
         config["polarized"] = configs.polarized
         config["time_like"] = configs.time_like
+        config["matching_order"] = matching_order
 
         if order == (1, 0) and method != "iterate-exact":
             logger.warning("Evolution: In LO we use the exact solution always!")
diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py
index 5696ead19..c25252675 100644
--- a/src/eko/evolution_operator/operator_matrix_element.py
+++ b/src/eko/evolution_operator/operator_matrix_element.py
@@ -220,7 +220,7 @@ def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar):
         self.L = L
         self.is_msbar = is_msbar
         # Note for the moment only QCD matching is implemented
-        self.order = (self.order[0] - 1, self.order[1])
+        self.order = tuple(config["matching_order"])
 
     @property
     def labels(self):
diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py
index 41aef9445..56d5cfd1a 100644
--- a/src/eko/io/runcards.py
+++ b/src/eko/io/runcards.py
@@ -49,6 +49,8 @@ class TheoryCard(DictLike):
     """Ratio between factorization scale and process scale."""
     n3lo_ad_variation: N3LOAdVariation
     """|N3LO| anomalous dimension variation: ``(gg_var, gq_var, qg_var, qq_var)``."""
+    matching_order: Order
+    """Matching conditions perturbative order tuple, ``(QCD, QED)``."""
 
 
 @dataclass
@@ -210,6 +212,8 @@ def new_theory(self):
 
         new["xif"] = old["XIF"]
         new["n3lo_ad_variation"] = old.get("n3lo_ad_variation", (0, 0, 0, 0))
+        # here PTO: 0 means truly LO
+        new["matching_order"] = old.get("PTO_matching", [old["PTO"], old["QED"]])
 
         return TheoryCard.from_dict(new)
 
diff --git a/src/eko/runner/legacy.py b/src/eko/runner/legacy.py
index c964c3415..991ec2bd0 100644
--- a/src/eko/runner/legacy.py
+++ b/src/eko/runner/legacy.py
@@ -77,6 +77,7 @@ def __init__(
             couplings=cs,
             interpol_dispatcher=bfd,
             n3lo_ad_variation=new_theory.n3lo_ad_variation,
+            matching_order=new_theory.matching_order,
         )
 
         with EKO.create(path) as builder:
diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py
index e385b8f49..9df80a4a8 100644
--- a/src/eko/runner/parts.py
+++ b/src/eko/runner/parts.py
@@ -88,6 +88,7 @@ def evolve_configs(eko: EKO) -> dict:
         n_integration_cores=ocard.configs.n_integration_cores,
         ModSV=ocard.configs.scvar_method,
         n3lo_ad_variation=tcard.n3lo_ad_variation,
+        matching_order=tcard.matching_order,
     )
 
 
diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py
index 474e7d931..559790d30 100644
--- a/src/ekobox/cards.py
+++ b/src/ekobox/cards.py
@@ -27,6 +27,7 @@
     ),
     xif=1.0,
     n3lo_ad_variation=(0, 0, 0, 0),
+    matching_order=[0, 0],
 )
 
 _operator = dict(

From 9eb1dc439049435deb3b1f4996b8ca0a5568f606 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 15:47:11 +0200
Subject: [PATCH 06/41] Interface QED and N3LO

---
 .../unpolarized/space_like/__init__.py        |  9 +++
 .../unpolarized/space_like/as4/__init__.py    | 81 +++++++++++++++++++
 2 files changed, 90 insertions(+)

diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
index 6a3c6cf06..0fbceb805 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
@@ -155,6 +155,11 @@ def gamma_ns_qed(order, mode, n, nf):
             gamma_ns[3, 0] = as3.gamma_nsp(n, nf, cache)
         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)
     return gamma_ns
 
 
@@ -280,6 +285,8 @@ def gamma_singlet_qed(order, n, nf):
         gamma_s[0, 2] = aem2.gamma_singlet(n, nf, cache)
     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)
     return gamma_s
 
 
@@ -314,4 +321,6 @@ def gamma_valence_qed(order, n, nf):
         gamma_v[0, 2] = aem2.gamma_valence(n, nf, cache)
     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)
     return gamma_v
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
index 10c4cb57c..56c2c16ff 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
@@ -55,3 +55,84 @@ 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
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(3,0)}(N)`
+
+    """
+    gamma_np_p = gamma_nsp(N, nf, cache)
+    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):
+    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
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{V}^{(3,0)}(N)`
+
+    """
+    gamma_V = np.array(
+        [
+            [gamma_nsv(N, nf, cache), 0.0],
+            [0.0, gamma_nsm(N, nf, cache)],
+        ],
+        np.complex_,
+    )
+    return gamma_V

From 089391fa2e9a0f6c4c24736ab923f858bb634e34 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 15:57:03 +0200
Subject: [PATCH 07/41] Add n3lo_ad_variation key to missing functions

---
 src/eko/evolution_operator/__init__.py                      | 5 ++++-
 .../anomalous_dimensions/unpolarized/space_like/__init__.py | 6 ++++--
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index eab3c64de..27b8933dc 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -450,6 +450,7 @@ def quad_ker_qed(
     ev_op_max_order,
     sv_mode,
     is_threshold,
+    n3lo_ad_variation,
 ):
     """Raw evolution kernel inside quad.
 
@@ -489,6 +490,8 @@ def quad_ker_qed(
         scale variation mode, see `eko.scale_variations.Modes`
     is_threshold : boolean
         is this an itermediate threshold operator?
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -497,7 +500,7 @@ 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)
+        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
index 0fbceb805..6358a5af8 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
@@ -255,7 +255,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache):
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(order, n, nf):
+def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
     r"""
     Compute the grid of the QED singlet anomalous dimensions matrices.
 
@@ -267,6 +267,8 @@ def gamma_singlet_qed(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -286,7 +288,7 @@ def gamma_singlet_qed(order, n, nf):
     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)
+        gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation)
     return gamma_s
 
 

From 9db7ac84a8345db686bdd5dbabeb9f9b0caa5c4b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:05:28 +0200
Subject: [PATCH 08/41] Add n3lo_ad_variation key to missing functions. again

---
 src/eko/evolution_operator/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 27b8933dc..122bfbc69 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -307,6 +307,7 @@ def quad_ker(
             ev_op_max_order,
             sv_mode,
             is_threshold,
+            n3lo_ad_variation,
         )
 
     # recombine everything

From 24bbfd30b5a5d0024abcda90b77bb7119591078a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:14:18 +0200
Subject: [PATCH 09/41] Fix test_kernels_QEDsinglet.py

---
 tests/eko/kernels/test_kernels_QEDsinglet.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/eko/kernels/test_kernels_QEDsinglet.py b/tests/eko/kernels/test_kernels_QEDsinglet.py
index 302bb941b..f5803750c 100644
--- a/tests/eko/kernels/test_kernels_QEDsinglet.py
+++ b/tests/eko/kernels/test_kernels_QEDsinglet.py
@@ -79,7 +79,7 @@ def test_zero_true_gamma(monkeypatch):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()
-            gamma_s = ad.gamma_singlet_qed(order, n, nf)
+            gamma_s = ad.gamma_singlet_qed(order, n, nf, (0, 0, 0, 0))
             # monkeypatch.setattr(
             #     ad,
             #     "exp_matrix_2D",

From 350871307f88cfebda1ac63446b6a00a23be918f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:30:38 +0200
Subject: [PATCH 10/41] Add order=4 in fixed_alphaem_exact

---
 src/eko/kernels/non_singlet_qed.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/eko/kernels/non_singlet_qed.py b/src/eko/kernels/non_singlet_qed.py
index 86ec1a9ff..36d4d3c35 100644
--- a/src/eko/kernels/non_singlet_qed.py
+++ b/src/eko/kernels/non_singlet_qed.py
@@ -120,6 +120,29 @@ def as3_exact(gamma_ns, a1, a0, beta):
     return ns.nnlo_exact(gamma_ns, a1, a0, beta)
 
 
+@nb.njit(cache=True)
+def as4_exact(gamma_ns, a1, a0, beta):
+    """O(as3aem1) non-singlet exact EKO.
+
+    Parameters
+    ----------
+    gamma_ns : numpy.ndarray
+        non-singlet anomalous dimensions
+    a1 : float
+        target coupling value
+    a0 : float
+        initial coupling value
+    beta : list
+        list of the values of the beta functions
+
+    Returns
+    -------
+    e_ns^3 : complex
+        O(as4aem1) non-singlet exact EKO
+    """
+    return ns.n3lo_exact(gamma_ns, a1, a0, beta)
+
+
 @nb.njit(cache=True)
 def dispatcher(
     order,
@@ -216,6 +239,8 @@ def fixed_alphaem_exact(order, gamma_ns, a1, a0, aem, nf, mu2_from, mu2_to):
         qcd_only = as2_exact(gamma_ns_list[1:], a1, a0, betalist)
     elif order[0] == 3:
         qcd_only = as3_exact(gamma_ns_list[1:], a1, a0, betalist)
+    elif order[0] == 4:
+        qcd_only = as4_exact(gamma_ns_list[1:], a1, a0, betalist)
     else:
         raise NotImplementedError("Selected order is not implemented")
     return qcd_only * apply_qed(gamma_ns_list[0], mu2_from, mu2_to)

From 0118d11807948113b59bc7c8149eb43990d52c04 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:45:39 +0200
Subject: [PATCH 11/41] Fix test_kernels_QEDns.py

---
 tests/eko/kernels/test_kernels_QEDns.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index e84956084..5574b1a12 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -304,7 +304,7 @@ def test_error():
     for running in [True, False]:
         with pytest.raises(NotImplementedError):
             ns.dispatcher(
-                (4, 2),
+                (5, 2),
                 "iterate-exact",
                 np.random.rand(4, 3),
                 [0.1, 0.2],

From 361afa8d2b4cef4843e459374cb6e45437991fba Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:56:16 +0200
Subject: [PATCH 12/41] Fix again test

---
 tests/eko/kernels/test_kernels_QEDns.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index 5574b1a12..f2bd45731 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -302,7 +302,7 @@ def test_ode_true_gamma():
 
 def test_error():
     for running in [True, False]:
-        with pytest.raises(NotImplementedError):
+        with pytest.raises(ValueError):
             ns.dispatcher(
                 (5, 2),
                 "iterate-exact",

From 17404d359beb475c67ea2a29a9f52be9ebf795c1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 17:09:22 +0200
Subject: [PATCH 13/41] Fix another test

---
 .../anomalous_dimensions/unpolarized/space_like/test_init.py    | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
index 2175eddcd..c7c955960 100644
--- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
@@ -208,7 +208,7 @@ def test_dim_singlet():
     nf = 3
     N = 2
     cache = h.cache.reset()
-    gamma_singlet = ad_us.gamma_singlet_qed((3, 2), N, nf)
+    gamma_singlet = ad_us.gamma_singlet_qed((3, 2), N, nf, (0, 0, 0, 0))
     assert gamma_singlet.shape == (4, 3, 4, 4)
     gamma_singlet_as1 = ad_us.as1.gamma_singlet_qed(N, cache, nf)
     assert gamma_singlet_as1.shape == (4, 4)

From 263f6063cca6e23535209feb2cc525f3bb14d40e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Sat, 7 Oct 2023 19:44:38 +0200
Subject: [PATCH 14/41] Change nf->beta into ns kernels at n3lo

---
 src/eko/kernels/non_singlet.py | 24 ++++++++----------------
 1 file changed, 8 insertions(+), 16 deletions(-)

diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py
index 24706a794..23b2ab47a 100644
--- a/src/eko/kernels/non_singlet.py
+++ b/src/eko/kernels/non_singlet.py
@@ -186,7 +186,7 @@ def nnlo_expanded(gamma_ns, a1, a0, beta):
 
 
 @nb.njit(cache=True)
-def n3lo_expanded(gamma_ns, a1, a0, nf):
+def n3lo_expanded(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet expanded EKO.
 
     Parameters
@@ -206,12 +206,8 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet expanded EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_expanded(a1, a0, beta0, b_list)
     j23 = as4_ei.j23_expanded(a1, a0, beta0, b_list)
@@ -225,7 +221,7 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
 
 
 @nb.njit(cache=True)
-def n3lo_exact(gamma_ns, a1, a0, nf):
+def n3lo_exact(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet exact EKO.
 
     Parameters
@@ -245,12 +241,8 @@ def n3lo_exact(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet exact EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     roots = as4_ei.roots(b_list)
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_exact(a1, a0, beta0, b_list, roots)
@@ -420,6 +412,6 @@ def dispatcher(
             "decompose-expanded",
             "perturbative-expanded",
         ]:
-            return n3lo_expanded(gamma_ns, a1, a0, nf)
-        return n3lo_exact(gamma_ns, a1, a0, nf)
+            return n3lo_expanded(gamma_ns, a1, a0, betalist)
+        return n3lo_exact(gamma_ns, a1, a0, betalist)
     raise NotImplementedError("Selected order is not implemented")

From 9393ea2c125cab80f970e62b9328f74059118458 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Sun, 8 Oct 2023 15:43:31 +0200
Subject: [PATCH 15/41] Add pto = 3 in some qed tests

---
 tests/eko/evolution_operator/test_init.py    |  2 +-
 tests/eko/kernels/test_kernels_QEDns.py      | 14 +++++++++-----
 tests/eko/kernels/test_kernels_QEDsinglet.py |  4 ++--
 tests/eko/kernels/test_kernels_QEDvalence.py |  4 ++--
 4 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index 642de5bb7..47d5e7001 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -382,7 +382,7 @@ def test_compute(self, monkeypatch, theory_ffns, operator_card, tmp_path):
                     err_msg=k,
                 )
 
-        for n in range(1, 3 + 1):
+        for n in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 g.config["order"] = (n, qed)
                 o1 = Operator(g.config, g.managers, Segment(2.0, 2.0, 3))
diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index f2bd45731..e6179ba99 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -25,7 +25,7 @@ def test_zero():
     nf = 3
     ev_op_iterations = 2
     running_alpha = [True, False]
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_ns = (
@@ -73,7 +73,7 @@ def test_zero_true_gamma():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        for qcd in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 n = np.random.rand()
@@ -136,7 +136,8 @@ def test_ode():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        max_qcd = 4
+        for qcd in range(1, max_qcd + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 sc = Couplings(
@@ -154,14 +155,17 @@ def test_ode():
                     )
                     a1 = sc.a_s(mu2_to)
                     gamma_ns = (
-                        np.random.rand(3 + 1, 2 + 1) + np.random.rand(3 + 1, 2 + 1) * 1j
+                        np.random.rand(max_qcd + 1, 2 + 1)
+                        + np.random.rand(max_qcd + 1, 2 + 1) * 1j
                     )
                     gamma_ns[0, 0] = 0.0
                     gamma_ns[2, 1] = 0.0
                     gamma_ns[3, 1] = 0.0
+                    gamma_ns[4, 1] = 0.0
                     gamma_ns[1, 2] = 0.0
                     gamma_ns[2, 2] = 0.0
                     gamma_ns[3, 2] = 0.0
+                    gamma_ns[4, 2] = 0.0
                     gammatot = 0.0
                     for i in range(0, order[0] + 1):
                         for j in range(0, order[1] + 1):
@@ -226,7 +230,7 @@ def test_ode_true_gamma():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        for qcd in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 sc = Couplings(
diff --git a/tests/eko/kernels/test_kernels_QEDsinglet.py b/tests/eko/kernels/test_kernels_QEDsinglet.py
index f5803750c..b953c053c 100644
--- a/tests/eko/kernels/test_kernels_QEDsinglet.py
+++ b/tests/eko/kernels/test_kernels_QEDsinglet.py
@@ -21,7 +21,7 @@ def test_zero(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_s = (
@@ -75,7 +75,7 @@ def test_zero_true_gamma(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()
diff --git a/tests/eko/kernels/test_kernels_QEDvalence.py b/tests/eko/kernels/test_kernels_QEDvalence.py
index 95d0cbc4c..a27bc6ef6 100644
--- a/tests/eko/kernels/test_kernels_QEDvalence.py
+++ b/tests/eko/kernels/test_kernels_QEDvalence.py
@@ -21,7 +21,7 @@ def test_zero(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_v = (
@@ -62,7 +62,7 @@ def test_zero_true_gamma(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()

From 66f02cc0f44482b2aab77e55d22ff275151d4719 Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Mon, 9 Oct 2023 10:19:23 +0200
Subject: [PATCH 16/41] add a failing bench

---
 benchmarks/eko/benchmark_iveverse_matching.py | 107 ++++++++++++++++++
 1 file changed, 107 insertions(+)
 create mode 100644 benchmarks/eko/benchmark_iveverse_matching.py

diff --git a/benchmarks/eko/benchmark_iveverse_matching.py b/benchmarks/eko/benchmark_iveverse_matching.py
new file mode 100644
index 000000000..18f715c10
--- /dev/null
+++ b/benchmarks/eko/benchmark_iveverse_matching.py
@@ -0,0 +1,107 @@
+import pathlib
+
+import numpy as np
+import pytest
+from banana import toy
+
+from eko import EKO
+from eko.io import runcards
+from eko.io.types import ReferenceRunning
+from eko.runner.managed import solve
+from ekobox import apply
+
+here = pathlib.Path(__file__).parent.absolute()
+MC = 1.51
+C_PID = 4
+
+# theory settings
+th_raw = dict(
+    order=[3, 0],
+    couplings=dict(
+        alphas=0.118,
+        alphaem=0.007496252,
+        scale=91.2,
+        num_flavs_ref=5,
+        max_num_flavs=6,
+    ),
+    heavy=dict(
+        num_flavs_init=4,
+        num_flavs_max_pdf=6,
+        intrinsic_flavors=[],
+        masses=[ReferenceRunning([mq, np.nan]) for mq in (MC, 4.92, 172.5)],
+        masses_scheme="POLE",
+        matching_ratios=[1.0, 1.0, np.inf],
+    ),
+    xif=1.0,
+    n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0),
+    matching_order=[2, 0],
+)
+
+# operator settings
+op_raw = dict(
+    mu0=1.65,
+    xgrid=[0.0001, 0.001, 0.01, 0.1, 1],
+    mugrid=[(MC, 3), (MC, 4)],
+    configs=dict(
+        evolution_method="truncated",
+        ev_op_max_order=[1, 0],
+        ev_op_iterations=1,
+        interpolation_polynomial_degree=4,
+        interpolation_is_log=True,
+        scvar_method="exponentiated",
+        inversion_method="exact",
+        n_integration_cores=0,
+        polarized=False,
+        time_like=False,
+    ),
+    debug=dict(
+        skip_singlet=False,
+        skip_non_singlet=False,
+    ),
+)
+
+
+@pytest.mark.isolated
+def BenchmarkInverseMatching():
+    th_card = runcards.TheoryCard.from_dict(th_raw)
+    op_card = runcards.OperatorCard.from_dict(op_raw)
+
+    eko_path2 = here / "test2.tar"
+    eko_path2.unlink(missing_ok=True)
+    solve(th_card, op_card, eko_path2)
+
+    th_card.matching_order = [1, 0]
+    eko_path1 = here / "test1.tar"
+    eko_path1.unlink(missing_ok=True)
+    solve(th_card, op_card, eko_path1)
+
+    eko_output1 = EKO.read(eko_path1)
+    eko_output2 = EKO.read(eko_path2)
+    op1_nf3 = eko_output1[(MC**2, 3)]
+    op2_nf3 = eko_output2[(MC**2, 3)]
+    op1_nf4 = eko_output1[(MC**2, 4)]
+    op2_nf4 = eko_output2[(MC**2, 4)]
+
+    # test that nf=4 operators are the same
+    np.testing.assert_allclose(op1_nf4.operator, op2_nf4.operator)
+
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(op2_nf3.operator, op2_nf4.operator)
+
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(op1_nf3.operator, op1_nf4.operator)
+
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(op1_nf3.operator, op2_nf3.operator)
+
+    pdf1 = apply.apply_pdf(eko_output1, toy.mkPDF("ToyLH", 0))
+    pdf2 = apply.apply_pdf(eko_output2, toy.mkPDF("ToyLH", 0))
+
+    # test that different PTO matching is applied correctly
+    np.testing.assert_allclose(
+        pdf1[(MC**2, 4)]["pdfs"][C_PID], pdf2[(MC**2, 4)]["pdfs"][C_PID]
+    )
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(
+            pdf1[(MC**2, 3)]["pdfs"][C_PID], pdf2[(MC**2, 3)]["pdfs"][C_PID]
+        )

From c712cad1c1c30244411db5f29bee0ceef1bb617a Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Mon, 9 Oct 2023 10:38:08 +0200
Subject: [PATCH 17/41] fix bench name

---
 benchmarks/eko/benchmark_iveverse_matching.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmarks/eko/benchmark_iveverse_matching.py b/benchmarks/eko/benchmark_iveverse_matching.py
index 18f715c10..7ca3b1a6a 100644
--- a/benchmarks/eko/benchmark_iveverse_matching.py
+++ b/benchmarks/eko/benchmark_iveverse_matching.py
@@ -62,7 +62,7 @@
 
 
 @pytest.mark.isolated
-def BenchmarkInverseMatching():
+def benchmark_inverse_matching():
     th_card = runcards.TheoryCard.from_dict(th_raw)
     op_card = runcards.OperatorCard.from_dict(op_raw)
 

From 9d74a9831f8729bd69eb2b8d415c5bd140f7d956 Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Mon, 9 Oct 2023 14:09:44 +0200
Subject: [PATCH 18/41] separate nf black in evolve_pdfs

---
 src/ekobox/evol_pdf.py | 40 ++++++++++++++++++++++++----------------
 1 file changed, 24 insertions(+), 16 deletions(-)

diff --git a/src/ekobox/evol_pdf.py b/src/ekobox/evol_pdf.py
index e6c765187..1f6b4b5df 100644
--- a/src/ekobox/evol_pdf.py
+++ b/src/ekobox/evol_pdf.py
@@ -1,5 +1,6 @@
 """Tools to evolve actual PDFs."""
 import pathlib
+from collections import defaultdict
 
 import numpy as np
 
@@ -78,24 +79,31 @@ def evolve_pdfs(
     )
     all_member_blocks = []
     targetlist = targetgrid.raw.tolist()
+
+    # separate by nf the evolgrid (and order per nf/q)
+    by_nf = defaultdict(list)
+    for q, nf in sorted(eko_output.evolgrid, key=lambda ep: ep[1]):
+        by_nf[nf].append(q)
+    q2block_per_nf = {nf: sorted(qs) for nf, qs in by_nf.items()}
+
+    # loop on replicas
     for evolved_PDF in evolved_PDF_list:
         all_blocks = []
-        evolved_PDF_q2 = {q2: val for (q2, _), val in evolved_PDF.items()}
-        sorted_q2grid = [
-            float(q2) for q2, _ in np.sort(operators_card.evolgrid, axis=0)
-        ]
-        block = genpdf.generate_block(
-            lambda pid, x, Q2, evolved_PDF=evolved_PDF_q2: targetlist[
-                targetlist.index(x)
-            ]
-            * evolved_PDF[Q2]["pdfs"][pid][targetlist.index(x)],
-            xgrid=targetlist,
-            sorted_q2grid=sorted_q2grid,
-            pids=np.array(br.flavor_basis_pids),
-        )
-        # all_blocks will be useful in case there will be necessity to dump many blocks
-        # for a single member
-        all_blocks.append(block)
+
+        # loop on nf patches
+        for nf, q2grid in q2block_per_nf.items():
+
+            def pdf_xq2(pid, x, Q2):
+                x_idx = targetlist.index(x)
+                return x * evolved_PDF[(Q2, nf)]["pdfs"][pid][x_idx]
+
+            block = genpdf.generate_block(
+                pdf_xq2,
+                xgrid=targetlist,
+                sorted_q2grid=q2grid,
+                pids=br.flavor_basis_pids,
+            )
+            all_blocks.append(block)
         all_member_blocks.append(all_blocks)
 
     genpdf.export.dump_set(name, info, all_member_blocks)

From df26b330a0690f92b2b6afe99905d1cddb8d6129 Mon Sep 17 00:00:00 2001
From: Giacomo Magni <39065935+giacomomagni@users.noreply.github.com>
Date: Mon, 9 Oct 2023 16:04:06 +0200
Subject: [PATCH 19/41] Update src/eko/io/runcards.py

Co-authored-by: Alessandro Candido <candido.ale@gmail.com>
---
 src/eko/io/runcards.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py
index 56d5cfd1a..e9260a6eb 100644
--- a/src/eko/io/runcards.py
+++ b/src/eko/io/runcards.py
@@ -49,7 +49,7 @@ class TheoryCard(DictLike):
     """Ratio between factorization scale and process scale."""
     n3lo_ad_variation: N3LOAdVariation
     """|N3LO| anomalous dimension variation: ``(gg_var, gq_var, qg_var, qq_var)``."""
-    matching_order: Order
+    matching_order: Optional[Order] = None
     """Matching conditions perturbative order tuple, ``(QCD, QED)``."""
 
 

From 6d1452978e20bb98f27edfa830e698069aa1d617 Mon Sep 17 00:00:00 2001
From: Giacomo Magni <39065935+giacomomagni@users.noreply.github.com>
Date: Mon, 9 Oct 2023 16:04:13 +0200
Subject: [PATCH 20/41] Update src/eko/runner/parts.py

Co-authored-by: Alessandro Candido <candido.ale@gmail.com>
---
 src/eko/runner/parts.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py
index 9df80a4a8..1bb69aab9 100644
--- a/src/eko/runner/parts.py
+++ b/src/eko/runner/parts.py
@@ -88,7 +88,7 @@ def evolve_configs(eko: EKO) -> dict:
         n_integration_cores=ocard.configs.n_integration_cores,
         ModSV=ocard.configs.scvar_method,
         n3lo_ad_variation=tcard.n3lo_ad_variation,
-        matching_order=tcard.matching_order,
+        matching_order=tcard.matching_order if tcard.matching_order is not None else tcard.order,
     )
 
 

From 8e9cf99ac03d176d46c1d3f1ab03f180b580ca4c Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Mon, 9 Oct 2023 16:07:24 +0200
Subject: [PATCH 21/41] add default matching_order to legacy runner

---
 src/eko/runner/legacy.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/eko/runner/legacy.py b/src/eko/runner/legacy.py
index 991ec2bd0..ff1bb1e7b 100644
--- a/src/eko/runner/legacy.py
+++ b/src/eko/runner/legacy.py
@@ -77,7 +77,9 @@ def __init__(
             couplings=cs,
             interpol_dispatcher=bfd,
             n3lo_ad_variation=new_theory.n3lo_ad_variation,
-            matching_order=new_theory.matching_order,
+            matching_order=new_theory.matching_order
+            if new_theory.matching_order is not None
+            else new_theory.order,
         )
 
         with EKO.create(path) as builder:

From e0d74a75b6ab1f7606e48e9d2568cd8f32ad34ca Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Mon, 9 Oct 2023 16:13:24 +0200
Subject: [PATCH 22/41] add fix and correct file name

---
 ...mark_iveverse_matching.py => benchmark_inverse_matching.py} | 0
 src/eko/runner/parts.py                                        | 3 +--
 2 files changed, 1 insertion(+), 2 deletions(-)
 rename benchmarks/eko/{benchmark_iveverse_matching.py => benchmark_inverse_matching.py} (100%)

diff --git a/benchmarks/eko/benchmark_iveverse_matching.py b/benchmarks/eko/benchmark_inverse_matching.py
similarity index 100%
rename from benchmarks/eko/benchmark_iveverse_matching.py
rename to benchmarks/eko/benchmark_inverse_matching.py
diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py
index 9df80a4a8..ce540d861 100644
--- a/src/eko/runner/parts.py
+++ b/src/eko/runner/parts.py
@@ -150,9 +150,8 @@ def match(eko: EKO, recipe: Matching) -> Operator:
     op.compute()
 
     binfo = blowup_info(eko)
-    nf_match = op.nf - 1 if recipe.inverse else op.nf
     res, err = matching_condition.MatchingCondition.split_ad_to_evol_map(
-        op.op_members, nf_match, recipe.scale, **binfo
+        op.op_members, op.nf, recipe.scale, **binfo
     ).to_flavor_basis_tensor(qed=binfo["qed"])
 
     return Operator(res, err)

From ce71c5d5fbd12f8c911502a61cf3f8b053d7b21c Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Mon, 9 Oct 2023 16:44:07 +0200
Subject: [PATCH 23/41] run pre-commit

---
 src/eko/runner/parts.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py
index 1bb69aab9..e2c09bf09 100644
--- a/src/eko/runner/parts.py
+++ b/src/eko/runner/parts.py
@@ -88,7 +88,9 @@ def evolve_configs(eko: EKO) -> dict:
         n_integration_cores=ocard.configs.n_integration_cores,
         ModSV=ocard.configs.scvar_method,
         n3lo_ad_variation=tcard.n3lo_ad_variation,
-        matching_order=tcard.matching_order if tcard.matching_order is not None else tcard.order,
+        matching_order=tcard.matching_order
+        if tcard.matching_order is not None
+        else tcard.order,
     )
 
 

From a05c58ca3300469f3cb2d6d7af03763d8ba27d52 Mon Sep 17 00:00:00 2001
From: Felix Hekhorn <felix.hekhorn@uni-tuebingen.de>
Date: Mon, 9 Oct 2023 17:57:14 +0300
Subject: [PATCH 24/41] Tear evolve_pdfs apart

---
 src/ekobox/evol_pdf.py | 72 ++++++++++++++++++++++++++++--------------
 1 file changed, 48 insertions(+), 24 deletions(-)

diff --git a/src/ekobox/evol_pdf.py b/src/ekobox/evol_pdf.py
index 1f6b4b5df..49cda53b6 100644
--- a/src/ekobox/evol_pdf.py
+++ b/src/ekobox/evol_pdf.py
@@ -58,6 +58,7 @@ def evolve_pdfs(
         eko.solve(theory_card, operators_card, path=store_path)
         eko_path = store_path
 
+    # apply PDF to eko
     evolved_PDF_list = []
     with EKO.read(eko_path) as eko_output:
         for initial_PDF in initial_PDF_list:
@@ -65,6 +66,7 @@ def evolve_pdfs(
                 apply.apply_pdf(eko_output, initial_PDF, targetgrid)
             )
 
+    # update info file
     if targetgrid is None:
         targetgrid = operators_card.xgrid
     if info_update is None:
@@ -77,36 +79,58 @@ def evolve_pdfs(
         len(evolved_PDF_list),
         info_update=info_update,
     )
-    all_member_blocks = []
-    targetlist = targetgrid.raw.tolist()
 
     # separate by nf the evolgrid (and order per nf/q)
-    by_nf = defaultdict(list)
-    for q, nf in sorted(eko_output.evolgrid, key=lambda ep: ep[1]):
-        by_nf[nf].append(q)
-    q2block_per_nf = {nf: sorted(qs) for nf, qs in by_nf.items()}
+    q2block_per_nf = regroup_evolgrid(eko_output.evolgrid)
 
-    # loop on replicas
+    # write all replicas
+    all_member_blocks = []
     for evolved_PDF in evolved_PDF_list:
-        all_blocks = []
-
-        # loop on nf patches
-        for nf, q2grid in q2block_per_nf.items():
-
-            def pdf_xq2(pid, x, Q2):
-                x_idx = targetlist.index(x)
-                return x * evolved_PDF[(Q2, nf)]["pdfs"][pid][x_idx]
-
-            block = genpdf.generate_block(
-                pdf_xq2,
-                xgrid=targetlist,
-                sorted_q2grid=q2grid,
-                pids=br.flavor_basis_pids,
-            )
-            all_blocks.append(block)
+        all_blocks = collect_blocks(
+            evolved_PDF, q2block_per_nf, targetgrid.raw.tolist()
+        )
         all_member_blocks.append(all_blocks)
-
     genpdf.export.dump_set(name, info, all_member_blocks)
 
     if install:
         genpdf.install_pdf(name)
+
+
+def regroup_evolgrid(evolgrid: list):
+    """Split evolution points by nf and sort by scale."""
+    by_nf = defaultdict(list)
+    for q, nf in sorted(evolgrid, key=lambda ep: ep[1]):
+        by_nf[nf].append(q)
+    return {nf: sorted(qs) for nf, qs in by_nf.items()}
+
+
+def collect_blocks(evolved_PDF: dict, q2block_per_nf: dict, xgrid: list):
+    """Collect all LHAPDF blocks for a given replica.
+
+    Parameters
+    ----------
+    evolved_PDF :
+        PDF evaluated at grid
+    q2block_per_nf :
+        block coordinates
+    xgrid :
+        x grid
+
+    """
+    all_blocks = []
+
+    # fake xfxQ2
+    def pdf_xq2(pid, x, Q2):
+        x_idx = xgrid.index(x)
+        return x * evolved_PDF[(Q2, nf)]["pdfs"][pid][x_idx]
+
+    # loop on nf patches
+    for nf, q2grid in q2block_per_nf.items():
+        block = genpdf.generate_block(
+            pdf_xq2,
+            xgrid=xgrid,
+            sorted_q2grid=q2grid,
+            pids=br.flavor_basis_pids,
+        )
+        all_blocks.append(block)
+    return all_blocks

From a7f5d7a91ee7d0afc08f88cb3e3e517c67f5ba2b Mon Sep 17 00:00:00 2001
From: Felix Hekhorn <felix.hekhorn@uni-tuebingen.de>
Date: Mon, 9 Oct 2023 18:16:50 +0300
Subject: [PATCH 25/41] Add more unit tests for evol_gri

---
 tests/ekobox/test_evol_pdf.py | 57 +++++++++++++++++++++++++++++++++++
 1 file changed, 57 insertions(+)

diff --git a/tests/ekobox/test_evol_pdf.py b/tests/ekobox/test_evol_pdf.py
index 1c51dfdad..54bce6623 100644
--- a/tests/ekobox/test_evol_pdf.py
+++ b/tests/ekobox/test_evol_pdf.py
@@ -1,7 +1,9 @@
+import numpy as np
 from banana import toy
 
 import eko
 from eko import EKO
+from eko import basis_rotation as br
 from eko.interpolation import XGrid
 from ekobox import cards
 from ekobox import evol_pdf as ev_p
@@ -58,3 +60,58 @@ def test_evolve_pdfs_dump_file(fake_lhapdf, cd):
         ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=peko)
     p = fake_lhapdf / n
     assert p.exists()
+
+
+def test_regroup_evolgrid():
+    # basic
+    i = [(3.0, 3), (4.0, 3)]
+    o = ev_p.regroup_evolgrid(i)
+    assert len(o.keys()) == 1
+    assert 3 in o
+    assert len(o[3]) == 2
+    # more advanced
+    i = [(4.0, 3), (3.0, 3), (4.0, 4), (3.0, 4)]
+    o = ev_p.regroup_evolgrid(i)
+    assert len(o.keys()) == 2
+    assert 3 in o
+    assert 4 in o
+    assert len(o[3]) == 2
+    assert len(o[4]) == 2
+    np.testing.assert_allclose(o[3], o[4])
+    # messed up
+    i = [(4.0, 3), (4.0, 4), (3.0, 4), (3.0, 3), (5.0, 5)]
+    o = ev_p.regroup_evolgrid(i)
+    assert len(o.keys()) == 3
+    assert 3 in o
+    assert 4 in o
+    assert 5 in o
+    assert len(o[3]) == 2
+    assert len(o[4]) == 2
+    assert len(o[5]) == 1
+    np.testing.assert_allclose(o[3], o[4])
+
+
+def test_collect_blocks():
+    xgrid = [0.1, 0.5, 0.1]
+
+    def mk(eps):
+        f = {}
+        for ep in eps:
+            f[ep] = {
+                "pdfs": {
+                    pid: np.random.rand(len(xgrid)) for pid in br.flavor_basis_pids
+                }
+            }
+        return f
+
+    # basic
+    eps = [(3.0, 3), (4.0, 3)]
+    bs = ev_p.collect_blocks(mk(eps), ev_p.regroup_evolgrid(eps), xgrid)
+    assert len(bs) == 1
+    np.testing.assert_allclose(bs[0]["mu2grid"], (3.0, 4.0))
+    # more advanced
+    eps = [(4.0, 3), (3.0, 3), (5.0, 4), (3.0, 4)]
+    bs = ev_p.collect_blocks(mk(eps), ev_p.regroup_evolgrid(eps), xgrid)
+    assert len(bs) == 2
+    np.testing.assert_allclose(bs[0]["mu2grid"], (3.0, 4.0))
+    np.testing.assert_allclose(bs[1]["mu2grid"], (3.0, 5.0))

From e0e4ea6d6ae4e429918625b69cc0ebd57585af5e Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 9 Oct 2023 18:04:49 +0000
Subject: [PATCH 26/41] [pre-commit.ci] pre-commit autoupdate
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

updates:
- [github.com/pre-commit/pre-commit-hooks: v4.4.0 → v4.5.0](https://github.com/pre-commit/pre-commit-hooks/compare/v4.4.0...v4.5.0)
- [github.com/asottile/pyupgrade: v3.13.0 → v3.15.0](https://github.com/asottile/pyupgrade/compare/v3.13.0...v3.15.0)
---
 .pre-commit-config.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index afebc737d..985db8a0d 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -4,7 +4,7 @@ ci:
   autofix_prs: false
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0
+    rev: v4.5.0
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer
@@ -31,7 +31,7 @@ repos:
       - id: isort
         args: ["--profile", "black"]
   - repo: https://github.com/asottile/pyupgrade
-    rev: v3.13.0
+    rev: v3.15.0
     hooks:
       - id: pyupgrade
   - repo: https://github.com/pycqa/pydocstyle

From 91e59c728d8a563bebb9fefcf5667173afcd5db1 Mon Sep 17 00:00:00 2001
From: giacomomagni <giac.magni@gmail.com>
Date: Tue, 10 Oct 2023 15:52:50 +0200
Subject: [PATCH 27/41] minor fix on docstring

---
 src/eko/evolution_operator/operator_matrix_element.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py
index c25252675..1a287b9aa 100644
--- a/src/eko/evolution_operator/operator_matrix_element.py
+++ b/src/eko/evolution_operator/operator_matrix_element.py
@@ -185,8 +185,12 @@ class OperatorMatrixElement(Operator):
         configuration
     managers : dict
         managers
-    segment: Segment
-        path segment
+    nf: int
+        number of active flavor below threshold
+    q2: float
+        squared matching scale
+    is_backward: bool
+        True for backward matching
     L: float
         :math:`\ln(\mu_F^2 / m_h^2)`
     is_msbar: bool

From 569a0fe8f972eaf0b71549c0f8969542b9e38336 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 10 Oct 2023 16:42:26 +0200
Subject: [PATCH 28/41] restructure self.compute_aem_list

---
 src/eko/evolution_operator/__init__.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 122bfbc69..29e67c19c 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -626,7 +626,7 @@ def __init__(
         self.alphaem_running = self.managers["couplings"].alphaem_running
         if self.log_label == "Evolution":
             self.a = self.compute_a()
-            self.compute_aem_list()
+            self.as_list, self.a_half_list = self.compute_aem_list()
 
     @property
     def n_pools(self):
@@ -701,8 +701,8 @@ def compute_aem_list(self):
         """
         ev_op_iterations = self.config["ev_op_iterations"]
         if self.order[1] == 0:
-            self.as_list = np.array([self.a_s[0], self.a_s[1]])
-            self.a_half_list = np.zeros((ev_op_iterations, 2))
+            as_list = np.array([self.a_s[0], self.a_s[1]])
+            a_half = np.zeros((ev_op_iterations, 2))
         else:
             as0 = self.a_s[0]
             as1 = self.a_s[1]
@@ -722,7 +722,7 @@ def compute_aem_list(self):
             couplings = self.managers["couplings"]
             mu2_steps = utils.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations)
             mu2_l = mu2_steps[0]
-            self.as_list = np.array(
+            as_list = np.array(
                 [
                     couplings.compute(
                         a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2
@@ -738,7 +738,7 @@ def compute_aem_list(self):
                 )
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
-            self.a_half_list = a_half
+        return as_list, a_half
 
     @property
     def labels(self):

From 2f65ee975f869bad50e79c17fea3a3a9a13427e8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 10 Oct 2023 17:41:23 +0200
Subject: [PATCH 29/41] Remove unused default parameters and remove qed
 argument from apply_pdf

---
 src/eko/basis_rotation.py                        |  4 ++--
 src/eko/evolution_operator/flavors.py            |  8 ++++----
 src/eko/evolution_operator/matching_condition.py |  2 +-
 src/eko/evolution_operator/physical.py           |  2 +-
 src/ekobox/apply.py                              | 11 ++++-------
 src/ekomark/benchmark/runner.py                  |  1 -
 tests/eko/evolution_operator/test_flavors.py     | 16 +++++++++-------
 .../test_matching_condition.py                   | 10 +++++-----
 tests/eko/test_basis_rotation.py                 |  8 ++++----
 9 files changed, 30 insertions(+), 32 deletions(-)

diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py
index 4b972592c..e31a7beb6 100644
--- a/src/eko/basis_rotation.py
+++ b/src/eko/basis_rotation.py
@@ -275,7 +275,7 @@
 }
 
 
-def ad_projector(ad_lab, nf, qed=False):
+def ad_projector(ad_lab, nf, qed):
     """
     Build a projector (as a numpy array) for the given anomalous dimension sector.
 
@@ -355,7 +355,7 @@ def select_light_flavors_uni_ev(ad_lab, nf):
         return map_ad_to_evolution[ad_lab]
 
 
-def ad_projectors(nf, qed=False):
+def ad_projectors(nf, qed):
     """
     Build projectors tensor (as a numpy array), collecting all the individual sector projectors.
 
diff --git a/src/eko/evolution_operator/flavors.py b/src/eko/evolution_operator/flavors.py
index 7ad7b5801..fac94a608 100644
--- a/src/eko/evolution_operator/flavors.py
+++ b/src/eko/evolution_operator/flavors.py
@@ -52,7 +52,7 @@ def pids_from_intrinsic_evol(label, nf, normalize):
     return weights
 
 
-def get_range(evol_labels, qed=False):
+def get_range(evol_labels, qed):
     """Determine the number of light and heavy flavors participating in the input and output.
 
     Here, we assume that the T distributions (e.g. T15) appears *always*
@@ -73,7 +73,7 @@ def get_range(evol_labels, qed=False):
     nf_in = 3
     nf_out = 3
 
-    def update(label, qed=False):
+    def update(label, qed):
         nf = 3
         if label[0] == "T":
             if not qed:
@@ -129,7 +129,7 @@ def rotate_pm_to_flavor(label):
     return l
 
 
-def rotate_matching(nf, qed=False, inverse=False):
+def rotate_matching(nf, qed, inverse=False):
     """Rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
@@ -206,7 +206,7 @@ def rotate_matching(nf, qed=False, inverse=False):
     return l
 
 
-def rotate_matching_inverse(nf, qed=False):
+def rotate_matching_inverse(nf, qed):
     """Inverse rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
diff --git a/src/eko/evolution_operator/matching_condition.py b/src/eko/evolution_operator/matching_condition.py
index 059617fbf..8df4775ff 100644
--- a/src/eko/evolution_operator/matching_condition.py
+++ b/src/eko/evolution_operator/matching_condition.py
@@ -21,7 +21,7 @@ def split_ad_to_evol_map(
         nf,
         q2_thr,
         intrinsic_range,
-        qed=False,
+        qed,
     ):
         """
         Create the instance from the |OME|.
diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py
index c963e2a26..1ee501a2a 100644
--- a/src/eko/evolution_operator/physical.py
+++ b/src/eko/evolution_operator/physical.py
@@ -22,7 +22,7 @@ class PhysicalOperator(member.OperatorBase):
     """
 
     @classmethod
-    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False):
+    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed):
         """
         Obtain map between the 3-dimensional anomalous dimension basis and the 4-dimensional evolution basis.
 
diff --git a/src/ekobox/apply.py b/src/ekobox/apply.py
index c392ffe78..f59928186 100644
--- a/src/ekobox/apply.py
+++ b/src/ekobox/apply.py
@@ -12,7 +12,6 @@ def apply_pdf(
     lhapdf_like,
     targetgrid=None,
     rotate_to_evolution_basis=False,
-    qed=False,
 ):
     """
     Apply all available operators to the input PDFs.
@@ -34,13 +33,14 @@ def apply_pdf(
         out_grid : dict
             output PDFs and their associated errors for the computed mu2grid
     """
+    qed = eko.theory_card.order[1] > 0
     if rotate_to_evolution_basis:
         if not qed:
             rotate_flavor_to_evolution = br.rotate_flavor_to_evolution
         else:
             rotate_flavor_to_evolution = br.rotate_flavor_to_unified_evolution
         return apply_pdf_flavor(
-            eko, lhapdf_like, targetgrid, rotate_flavor_to_evolution, qed
+            eko, lhapdf_like, targetgrid, rotate_flavor_to_evolution
         )
     return apply_pdf_flavor(eko, lhapdf_like, targetgrid)
 
@@ -48,9 +48,7 @@ def apply_pdf(
 CONTRACTION = "ajbk,bk"
 
 
-def apply_pdf_flavor(
-    eko: EKO, lhapdf_like, targetgrid=None, flavor_rotation=None, qed=False
-):
+def apply_pdf_flavor(eko: EKO, lhapdf_like, targetgrid=None, flavor_rotation=None):
     """
     Apply all available operators to the input PDFs.
 
@@ -65,8 +63,6 @@ def apply_pdf_flavor(
             if given, interpolates to the pdfs given at targetgrid (instead of xgrid)
         flavor_rotation : np.ndarray
             Rotation matrix in flavor space
-        qed : bool
-            activate qed
 
     Returns
     -------
@@ -97,6 +93,7 @@ def apply_pdf_flavor(
         if error_final is not None:
             out_grid[ep]["errors"] = dict(zip(eko.bases.targetpids, error_final))
 
+    qed = eko.theory_card.order[1] > 0
     # rotate to evolution basis
     if flavor_rotation is not None:
         for q2, op in out_grid.items():
diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py
index 8c07e7153..e9ae298ee 100644
--- a/src/ekomark/benchmark/runner.py
+++ b/src/ekomark/benchmark/runner.py
@@ -243,7 +243,6 @@ def log(self, theory, _, pdf, me, ext):
                 pdf,
                 xgrid,
                 flavor_rotation=rotate_to_evolution,
-                qed=qed,
             )
         for q2, ref_pdfs in ext["values"].items():
             log_tab = dfdict.DFdict()
diff --git a/tests/eko/evolution_operator/test_flavors.py b/tests/eko/evolution_operator/test_flavors.py
index 81bb35fc4..2f3761404 100644
--- a/tests/eko/evolution_operator/test_flavors.py
+++ b/tests/eko/evolution_operator/test_flavors.py
@@ -61,13 +61,15 @@ def get(d):
 
 
 def test_get_range():
-    assert (3, 3) == flavors.get_range([])
-    assert (3, 3) == flavors.get_range([member.MemberName(n) for n in ["S.S", "V3.V3"]])
+    assert (3, 3) == flavors.get_range([], False)
+    assert (3, 3) == flavors.get_range(
+        [member.MemberName(n) for n in ["S.S", "V3.V3"]], False
+    )
     assert (3, 4) == flavors.get_range(
-        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.S"]]
+        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.S"]], False
     )
     assert (4, 4) == flavors.get_range(
-        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.T15"]]
+        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.T15"]], False
     )
     assert (3, 3) == flavors.get_range(
         [member.MemberName(n) for n in ["S.S", "Td3.Td3"]], True
@@ -113,7 +115,7 @@ def test_rotate_pm_to_flavor():
 
 
 def test_rotate_matching():
-    m = flavors.rotate_matching(4)
+    m = flavors.rotate_matching(4, False)
     assert len(list(filter(lambda e: "c+" in e, m.keys()))) == 2
     assert len(list(filter(lambda e: "b-" in e, m.keys()))) == 1
 
@@ -156,8 +158,8 @@ def load(m):
         return mm
 
     for nf in range(4, 6 + 1):
-        m = load(flavors.rotate_matching(nf))
-        minv = load(flavors.rotate_matching_inverse(nf))
+        m = load(flavors.rotate_matching(nf, False))
+        minv = load(flavors.rotate_matching_inverse(nf, False))
         np.testing.assert_allclose(m @ minv, np.eye(len(br.evol_basis)), atol=1e-10)
 
 
diff --git a/tests/eko/evolution_operator/test_matching_condition.py b/tests/eko/evolution_operator/test_matching_condition.py
index 8cbd6b537..c854fb90c 100644
--- a/tests/eko/evolution_operator/test_matching_condition.py
+++ b/tests/eko/evolution_operator/test_matching_condition.py
@@ -38,7 +38,7 @@ def update_intrinsic_OME(self, ome):
 
     def test_split_ad_to_evol_map(self):
         ome = self.mkOME()
-        a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [])
+        a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False)
         triv_keys = [
             "V.V",
             "T3.T3",
@@ -64,7 +64,7 @@ def test_split_ad_to_evol_map(self):
             ome[(200, 200)].value,
         )
         # # if alpha is zero, nothing non-trivial should happen
-        b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [])
+        b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False)
         assert sorted(str(k) for k in b.op_members.keys()) == sorted(
             [*triv_keys, *keys3]
         )
@@ -74,7 +74,7 @@ def test_split_ad_to_evol_map(self):
         # )
         # nf=3 + IC
         self.update_intrinsic_OME(ome)
-        c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4])
+        c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4], False)
         assert sorted(str(k) for k in c.op_members.keys()) == sorted(
             [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"]
         )
@@ -83,7 +83,7 @@ def test_split_ad_to_evol_map(self):
             b.op_members[member.MemberName("V.V")].value,
         )
         # nf=3 + IB
-        d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5])
+        d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5], False)
         assert sorted(str(k) for k in d.op_members.keys()) == sorted(
             [*triv_keys, *keys3, "b+.b+", "b-.b-"]
         )
@@ -92,7 +92,7 @@ def test_split_ad_to_evol_map(self):
             np.eye(self.shape[0]),
         )
         # nf=4 + IB
-        d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5])
+        d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5], False)
         assert sorted(str(k) for k in d.op_members.keys()) == sorted(
             [
                 *triv_keys,
diff --git a/tests/eko/test_basis_rotation.py b/tests/eko/test_basis_rotation.py
index 680ef065d..194669d46 100644
--- a/tests/eko/test_basis_rotation.py
+++ b/tests/eko/test_basis_rotation.py
@@ -9,19 +9,19 @@ def test_ad_projector():
     g = br.rotate_flavor_to_evolution[2]
     v3 = br.rotate_flavor_to_evolution[br.evol_basis.index("V3")]
 
-    s_to_s = br.ad_projector((100, 100), nf=6)
+    s_to_s = br.ad_projector((100, 100), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ s_to_s, s)
     np.testing.assert_allclose(g @ s_to_s, 0.0)
     np.testing.assert_allclose(v3 @ s_to_s, 0.0)
 
-    g_to_s = br.ad_projector((21, 100), nf=6)
+    g_to_s = br.ad_projector((21, 100), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ g_to_s, 0.0)
     np.testing.assert_allclose(g @ g_to_s, s)
     np.testing.assert_allclose(v3 @ g_to_s, 0.0)
 
-    ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6)
+    ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15)
     np.testing.assert_allclose(g @ ns_m, 0.0)
@@ -74,7 +74,7 @@ def test_ad_projectors():
     for nf in range(3, 6 + 1):
         diag = np.array([0] * (1 + 6 - nf) + [1] * (1 + 2 * nf) + [0] * (6 - nf))
         identity = np.diag(diag)
-        projs = br.ad_projectors(nf)
+        projs = br.ad_projectors(nf, qed=False)
 
         # sum over diagonal projectors form an identity
         np.testing.assert_allclose(

From 0d52dac77af2385e43edeae223132c5afca8b906 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 15:47:11 +0200
Subject: [PATCH 30/41] Interface QED and N3LO

---
 .../unpolarized/space_like/__init__.py        |  9 +++
 .../unpolarized/space_like/as4/__init__.py    | 81 +++++++++++++++++++
 2 files changed, 90 insertions(+)

diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
index 6a3c6cf06..0fbceb805 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
@@ -155,6 +155,11 @@ def gamma_ns_qed(order, mode, n, nf):
             gamma_ns[3, 0] = as3.gamma_nsp(n, nf, cache)
         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)
     return gamma_ns
 
 
@@ -280,6 +285,8 @@ def gamma_singlet_qed(order, n, nf):
         gamma_s[0, 2] = aem2.gamma_singlet(n, nf, cache)
     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)
     return gamma_s
 
 
@@ -314,4 +321,6 @@ def gamma_valence_qed(order, n, nf):
         gamma_v[0, 2] = aem2.gamma_valence(n, nf, cache)
     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)
     return gamma_v
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
index 10c4cb57c..56c2c16ff 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
@@ -55,3 +55,84 @@ 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
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(3,0)}(N)`
+
+    """
+    gamma_np_p = gamma_nsp(N, nf, cache)
+    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):
+    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
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{V}^{(3,0)}(N)`
+
+    """
+    gamma_V = np.array(
+        [
+            [gamma_nsv(N, nf, cache), 0.0],
+            [0.0, gamma_nsm(N, nf, cache)],
+        ],
+        np.complex_,
+    )
+    return gamma_V

From c6ad7dee90cb89b7f841b622e7d2e7f195eba8f6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 15:57:03 +0200
Subject: [PATCH 31/41] Add n3lo_ad_variation key to missing functions

---
 src/eko/evolution_operator/__init__.py                      | 5 ++++-
 .../anomalous_dimensions/unpolarized/space_like/__init__.py | 6 ++++--
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index eab3c64de..27b8933dc 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -450,6 +450,7 @@ def quad_ker_qed(
     ev_op_max_order,
     sv_mode,
     is_threshold,
+    n3lo_ad_variation,
 ):
     """Raw evolution kernel inside quad.
 
@@ -489,6 +490,8 @@ def quad_ker_qed(
         scale variation mode, see `eko.scale_variations.Modes`
     is_threshold : boolean
         is this an itermediate threshold operator?
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -497,7 +500,7 @@ 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)
+        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
index 0fbceb805..6358a5af8 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
@@ -255,7 +255,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache):
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(order, n, nf):
+def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
     r"""
     Compute the grid of the QED singlet anomalous dimensions matrices.
 
@@ -267,6 +267,8 @@ def gamma_singlet_qed(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -286,7 +288,7 @@ def gamma_singlet_qed(order, n, nf):
     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)
+        gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation)
     return gamma_s
 
 

From b7eedfc3b62a5ebf844acd7612e3eb377a13883b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:05:28 +0200
Subject: [PATCH 32/41] Add n3lo_ad_variation key to missing functions. again

---
 src/eko/evolution_operator/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 27b8933dc..122bfbc69 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -307,6 +307,7 @@ def quad_ker(
             ev_op_max_order,
             sv_mode,
             is_threshold,
+            n3lo_ad_variation,
         )
 
     # recombine everything

From 564675382ccc14fa3575065240586cdb0092095c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:14:18 +0200
Subject: [PATCH 33/41] Fix test_kernels_QEDsinglet.py

---
 tests/eko/kernels/test_kernels_QEDsinglet.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/eko/kernels/test_kernels_QEDsinglet.py b/tests/eko/kernels/test_kernels_QEDsinglet.py
index 302bb941b..f5803750c 100644
--- a/tests/eko/kernels/test_kernels_QEDsinglet.py
+++ b/tests/eko/kernels/test_kernels_QEDsinglet.py
@@ -79,7 +79,7 @@ def test_zero_true_gamma(monkeypatch):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()
-            gamma_s = ad.gamma_singlet_qed(order, n, nf)
+            gamma_s = ad.gamma_singlet_qed(order, n, nf, (0, 0, 0, 0))
             # monkeypatch.setattr(
             #     ad,
             #     "exp_matrix_2D",

From b67c9d008b2a5ca2a6b8e3dbac84acbeac51a10b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:30:38 +0200
Subject: [PATCH 34/41] Add order=4 in fixed_alphaem_exact

---
 src/eko/kernels/non_singlet_qed.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/eko/kernels/non_singlet_qed.py b/src/eko/kernels/non_singlet_qed.py
index 86ec1a9ff..36d4d3c35 100644
--- a/src/eko/kernels/non_singlet_qed.py
+++ b/src/eko/kernels/non_singlet_qed.py
@@ -120,6 +120,29 @@ def as3_exact(gamma_ns, a1, a0, beta):
     return ns.nnlo_exact(gamma_ns, a1, a0, beta)
 
 
+@nb.njit(cache=True)
+def as4_exact(gamma_ns, a1, a0, beta):
+    """O(as3aem1) non-singlet exact EKO.
+
+    Parameters
+    ----------
+    gamma_ns : numpy.ndarray
+        non-singlet anomalous dimensions
+    a1 : float
+        target coupling value
+    a0 : float
+        initial coupling value
+    beta : list
+        list of the values of the beta functions
+
+    Returns
+    -------
+    e_ns^3 : complex
+        O(as4aem1) non-singlet exact EKO
+    """
+    return ns.n3lo_exact(gamma_ns, a1, a0, beta)
+
+
 @nb.njit(cache=True)
 def dispatcher(
     order,
@@ -216,6 +239,8 @@ def fixed_alphaem_exact(order, gamma_ns, a1, a0, aem, nf, mu2_from, mu2_to):
         qcd_only = as2_exact(gamma_ns_list[1:], a1, a0, betalist)
     elif order[0] == 3:
         qcd_only = as3_exact(gamma_ns_list[1:], a1, a0, betalist)
+    elif order[0] == 4:
+        qcd_only = as4_exact(gamma_ns_list[1:], a1, a0, betalist)
     else:
         raise NotImplementedError("Selected order is not implemented")
     return qcd_only * apply_qed(gamma_ns_list[0], mu2_from, mu2_to)

From 1f953265a8bae73af8fd1844e0a867d30daacc9c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:45:39 +0200
Subject: [PATCH 35/41] Fix test_kernels_QEDns.py

---
 tests/eko/kernels/test_kernels_QEDns.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index e84956084..5574b1a12 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -304,7 +304,7 @@ def test_error():
     for running in [True, False]:
         with pytest.raises(NotImplementedError):
             ns.dispatcher(
-                (4, 2),
+                (5, 2),
                 "iterate-exact",
                 np.random.rand(4, 3),
                 [0.1, 0.2],

From 607d0ebe0baa9fe040ce268a6da4e2ccf11f70b8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 16:56:16 +0200
Subject: [PATCH 36/41] Fix again test

---
 tests/eko/kernels/test_kernels_QEDns.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index 5574b1a12..f2bd45731 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -302,7 +302,7 @@ def test_ode_true_gamma():
 
 def test_error():
     for running in [True, False]:
-        with pytest.raises(NotImplementedError):
+        with pytest.raises(ValueError):
             ns.dispatcher(
                 (5, 2),
                 "iterate-exact",

From 0009b3f1ac2b37ddefb15cfdf0341efc99d258e1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Fri, 6 Oct 2023 17:09:22 +0200
Subject: [PATCH 37/41] Fix another test

---
 .../anomalous_dimensions/unpolarized/space_like/test_init.py    | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
index 2175eddcd..c7c955960 100644
--- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
@@ -208,7 +208,7 @@ def test_dim_singlet():
     nf = 3
     N = 2
     cache = h.cache.reset()
-    gamma_singlet = ad_us.gamma_singlet_qed((3, 2), N, nf)
+    gamma_singlet = ad_us.gamma_singlet_qed((3, 2), N, nf, (0, 0, 0, 0))
     assert gamma_singlet.shape == (4, 3, 4, 4)
     gamma_singlet_as1 = ad_us.as1.gamma_singlet_qed(N, cache, nf)
     assert gamma_singlet_as1.shape == (4, 4)

From 233e1de3ee16562d429c5bcf98f8c9e463d6f477 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Sat, 7 Oct 2023 19:44:38 +0200
Subject: [PATCH 38/41] Change nf->beta into ns kernels at n3lo

---
 src/eko/kernels/non_singlet.py | 24 ++++++++----------------
 1 file changed, 8 insertions(+), 16 deletions(-)

diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py
index 24706a794..23b2ab47a 100644
--- a/src/eko/kernels/non_singlet.py
+++ b/src/eko/kernels/non_singlet.py
@@ -186,7 +186,7 @@ def nnlo_expanded(gamma_ns, a1, a0, beta):
 
 
 @nb.njit(cache=True)
-def n3lo_expanded(gamma_ns, a1, a0, nf):
+def n3lo_expanded(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet expanded EKO.
 
     Parameters
@@ -206,12 +206,8 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet expanded EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_expanded(a1, a0, beta0, b_list)
     j23 = as4_ei.j23_expanded(a1, a0, beta0, b_list)
@@ -225,7 +221,7 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
 
 
 @nb.njit(cache=True)
-def n3lo_exact(gamma_ns, a1, a0, nf):
+def n3lo_exact(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet exact EKO.
 
     Parameters
@@ -245,12 +241,8 @@ def n3lo_exact(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet exact EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     roots = as4_ei.roots(b_list)
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_exact(a1, a0, beta0, b_list, roots)
@@ -420,6 +412,6 @@ def dispatcher(
             "decompose-expanded",
             "perturbative-expanded",
         ]:
-            return n3lo_expanded(gamma_ns, a1, a0, nf)
-        return n3lo_exact(gamma_ns, a1, a0, nf)
+            return n3lo_expanded(gamma_ns, a1, a0, betalist)
+        return n3lo_exact(gamma_ns, a1, a0, betalist)
     raise NotImplementedError("Selected order is not implemented")

From 69ddbde604e7c3121895028da13293394d53b2bf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Sun, 8 Oct 2023 15:43:31 +0200
Subject: [PATCH 39/41] Add pto = 3 in some qed tests

---
 tests/eko/evolution_operator/test_init.py    |  2 +-
 tests/eko/kernels/test_kernels_QEDns.py      | 14 +++++++++-----
 tests/eko/kernels/test_kernels_QEDsinglet.py |  4 ++--
 tests/eko/kernels/test_kernels_QEDvalence.py |  4 ++--
 4 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index 642de5bb7..47d5e7001 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -382,7 +382,7 @@ def test_compute(self, monkeypatch, theory_ffns, operator_card, tmp_path):
                     err_msg=k,
                 )
 
-        for n in range(1, 3 + 1):
+        for n in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 g.config["order"] = (n, qed)
                 o1 = Operator(g.config, g.managers, Segment(2.0, 2.0, 3))
diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index f2bd45731..e6179ba99 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -25,7 +25,7 @@ def test_zero():
     nf = 3
     ev_op_iterations = 2
     running_alpha = [True, False]
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_ns = (
@@ -73,7 +73,7 @@ def test_zero_true_gamma():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        for qcd in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 n = np.random.rand()
@@ -136,7 +136,8 @@ def test_ode():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        max_qcd = 4
+        for qcd in range(1, max_qcd + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 sc = Couplings(
@@ -154,14 +155,17 @@ def test_ode():
                     )
                     a1 = sc.a_s(mu2_to)
                     gamma_ns = (
-                        np.random.rand(3 + 1, 2 + 1) + np.random.rand(3 + 1, 2 + 1) * 1j
+                        np.random.rand(max_qcd + 1, 2 + 1)
+                        + np.random.rand(max_qcd + 1, 2 + 1) * 1j
                     )
                     gamma_ns[0, 0] = 0.0
                     gamma_ns[2, 1] = 0.0
                     gamma_ns[3, 1] = 0.0
+                    gamma_ns[4, 1] = 0.0
                     gamma_ns[1, 2] = 0.0
                     gamma_ns[2, 2] = 0.0
                     gamma_ns[3, 2] = 0.0
+                    gamma_ns[4, 2] = 0.0
                     gammatot = 0.0
                     for i in range(0, order[0] + 1):
                         for j in range(0, order[1] + 1):
@@ -226,7 +230,7 @@ def test_ode_true_gamma():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        for qcd in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 sc = Couplings(
diff --git a/tests/eko/kernels/test_kernels_QEDsinglet.py b/tests/eko/kernels/test_kernels_QEDsinglet.py
index f5803750c..b953c053c 100644
--- a/tests/eko/kernels/test_kernels_QEDsinglet.py
+++ b/tests/eko/kernels/test_kernels_QEDsinglet.py
@@ -21,7 +21,7 @@ def test_zero(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_s = (
@@ -75,7 +75,7 @@ def test_zero_true_gamma(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()
diff --git a/tests/eko/kernels/test_kernels_QEDvalence.py b/tests/eko/kernels/test_kernels_QEDvalence.py
index 95d0cbc4c..a27bc6ef6 100644
--- a/tests/eko/kernels/test_kernels_QEDvalence.py
+++ b/tests/eko/kernels/test_kernels_QEDvalence.py
@@ -21,7 +21,7 @@ def test_zero(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_v = (
@@ -62,7 +62,7 @@ def test_zero_true_gamma(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()

From ed4dedf67326fe5231316743007de7da4c9777ec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 10 Oct 2023 16:42:26 +0200
Subject: [PATCH 40/41] restructure self.compute_aem_list

---
 src/eko/evolution_operator/__init__.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 122bfbc69..29e67c19c 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -626,7 +626,7 @@ def __init__(
         self.alphaem_running = self.managers["couplings"].alphaem_running
         if self.log_label == "Evolution":
             self.a = self.compute_a()
-            self.compute_aem_list()
+            self.as_list, self.a_half_list = self.compute_aem_list()
 
     @property
     def n_pools(self):
@@ -701,8 +701,8 @@ def compute_aem_list(self):
         """
         ev_op_iterations = self.config["ev_op_iterations"]
         if self.order[1] == 0:
-            self.as_list = np.array([self.a_s[0], self.a_s[1]])
-            self.a_half_list = np.zeros((ev_op_iterations, 2))
+            as_list = np.array([self.a_s[0], self.a_s[1]])
+            a_half = np.zeros((ev_op_iterations, 2))
         else:
             as0 = self.a_s[0]
             as1 = self.a_s[1]
@@ -722,7 +722,7 @@ def compute_aem_list(self):
             couplings = self.managers["couplings"]
             mu2_steps = utils.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations)
             mu2_l = mu2_steps[0]
-            self.as_list = np.array(
+            as_list = np.array(
                 [
                     couplings.compute(
                         a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2
@@ -738,7 +738,7 @@ def compute_aem_list(self):
                 )
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
-            self.a_half_list = a_half
+        return as_list, a_half
 
     @property
     def labels(self):

From b13a26057961eb9dfb4d8d48aad5249c382de290 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 10 Oct 2023 17:41:23 +0200
Subject: [PATCH 41/41] Remove unused default parameters and remove qed
 argument from apply_pdf

---
 src/eko/basis_rotation.py                        |  4 ++--
 src/eko/evolution_operator/flavors.py            |  8 ++++----
 src/eko/evolution_operator/matching_condition.py |  2 +-
 src/eko/evolution_operator/physical.py           |  2 +-
 src/ekobox/apply.py                              |  3 ++-
 tests/eko/evolution_operator/test_flavors.py     | 16 +++++++++-------
 .../test_matching_condition.py                   | 10 +++++-----
 tests/eko/test_basis_rotation.py                 |  8 ++++----
 8 files changed, 28 insertions(+), 25 deletions(-)

diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py
index 4b972592c..e31a7beb6 100644
--- a/src/eko/basis_rotation.py
+++ b/src/eko/basis_rotation.py
@@ -275,7 +275,7 @@
 }
 
 
-def ad_projector(ad_lab, nf, qed=False):
+def ad_projector(ad_lab, nf, qed):
     """
     Build a projector (as a numpy array) for the given anomalous dimension sector.
 
@@ -355,7 +355,7 @@ def select_light_flavors_uni_ev(ad_lab, nf):
         return map_ad_to_evolution[ad_lab]
 
 
-def ad_projectors(nf, qed=False):
+def ad_projectors(nf, qed):
     """
     Build projectors tensor (as a numpy array), collecting all the individual sector projectors.
 
diff --git a/src/eko/evolution_operator/flavors.py b/src/eko/evolution_operator/flavors.py
index 7ad7b5801..fac94a608 100644
--- a/src/eko/evolution_operator/flavors.py
+++ b/src/eko/evolution_operator/flavors.py
@@ -52,7 +52,7 @@ def pids_from_intrinsic_evol(label, nf, normalize):
     return weights
 
 
-def get_range(evol_labels, qed=False):
+def get_range(evol_labels, qed):
     """Determine the number of light and heavy flavors participating in the input and output.
 
     Here, we assume that the T distributions (e.g. T15) appears *always*
@@ -73,7 +73,7 @@ def get_range(evol_labels, qed=False):
     nf_in = 3
     nf_out = 3
 
-    def update(label, qed=False):
+    def update(label, qed):
         nf = 3
         if label[0] == "T":
             if not qed:
@@ -129,7 +129,7 @@ def rotate_pm_to_flavor(label):
     return l
 
 
-def rotate_matching(nf, qed=False, inverse=False):
+def rotate_matching(nf, qed, inverse=False):
     """Rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
@@ -206,7 +206,7 @@ def rotate_matching(nf, qed=False, inverse=False):
     return l
 
 
-def rotate_matching_inverse(nf, qed=False):
+def rotate_matching_inverse(nf, qed):
     """Inverse rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
diff --git a/src/eko/evolution_operator/matching_condition.py b/src/eko/evolution_operator/matching_condition.py
index 059617fbf..8df4775ff 100644
--- a/src/eko/evolution_operator/matching_condition.py
+++ b/src/eko/evolution_operator/matching_condition.py
@@ -21,7 +21,7 @@ def split_ad_to_evol_map(
         nf,
         q2_thr,
         intrinsic_range,
-        qed=False,
+        qed,
     ):
         """
         Create the instance from the |OME|.
diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py
index c963e2a26..1ee501a2a 100644
--- a/src/eko/evolution_operator/physical.py
+++ b/src/eko/evolution_operator/physical.py
@@ -22,7 +22,7 @@ class PhysicalOperator(member.OperatorBase):
     """
 
     @classmethod
-    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False):
+    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed):
         """
         Obtain map between the 3-dimensional anomalous dimension basis and the 4-dimensional evolution basis.
 
diff --git a/src/ekobox/apply.py b/src/ekobox/apply.py
index 354461c0a..336900798 100644
--- a/src/ekobox/apply.py
+++ b/src/ekobox/apply.py
@@ -12,7 +12,6 @@ def apply_pdf(
     lhapdf_like,
     targetgrid=None,
     rotate_to_evolution_basis=False,
-    qed=False,
 ):
     """
     Apply all available operators to the input PDFs.
@@ -34,6 +33,7 @@ def apply_pdf(
         out_grid : dict
             output PDFs and their associated errors for the computed mu2grid
     """
+    qed = eko.theory_card.order[1] > 0
     if rotate_to_evolution_basis:
         if not qed:
             rotate_flavor_to_evolution = br.rotate_flavor_to_evolution
@@ -99,6 +99,7 @@ def apply_pdf_flavor(
         if error_final is not None:
             out_grid[ep]["errors"] = dict(zip(eko.bases.targetpids, error_final))
 
+    qed = eko.theory_card.order[1] > 0
     # rotate to evolution basis
     if flavor_rotation is not None:
         for q2, op in out_grid.items():
diff --git a/tests/eko/evolution_operator/test_flavors.py b/tests/eko/evolution_operator/test_flavors.py
index 81bb35fc4..2f3761404 100644
--- a/tests/eko/evolution_operator/test_flavors.py
+++ b/tests/eko/evolution_operator/test_flavors.py
@@ -61,13 +61,15 @@ def get(d):
 
 
 def test_get_range():
-    assert (3, 3) == flavors.get_range([])
-    assert (3, 3) == flavors.get_range([member.MemberName(n) for n in ["S.S", "V3.V3"]])
+    assert (3, 3) == flavors.get_range([], False)
+    assert (3, 3) == flavors.get_range(
+        [member.MemberName(n) for n in ["S.S", "V3.V3"]], False
+    )
     assert (3, 4) == flavors.get_range(
-        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.S"]]
+        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.S"]], False
     )
     assert (4, 4) == flavors.get_range(
-        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.T15"]]
+        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.T15"]], False
     )
     assert (3, 3) == flavors.get_range(
         [member.MemberName(n) for n in ["S.S", "Td3.Td3"]], True
@@ -113,7 +115,7 @@ def test_rotate_pm_to_flavor():
 
 
 def test_rotate_matching():
-    m = flavors.rotate_matching(4)
+    m = flavors.rotate_matching(4, False)
     assert len(list(filter(lambda e: "c+" in e, m.keys()))) == 2
     assert len(list(filter(lambda e: "b-" in e, m.keys()))) == 1
 
@@ -156,8 +158,8 @@ def load(m):
         return mm
 
     for nf in range(4, 6 + 1):
-        m = load(flavors.rotate_matching(nf))
-        minv = load(flavors.rotate_matching_inverse(nf))
+        m = load(flavors.rotate_matching(nf, False))
+        minv = load(flavors.rotate_matching_inverse(nf, False))
         np.testing.assert_allclose(m @ minv, np.eye(len(br.evol_basis)), atol=1e-10)
 
 
diff --git a/tests/eko/evolution_operator/test_matching_condition.py b/tests/eko/evolution_operator/test_matching_condition.py
index 8cbd6b537..c854fb90c 100644
--- a/tests/eko/evolution_operator/test_matching_condition.py
+++ b/tests/eko/evolution_operator/test_matching_condition.py
@@ -38,7 +38,7 @@ def update_intrinsic_OME(self, ome):
 
     def test_split_ad_to_evol_map(self):
         ome = self.mkOME()
-        a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [])
+        a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False)
         triv_keys = [
             "V.V",
             "T3.T3",
@@ -64,7 +64,7 @@ def test_split_ad_to_evol_map(self):
             ome[(200, 200)].value,
         )
         # # if alpha is zero, nothing non-trivial should happen
-        b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [])
+        b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False)
         assert sorted(str(k) for k in b.op_members.keys()) == sorted(
             [*triv_keys, *keys3]
         )
@@ -74,7 +74,7 @@ def test_split_ad_to_evol_map(self):
         # )
         # nf=3 + IC
         self.update_intrinsic_OME(ome)
-        c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4])
+        c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4], False)
         assert sorted(str(k) for k in c.op_members.keys()) == sorted(
             [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"]
         )
@@ -83,7 +83,7 @@ def test_split_ad_to_evol_map(self):
             b.op_members[member.MemberName("V.V")].value,
         )
         # nf=3 + IB
-        d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5])
+        d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5], False)
         assert sorted(str(k) for k in d.op_members.keys()) == sorted(
             [*triv_keys, *keys3, "b+.b+", "b-.b-"]
         )
@@ -92,7 +92,7 @@ def test_split_ad_to_evol_map(self):
             np.eye(self.shape[0]),
         )
         # nf=4 + IB
-        d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5])
+        d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5], False)
         assert sorted(str(k) for k in d.op_members.keys()) == sorted(
             [
                 *triv_keys,
diff --git a/tests/eko/test_basis_rotation.py b/tests/eko/test_basis_rotation.py
index 680ef065d..194669d46 100644
--- a/tests/eko/test_basis_rotation.py
+++ b/tests/eko/test_basis_rotation.py
@@ -9,19 +9,19 @@ def test_ad_projector():
     g = br.rotate_flavor_to_evolution[2]
     v3 = br.rotate_flavor_to_evolution[br.evol_basis.index("V3")]
 
-    s_to_s = br.ad_projector((100, 100), nf=6)
+    s_to_s = br.ad_projector((100, 100), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ s_to_s, s)
     np.testing.assert_allclose(g @ s_to_s, 0.0)
     np.testing.assert_allclose(v3 @ s_to_s, 0.0)
 
-    g_to_s = br.ad_projector((21, 100), nf=6)
+    g_to_s = br.ad_projector((21, 100), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ g_to_s, 0.0)
     np.testing.assert_allclose(g @ g_to_s, s)
     np.testing.assert_allclose(v3 @ g_to_s, 0.0)
 
-    ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6)
+    ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15)
     np.testing.assert_allclose(g @ ns_m, 0.0)
@@ -74,7 +74,7 @@ def test_ad_projectors():
     for nf in range(3, 6 + 1):
         diag = np.array([0] * (1 + 6 - nf) + [1] * (1 + 2 * nf) + [0] * (6 - nf))
         identity = np.diag(diag)
-        projs = br.ad_projectors(nf)
+        projs = br.ad_projectors(nf, qed=False)
 
         # sum over diagonal projectors form an identity
         np.testing.assert_allclose(