Fix apfel benchmarks

benchmarks/CT14_bench.py

@@ -5,8 +5,6 @@
   whereas the lo family uses LO splitting functions with NLO alphas evolution
 """
 
-from math import nan
-
 from banana import register
 
 from ekomark.benchmark.runner import Runner
@@ -16,7 +14,7 @@
 
 base_theory = {
     "Qref": 91.1876,
-    "QrefQED": nan,
+    "Qedref": 0.0,
     "mc": 1.3,
     "mb": 4.75,
     "mt": 172,

benchmarks/CT18_bench.py

@@ -12,7 +12,7 @@
 
 base_theory = {
     "Qref": 91.1870,
-    "QrefQED": 91.1870,
+    "Qedref": 91.1870,
     "mc": 1.3,
     "mb": 4.75,
     "mt": 172.0,

benchmarks/HERA20_bench.py

@@ -2,7 +2,6 @@
 Benchmark HERAPDF2.0 pdf family
 
 """
-from math import nan
 
 from banana import register
 
@@ -14,7 +13,7 @@
 
 base_theory = {
     "Qref": 91.1876,
-    "QrefQED": nan,
+    "Qedref": 0.0,
     "mc": 1.43,
     "mb": 4.5,
     "mt": 173.0,

benchmarks/NNPDF_bench.py

@@ -72,7 +72,7 @@ def benchmark_nnlo(self, Q0=1.65, mugrid=(10,)):
             "QED": 2,
             "Q0": Q0,
         }
-        theory_card.update({"ModEv": "iterate-exact", "FNS": "VFNS", "QrefQED": 91.2})
+        theory_card.update({"ModEv": "iterate-exact", "FNS": "VFNS", "Qedref": 91.2})
 
         self.skip_pdfs = lambda _theory: [
             -6,

benchmarks/apfel_bench.py

@@ -202,7 +202,7 @@ class BenchmarkFFNS_qed(ApfelBenchmark):
 
     ffns_theory = {
         "Qref": 91.1870,
-        "QrefQED": 91.1870,
+        "Qedref": 91.1870,
         "mc": 1.3,
         "mb": 4.75,
         "mt": 172.0,
@@ -281,7 +281,7 @@ class BenchmarkVFNS_qed(ApfelBenchmark):
 
     vfns_theory = {
         "Qref": 91.1870,
-        "QrefQED": 91.1870,
+        "Qedref": 91.1870,
         "mc": 1.3,
         "mb": 4.75,
         "mt": 172.0,

benchmarks/lha_paper_bench.py

@@ -2,7 +2,6 @@
 Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO).
 """
 import os
-from math import nan
 
 import numpy as np
 import pytest
@@ -29,7 +28,7 @@
     "alphas": 0.35,  # Eq. (4.55) :cite:`Dittmar:2005ed`
     "alphaqed": 0.007496,
     "QED": 0,
-    "QrefQED": nan,
+    "Qedref": 0.0,
 }
 """Global theory settings"""
 

benchmarks/pegasus_bench.py

@@ -1,7 +1,6 @@
 """
 Benchmark to Pegasus :cite:`Vogt:2004ns`
 """
-from math import nan
 
 import numpy as np
 from banana import register
@@ -56,7 +55,7 @@ class BenchmarkVFNS(PegasusBenchmark):
         "kbThr": 1.0,
         "ktThr": 1.0,
         "Qref": np.sqrt(2.0),
-        "QrefQED": nan,
+        "Qedref": 0.0,
         "alphas": 0.35,
         "alphaqed": 0.007496,
         "QED": 0,
@@ -110,7 +109,7 @@ class BenchmarkFFNS(PegasusBenchmark):
         "kbThr": np.inf,
         "ktThr": np.inf,
         "Qref": np.sqrt(2.0),
-        "QrefQED": nan,
+        "Qedref": 0.0,
         "alphas": 0.35,
         "alphaqed": 0.007496,
         "Q0": np.sqrt(2.0),

benchmarks/sandbox.py

@@ -1,5 +1,4 @@
 # pylint: skip-file
-from math import nan
 
 import numpy as np
 from banana import register
@@ -27,7 +26,7 @@
 }
 nnpdf_base_theory = {
     "Qref": 91.2,
-    "QrefQED": nan,
+    "Qedref": 0.0,
     "mc": 1.51,
     "mb": 4.92,
     "mt": 172.5,

doc/source/code/IO.rst

@@ -69,7 +69,7 @@ and environment. The benchmark settings are available at :mod:`banana.data.theor
   * - ``alphaqed``
     - :py:obj:`float`
     - Reference value of the electromagnetic coupling :math:`\alpha_{em}`.
-  * - ``QrefQED``
+  * - ``Qedref``
     - :py:obj:`float`
     - Reference scale at which the ``alphaqed`` value is given (in GeV).
   * - ``HQ``

pyproject.toml

@@ -44,7 +44,7 @@ PyYAML = "^6.0"
 lz4 = "^4.0.2"
 numba = "^0.57.0"
 # ekomark
-banana-hep = { version = "^0.6.8", optional = true }
+banana-hep = { version = "^0.6.9", optional = true }
 sqlalchemy = { version = "^1.4.21", optional = true }
 pandas = { version = "^1.3.0", optional = true }
 matplotlib = { version = "^3.5.1", optional = true }

src/eko/io/runcards.py

@@ -173,11 +173,16 @@ def new_theory(self):
 
         new["order"] = [old["PTO"] + 1, old["QED"]]
         alphaem = self.fallback(old.get("alphaqed"), old.get("alphaem"), default=0.0)
+        if "Qedref" in old:
+            em_running = bool(np.isclose(old["Qedref"], old["Qref"]))
+        else:
+            em_running = False
         ms = self.heavies("m%s", old)
         ks = self.heavies("k%sThr", old)
         new["couplings"] = dict(
             alphas=old["alphas"],
             alphaem=alphaem,
+            em_running=em_running,
             scale=old["Qref"],
             num_flavs_ref=old["nfref"],
             max_num_flavs=old["MaxNfAs"],
@@ -220,7 +225,7 @@ def new_operator(self):
             mugrid = old["mugrid"]
         else:
             mu2grid = old["Q2grid"] if "Q2grid" in old else old["mu2grid"]
-            mugrid = np.sqrt(mu2grid)
+            mugrid = np.sqrt(mu2grid).tolist()
         new["mugrid"] = flavored_mugrid(
             mugrid,
             list(self.heavies("m%s", old_th)),

src/ekobox/apply.py

@@ -84,18 +84,18 @@ def apply_pdf_flavor(
 
     # build output
     out_grid = {}
-    for (mu2, nf), elem in eko.items():
+    for ep, elem in eko.items():
         pdf_final = np.einsum(CONTRACTION, elem.operator, pdfs, optimize="optimal")
         if elem.error is not None:
             error_final = np.einsum(CONTRACTION, elem.error, pdfs, optimize="optimal")
         else:
             error_final = None
-        out_grid[mu2] = {
+        out_grid[ep] = {
             "pdfs": dict(zip(eko.bases.targetpids, pdf_final)),
             "errors": None,
         }
         if error_final is not None:
-            out_grid[mu2]["errors"] = dict(zip(eko.bases.targetpids, error_final))
+            out_grid[ep]["errors"] = dict(zip(eko.bases.targetpids, error_final))
 
     # rotate to evolution basis
     if flavor_rotation is not None:

src/ekobox/evol_pdf.py

@@ -80,8 +80,11 @@ def evolve_pdfs(
     targetlist = targetgrid.raw.tolist()
     for evolved_PDF in evolved_PDF_list:
         all_blocks = []
+        evolved_PDF_q2 = {q2: val for (q2, _), val in evolved_PDF.items()}
         block = genpdf.generate_block(
-            lambda pid, x, Q2, evolved_PDF=evolved_PDF: targetlist[targetlist.index(x)]
+            lambda pid, x, Q2, evolved_PDF=evolved_PDF_q2: targetlist[
+                targetlist.index(x)
+            ]
             * evolved_PDF[Q2]["pdfs"][pid][targetlist.index(x)],
             xgrid=targetlist,
             evolgrid=operators_card.evolgrid,

src/ekomark/benchmark/external/pegasus_utils.py

@@ -82,12 +82,12 @@ def compute_pegasus_data(theory, operators, skip_pdfs, rotate_to_evolution_basis
 
     # run pegaus
     out_tabs = {}
-    for mu2 in np.array(operators["mugrid"]) ** 2:
+    for mu in operators["mugrid"]:
         tab = {}
         for x in target_xgrid:
             # last two numbers are the min and max pid to calculate,
             # keep everthing for simplicity.
-            xf, _ = pegasus.xparton(x, mu2, -6, 6)
+            xf, _ = pegasus.xparton(x, mu**2, -6, 6)
             temp = dict(zip(labels, xf))
             for pid in labels:
                 if pid in skip_pdfs:
@@ -114,7 +114,7 @@ def compute_pegasus_data(theory, operators, skip_pdfs, rotate_to_evolution_basis
             evol_pdf = rotate_flavor_to_evolution @ pdfs
             tab = dict(zip(evol_basis, evol_pdf))
 
-        out_tabs[mu2] = tab
+        out_tabs[mu**2] = tab
 
     ref = {"target_xgrid": target_xgrid, "values": out_tabs}
 

src/ekomark/benchmark/runner.py

@@ -4,6 +4,7 @@
 import os
 import sys
 
+import numpy as np
 import pandas as pd
 from banana import cfg as banana_cfg
 from banana.benchmark.runner import BenchmarkRunner
@@ -219,7 +220,6 @@ def log(self, theory, _, pdf, me, ext):
         # return a proper log table
         log_tabs = {}
         xgrid = ext["target_xgrid"]
-        q2s = list(ext["values"].keys())
 
         # LHA NNLO VFNS needs a special treatment
         # Valence contains only u and d
@@ -245,10 +245,15 @@ def log(self, theory, _, pdf, me, ext):
                 flavor_rotation=rotate_to_evolution,
                 qed=qed,
             )
-        for q2 in q2s:
+        for q2, ref_pdfs in ext["values"].items():
             log_tab = dfdict.DFdict()
-            ref_pdfs = ext["values"][q2]
-            res = pdf_grid[q2]
+            # find the first hit, regardless of the flavor, since others can not deal with them appearing in parallel
+            eps = [ep for ep in pdf_grid.keys() if np.isclose(ep[0], q2)]
+            if len(eps) != 1:
+                if len(eps) == 0:
+                    raise KeyError(f"PDF at Q2={q2} not found")
+                raise KeyError(f"More than one evolution point found for Q2={q2}")
+            res = pdf_grid[eps[0]]
             my_pdfs = res["pdfs"]
             my_pdf_errs = res["errors"]
 

tests/eko/io/test_runcards.py

@@ -47,6 +47,7 @@ def test_runcards_opcard():
     np.testing.assert_allclose(no.mu0, tc["Q0"])
     np.testing.assert_allclose(no.mu20, tc["Q0"] ** 2.0)
     assert len(no.pids) == 14
+    check_dumpable(no)
     del oc["mugrid"]
     # or mu2grid
     oc["mu2grid"] = [9.0, 30.0, 32.0]
@@ -57,6 +58,7 @@ def test_runcards_opcard():
     assert no.evolgrid[0][-1] == 4
     assert no.evolgrid[1][-1] == 5
     assert no.evolgrid[2][-1] == 5
+    check_dumpable(no)
     del oc["mu2grid"]
     # or Q2grid
     oc["Q2grid"] = [15.0, 130.0, 140.0, 1e5]
@@ -68,6 +70,15 @@ def test_runcards_opcard():
     assert no.evolgrid[1][-1] == 5
     assert no.evolgrid[2][-1] == 5
     assert no.evolgrid[3][-1] == 6
+    check_dumpable(no)
+
+
+def check_dumpable(no):
+    """Check we can write and read to yaml."""
+    so = StringIO()
+    yaml.dump(no.raw, so)
+    so.seek(0)
+    noo = yaml.safe_load(so)
 
 
 def test_runcards_ekomark():

tests/ekobox/test_apply.py

@@ -11,13 +11,13 @@ def test_apply(self, eko_factory: EKOFactory, fake_pdf):
         # fake pdfs
         pdf_grid = apply.apply_pdf(eko, fake_pdf)
         assert len(pdf_grid) == len(eko.evolgrid)
-        pdfs = pdf_grid[ep_out[0]]["pdfs"]
+        pdfs = pdf_grid[ep_out]["pdfs"]
         assert list(pdfs.keys()) == list(eko.bases.targetpids)
         # rotate to target_grid
         target_grid = [0.75]
         pdf_grid = apply.apply_pdf(eko, fake_pdf, target_grid)
         assert len(pdf_grid) == 1
-        pdfs = pdf_grid[ep_out[0]]["pdfs"]
+        pdfs = pdf_grid[ep_out]["pdfs"]
         assert list(pdfs.keys()) == list(eko.bases.targetpids)
 
     def test_apply_flavor(self, eko_factory: EKOFactory, fake_pdf, monkeypatch):
@@ -37,5 +37,5 @@ def test_apply_flavor(self, eko_factory: EKOFactory, fake_pdf, monkeypatch):
         monkeypatch.setattr("eko.basis_rotation.evol_basis", fake_evol_basis)
         pdf_grid = apply.apply_pdf(eko, fake_pdf, rotate_to_evolution_basis=True)
         assert len(pdf_grid) == len(eko.evolgrid)
-        pdfs = pdf_grid[ep_out[0]]["pdfs"]
+        pdfs = pdf_grid[ep_out]["pdfs"]
         assert list(pdfs.keys()) == list(fake_evol_basis)