Merge pull request #2028 from NNPDF/single_truth

Provide an unique truth for all theory cards

n3fit/src/evolven3fit/eko_utils.py

@@ -7,7 +7,6 @@
 from eko.io import runcards
 from eko.matchings import Atlas, nf_default
 from eko.quantities.heavy_quarks import MatchingScales
-from validphys.loader import Loader
 
 from . import utils
 
@@ -22,7 +21,6 @@
 }
 
 EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = {
-    "ev_op_iterations": 30,
     "ev_op_max_order": (1, 0),
     "evolution_method": "iterate-exact",
     "inversion_method": "exact",
@@ -33,7 +31,7 @@
 
 
 def construct_eko_cards(
-    theoryID,
+    nnpdf_theory,
     q_fin,
     q_points,
     x_grid,
@@ -43,21 +41,21 @@ def construct_eko_cards(
 ):
     """
     Return the theory and operator cards used to construct the eko.
-    theoryID is the ID of the theory for which we are computing the theory and operator card.
+    nnpdf_theory is a NNPDF theory card for which we are computing the operator card and eko
     q_fin is the final point of the q grid while q_points is the number of points of the grid.
     x_grid is the x grid to be used.
     op_card_dict and theory_card_dict are optional updates that can be provided respectively to the
     operator card and to the theory card.
     """
-    theory, thresholds = load_theory(theoryID, theory_card_dict)
+    theory, thresholds = load_theory(nnpdf_theory, theory_card_dict)
 
     mu0 = theory["Q0"]
 
-    # Set nf_0 according to the fitting scale unless set explicitly
-    if "nf0" not in theory:
-        theory["nf0"] = find_nf(mu0, theory, thresholds)
+    # eko needs a value for Qedref and for max nf alphas
+    theory["Qedref"] = theory["Qref"]
+    theory["MaxNfAs"] = theory["MaxNfPdf"]
 
-    # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory
+    # The Legacy function is able to construct a theory card for eko starting from a NNPDF theory
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
 
@@ -106,7 +104,7 @@ def construct_eko_cards(
 
 
 def construct_eko_photon_cards(
-    theoryID,
+    nnpdf_theory,
     q_fin,
     x_grid,
     q_gamma,
@@ -115,28 +113,24 @@ def construct_eko_photon_cards(
 ):
     """
     Return the theory and operator cards used to construct the eko_photon.
-    theoryID is the ID of the theory for which we are computing the theory and operator card.
+    nnpdf_theory is a NNPDF theory card for which we are computing the operator card and eko
     q_fin is the final point of the q grid while q_points is the number of points of the grid.
     x_grid is the x grid to be used.
     op_card_dict and theory_card_dict are optional updates that can be provided respectively to the
     operator card and to the theory card.
     """
-    theory, thresholds = load_theory(theoryID, theory_card_dict)
+    theory, _ = load_theory(nnpdf_theory, theory_card_dict)
 
     # if is eko_photon then mu0 = q_gamma
     mu0 = float(q_gamma)
 
-    # Set nf_0 according to mu0 unless set explicitly
-    if "nf0" not in theory:
-        theory["nf0"] = find_nf(mu0, theory, thresholds)
-
-    # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory
+    # The Legacy function is able to construct a theory card for eko starting from a NNPDF theory
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
 
-    q_fin = float(theory["Q0"])
-
-    nf_fin = find_nf(q_fin, theory, thresholds)
+    # The photon needs to be evolved down to Q0
+    q_fin = theory["Q0"]
+    nf_fin = theory["nf0"]
 
     # construct mugrid
     mugrid = [(q_fin, nf_fin)]
@@ -147,12 +141,12 @@ def construct_eko_photon_cards(
     return theory_card, op_card
 
 
-def load_theory(theoryID, theory_card_dict):
+def load_theory(nnpdf_theory, theory_card_dict):
     """loads and returns the theory dictionary and the thresholds"""
     if theory_card_dict is None:
         theory_card_dict = {}
     # theory_card construction
-    theory = Loader().check_theoryID(theoryID).get_description()
+    theory = dict(nnpdf_theory)
     theory.pop("FNS")
     theory.update(theory_card_dict)
 
@@ -173,7 +167,9 @@ def load_theory(theoryID, theory_card_dict):
 
 
 def build_opcard(op_card_dict, theory, x_grid, mu0, mugrid):
-    """builds the opcard"""
+    """Build the operator card.
+    The user provided options should be given as part of ``op_card_dict``
+    """
     if op_card_dict is None:
         op_card_dict = {}
 
@@ -182,12 +178,17 @@ def build_opcard(op_card_dict, theory, x_grid, mu0, mugrid):
     op_card.update({"mu0": mu0, "mugrid": mugrid})
 
     op_card["xgrid"] = x_grid
-    # Specific defaults for evolven3fit evolution
-    if theory["ModEv"] == "TRN":
-        op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN)
-    if theory["ModEv"] == "EXA":
-        op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA)
-    # User can still change the configs via op_card_dict
+
+    # Specify the evolution options and defaults differently from TRN / EXA
+    configs = op_card["configs"]
+    if theory.get("ModEv") == "TRN":
+        configs.update(EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN)
+    elif theory.get("ModEv") == "EXA":
+        # Set the default from the theory card unless it was given in the input
+        op_card_dict.setdefault("ev_op_iterations", theory.get("IterEv"))
+
+        configs.update(EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA)
+        configs["ev_op_iterations"] = op_card_dict["ev_op_iterations"]
 
     # Note that every entry that is not a dictionary should not be
     # touched by the user and indeed an user cannot touch them
@@ -197,30 +198,5 @@ def build_opcard(op_card_dict, theory, x_grid, mu0, mugrid):
         elif key in op_card_dict:
             _logger.warning("Entry %s is not a dictionary and will be ignored", key)
 
-    # if no -e was given, take ev_op_iterations from EVOLVEN3FIT_CONFIGS_DEFAULTS_{TRN,EXA}
-    if op_card['configs']['ev_op_iterations'] is None:
-        if theory["ModEv"] == "TRN":
-            op_card['configs']['ev_op_iterations'] = EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN[
-                "ev_op_iterations"
-            ]
-        if theory["ModEv"] == "EXA":
-            op_card['configs']['ev_op_iterations'] = EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA[
-                "ev_op_iterations"
-            ]
-
     op_card = runcards.OperatorCard.from_dict(op_card)
-
     return op_card
-
-
-def find_nf(mu, theory, thresholds):
-    """compute nf for a given mu"""
-    if mu < theory["mc"] * thresholds["c"]:
-        nf = 3
-    elif mu < theory["mb"] * thresholds["b"]:
-        nf = 4
-    elif mu < theory["mt"] * thresholds["t"]:
-        nf = 5
-    else:
-        nf = 6
-    return nf

n3fit/src/n3fit/scripts/evolven3fit.py

@@ -1,6 +1,7 @@
 """
 This module contains the CLI for evolven3fit
 """
+
 from argparse import ArgumentParser
 import logging
 import pathlib
@@ -11,6 +12,7 @@
 
 from eko.runner.managed import solve
 from n3fit.io.writer import XGRID
+from validphys.loader import FallbackLoader
 
 _logger = logging.getLogger(__name__)
 
@@ -118,7 +120,7 @@ def main():
         "--ev-op-iterations",
         type=int,
         default=None,
-        help="ev_op_iterations for the EXA theory",
+        help="ev_op_iterations for the EXA theory. Overrides the settings given in the theory card.",
     )
     parser.add_argument(
         "--use-fhmruvv",
@@ -156,6 +158,12 @@ def main():
             args.force,
         )
     else:
+        # If we are in the business of producing an eko, do some checks before starting:
+        # 1. load the nnpdf theory early to check for inconsistent options and theory problems
+        nnpdf_theory = FallbackLoader().check_theoryID(args.theoryID).get_description()
+        if nnpdf_theory.get("ModEv") == "TRN" and args.ev_op_iterations is not None:
+            raise ValueError("ev_op_iterations is not accepted with ModEv=TRN solution")
+
         stdout_log = logging.StreamHandler(sys.stdout)
         stdout_log.setLevel(evolve.LOGGING_SETTINGS["level"])
         stdout_log.setFormatter(evolve.LOGGING_SETTINGS["formatter"])
@@ -178,7 +186,7 @@ def main():
             x_grid = np.geomspace(args.x_grid_ini, 1.0, args.x_grid_points)
         if args.actions == "produce_eko":
             tcard, opcard = eko_utils.construct_eko_cards(
-                args.theoryID,
+                nnpdf_theory,
                 args.q_fin,
                 args.q_points,
                 x_grid,
@@ -188,7 +196,7 @@ def main():
             )
         elif args.actions == "produce_eko_photon":
             tcard, opcard = eko_utils.construct_eko_photon_cards(
-                args.theoryID, args.q_fin, x_grid, args.q_gamma, op_card_info, theory_card_info
+                nnpdf_theory, args.q_fin, x_grid, args.q_gamma, op_card_info, theory_card_info
             )
         solve(tcard, opcard, args.dump)
 

n3fit/src/n3fit/tests/conftest.py

@@ -1,9 +1,22 @@
 """
     Add markers for pytest
 """
+
 import sys
+
 import pytest
 
+from validphys.loader import FallbackLoader
+
+THEORYID = 399
+
+
+@pytest.fixture(scope='module')
+def nnpdf_theory_card():
+    """Return a theory card already loaded as a dictionary"""
+    th = FallbackLoader().check_theoryID(THEORYID)
+    return th.get_description()
+
 
 def pytest_runtest_setup(item):
     ALL = {"darwin", "linux"}

n3fit/src/n3fit/tests/test_evolven3fit.py

@@ -118,16 +118,15 @@ def test_utils():
     assert ID == 162
 
 
-def test_eko_utils(tmp_path):
+def test_eko_utils(tmp_path, nnpdf_theory_card):
     # Testing construct eko cards
-    theoryID = 162
     q_fin = 100
     q_points = 5
     x_grid = [1.0e-3, 0.1, 1.0]
     pto = 2
     comments = "Test"
     t_card, op_card = eko_utils.construct_eko_cards(
-        theoryID,
+        nnpdf_theory_card,
         q_fin,
         q_points,
         x_grid,
@@ -140,7 +139,7 @@ def test_eko_utils(tmp_path):
         t_card_dict["order"][0] == pto + 1
     )  # This is due to a different convention in eko orders due to QED
     np.testing.assert_allclose(op_card_dict["xgrid"], x_grid)
-    # In theory 162 the charm threshold is at 1.51
+    # In theory 399 the charm threshold is at 1.51
     # and we should find two entries, one for nf=3 and another one for nf=4
     np.testing.assert_allclose(op_card_dict["mugrid"][0], (1.51, 3))
     np.testing.assert_allclose(op_card_dict["mugrid"][1], (1.51, 4))