Interface QED and N3LO

.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

benchmarks/eko/benchmark_inverse_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 benchmark_inverse_matching():
+    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]
+        )

extras/lh_bench_23/.gitignore

@@ -0,0 +1,2 @@
+*.tar
+*.csv

extras/lh_bench_23/cfg.py

@@ -0,0 +1,145 @@
+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),
+    matching_order=[2, 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

extras/lh_bench_23/run-nnlo.py

@@ -0,0 +1,118 @@
+import argparse
+import logging
+import pathlib
+import sys
+
+import numpy as np
+import pandas as pd
+import yaml
+from banana import toy
+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.runner.managed import solve
+from ekobox import apply
+from ekomark.benchmark.external.LHA_utils import here as there
+
+_sqrt2 = float(np.sqrt(2))
+
+
+# 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])
+
+
+# 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("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
+        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")
+
+    # eko path
+    p = pathlib.Path(f"{scheme}-{sv}.tar")
+
+    # recompute?
+    if not p.exists() or args.rerun:
+        print("(Re)running eko ...")
+        p.unlink(True)
+        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 = {}
+    with eko.EKO.read(p) as eko_:
+        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
+    pd.set_option("display.float_format", "{:.4e}".format)
+    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}"])
+    print()
+    print("rel. distance to reference")
+    print((me - ref) / ref)