Merge pull request #391 from NNPDF/alphas_by_nf_info

Use nf to generate alpha_s

src/ekobox/evol_pdf.py

@@ -1,13 +1,15 @@
 """Tools to evolve actual PDFs."""
 
 import pathlib
-from collections import defaultdict
+
+import numpy as np
 
 from eko import basis_rotation as br
 from eko.io import EKO
 from eko.runner import managed
 
 from . import apply, genpdf, info_file
+from .utils import regroup_evolgrid
 
 DEFAULT_NAME = "eko.tar"
 
@@ -46,6 +48,18 @@ def evolve_pdfs(
     info_update : dict
         dict of info to add or update to default info file
     """
+    # separate by nf the evolgrid (and order per nf/q)
+    q2block_per_nf = regroup_evolgrid(operators_card.mugrid)
+
+    # check we have disjoint scale ranges
+    nfs = list(q2block_per_nf.keys())
+    for j in range(len(nfs) - 1):
+        # equal points are allowed by LHAPDF
+        if q2block_per_nf[nfs[j]][-1] > q2block_per_nf[nfs[j + 1]][0]:
+            raise ValueError(
+                f"Last scale point for nf={nfs[j]} is bigger than first in nf={nfs[j+1]}"
+            )
+
     # update op and th cards
     if path is not None:
         eko_path = pathlib.Path(path)
@@ -79,9 +93,8 @@ def evolve_pdfs(
         info_update=info_update,
     )
 
-    # separate by nf the evolgrid (and order per nf/q)
-    q2block_per_nf = regroup_evolgrid(eko_output.evolgrid)
-
+    # in the eko scales are squared
+    q2block_per_nf = {nf: np.power(q2s, 2) for nf, q2s in q2block_per_nf.items()}
     # write all replicas
     all_member_blocks = []
     for evolved_PDF in evolved_PDF_list:
@@ -95,14 +108,6 @@ def evolve_pdfs(
         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.
 

src/ekobox/info_file.py

@@ -9,26 +9,27 @@
 from eko.io.runcards import OperatorCard, TheoryCard
 
 from .genpdf import load
+from .utils import regroup_evolgrid
 
 
 def build(
     theory_card: TheoryCard,
     operators_card: OperatorCard,
     num_members: int,
     info_update: dict,
-):
-    """Generate a lhapdf info file from theory and operators card.
+) -> dict:
+    """Generate a lhapdf info file.
 
     Parameters
     ----------
-    theory_card : dict
+    theory_card :
         theory card
-    operators_card : dict
-        operators_card
-    num_members : int
+    operators_card :
+        operators card
+    num_members :
         number of pdf set members
-    info_update : dict
-        info to update
+    info_update :
+        additional info to update
 
     Returns
     -------
@@ -56,8 +57,35 @@ def build(
     template_info["MCharm"] = theory_card.heavy.masses.c.value
     template_info["MBottom"] = theory_card.heavy.masses.b.value
     template_info["MTop"] = theory_card.heavy.masses.t.value
+    # dump alphas
+    template_info.update(build_alphas(theory_card, operators_card))
+    return template_info
+
+
+def build_alphas(
+    theory_card: TheoryCard,
+    operators_card: OperatorCard,
+) -> dict:
+    """Generate a couplings section of lhapdf info file.
+
+    Parameters
+    ----------
+    theory_card : dict
+        theory card
+    operators_card : dict
+        operators card
+
+    Returns
+    -------
+    dict
+        info file section in lhapdf format
+    """
+    # start with meta stuff
+    template_info = {}
     template_info["AlphaS_MZ"] = theory_card.couplings.alphas
     template_info["AlphaS_OrderQCD"] = theory_card.order[0] - 1
+    # prepare
+    evolgrid = regroup_evolgrid(operators_card.mugrid)
     evmod = couplings.couplings_mod_ev(operators_card.configs.evolution_method)
     quark_masses = [(x.value) ** 2 for x in theory_card.heavy.masses]
     sc = couplings.Couplings(
@@ -68,19 +96,14 @@ def build(
         hqm_scheme=theory_card.heavy.masses_scheme,
         thresholds_ratios=np.power(list(iter(theory_card.heavy.matching_ratios)), 2.0),
     )
-    alphas_values = np.array(
-        [
-            4.0
-            * np.pi
-            * sc.a_s(
-                muf2,
-            )
-            for muf2 in operators_card.mu2grid
-        ],
-        dtype=float,
-    )
-    template_info["AlphaS_Vals"] = alphas_values.tolist()
-    template_info["AlphaS_Qs"] = np.array(
-        [mu for mu, _ in operators_card.mugrid]
-    ).tolist()
+    # add actual values
+    alphas_values = []
+    alphas_qs = []
+    for nf, mus in evolgrid.items():
+        for mu in mus:
+            alphas_values.append(float(4.0 * np.pi * sc.a_s(mu * mu, nf_to=nf)))
+            alphas_qs.append(mu)
+
+    template_info["AlphaS_Vals"] = alphas_values
+    template_info["AlphaS_Qs"] = alphas_qs
     return template_info

src/ekobox/utils.py

@@ -1,6 +1,7 @@
 """Generic utilities to work with EKOs."""
 
 import os
+from collections import defaultdict
 from typing import Optional
 
 import numpy as np
@@ -78,3 +79,11 @@ def ekos_product(
 
     if path is not None:
         final_eko.close()
+
+
+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()}

tests/ekobox/test_evol_pdf.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 from banana import toy
 
 import eko
@@ -50,6 +51,21 @@ def test_evolve_pdfs_dump_path(fake_lhapdf, cd):
     assert p.exists()
 
 
+def test_evolve_pdfs_bad_scales(fake_lhapdf, cd):
+    """Bad scale configurations."""
+    theory, op = init_cards()
+    n = "test_evolve_pdfs_bad_scales"
+    op = cards.example.operator()
+    op.mugrid = [(5.0, 3), (15.0, 4), (10.0, 5), (100.0, 5)]
+    with pytest.raises(ValueError, match="is bigger"):
+        with cd(fake_lhapdf):
+            ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=fake_lhapdf)
+    op.mugrid = [(5.0, 3), (10.0, 3), (10.0, 4), (15.0, 4), (10.0, 5)]
+    with pytest.raises(ValueError, match="is bigger"):
+        with cd(fake_lhapdf):
+            ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=fake_lhapdf)
+
+
 def test_evolve_pdfs_dump_file(fake_lhapdf, cd):
     theory, op = init_cards()
     n = "test_evolve_pdfs_dump_file"

tests/ekobox/test_info_file.py

@@ -23,3 +23,27 @@ def test_build():
     np.testing.assert_allclose(info["QMin"], math.sqrt(op.mu2grid[0]), rtol=1e-5)
     assert info["XMin"] == op.xgrid.raw[0]
     assert info["XMax"] == op.xgrid.raw[-1] == 1.0
+
+
+def test_build_alphas_good():
+    """Good configurations."""
+    theory = cards.example.theory()
+    theory.order = (2, 0)
+    theory.couplings.alphas = 0.2
+    op = cards.example.operator()
+    op.mu0 = 1.0
+    # base case
+    op.mugrid = [(100.0, 5), (10.0, 5)]
+    info = info_file.build_alphas(theory, op)
+    assert len(info["AlphaS_Vals"]) == 2
+    np.testing.assert_allclose(info["AlphaS_Qs"], [10.0, 100.0])
+    # several nf
+    op.mugrid = [(5.0, 4), (10.0, 5), (1.0, 3), (5.0, 3), (10.0, 5), (100.0, 5)]
+    info = info_file.build_alphas(theory, op)
+    assert len(info["AlphaS_Vals"]) == 6
+    np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 5.0, 5.0, 10.0, 10.0, 100.0])
+    # several nf with gap
+    op.mugrid = [(1.0, 3), (10.0, 3), (10.0, 5), (100.0, 5)]
+    info = info_file.build_alphas(theory, op)
+    assert len(info["AlphaS_Vals"]) == 4
+    np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 10.0, 10.0, 100.0])