Provide an unique truth for all theory cards

n3fit/src/evolven3fit/eko_utils.py

@@ -58,6 +58,10 @@ def construct_eko_cards(
     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
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory

nnpdf_data/nnpdf_data/theory.py

@@ -0,0 +1,99 @@
+"""
+This module provides an unique source and definition for all the possible parameters
+that a theory card can contain.
+"""
+
+import dataclasses
+import logging
+
+DEPRECATED_KEYS = ["MaxNfAs", "SxRes", "SxOrd" "EScaleVar", "Qedref", "global_nx"]
+
+log = logging.getLogger(__name__)
+
+
+class TheoryCardError(Exception):
+    pass
+
+
+@dataclasses.dataclass(frozen=True)
+class TheoryCard:
+    ID: int  # ID number of the theory
+    PTO: int  # Perturbative order (0 = LO, 1 = NLO, 2 = NNLO ...)
+    FNS: str  # Flavor number scheme (e.g. FONLL-C)
+    DAMP: int  # Whether a damping function is applied or not for FONLL
+    IC: int  # 0 = perturbative charm only , 1 = intrinsic charm allowed
+    ModEv: str  # DGLAP evolution solution method (e.g. EXA or TRN)
+    XIR: float  # Renormalization scale over the hard scattering scale ratio
+    XIF: float  # Factorization scale over the hard scattering scale ratio
+    NfFF: int  # Number of active flavors, only for FFNS or FFN0 schemes
+    QED: int  # Max order of alpha_qed in the splitting functions
+    mc: float  # [GeV] charm mass
+    Qmc: float  # [GeV] MSbar mass reference scale of the charm
+    kcThr: float  # Threshold ratio of the charm
+    mb: float  # # [GeV] bottom mass
+    Qmb: float  # [GeV] MSbar mass reference scale of the bottom
+    kbThr: float  # Threshold ratio of the bottom
+    mt: float  # # [GeV] top mass
+    Qmt: float  # [GeV] MSbar mass reference scale of the top
+    ktThr: float  # Threshold ratio of the top
+    CKM: list[float]  # CKM matrix elements (running on the columns first, i.e. CKM_11 is CKM[0] and CKM_12 is CKM[1] and so on)
+    MZ: float  # [GeV] Mass of Z
+    MW: float  # [GeV] Mass of W
+    GF: float  # Fermi constant
+    SIN2TW: float
+    TMC: int  # Include target mass corrections: 0 = disabled, 1 = leading twist, 2 = higher twist approximated, 3 = higher twist exact
+    MP: float  # [GeV] Mass of the proton
+    Comments: str  # Comments on the theory
+    MaxNfPdf: int = 5  # Used by pineko and the photon module to define the thresholds
+    # Fit theory parameters default
+    nf0: int = 4  # Number of active flavors at the parametrization scale Q0
+    Q0: float = 1.65  # [GeV] Parametrization scale
+    nfref: int = 5  # Number of active flavors at Qref
+    Qref: float = 91.2  # [GeV] Reference scale for alphas and alphaqed
+    alphas: float = 0.118  # Value of alpha_s at the scale Qref
+    alphaqed: float = 0.007496252  # Values of alpha QED at the scale Qref
+    # Evolution parameters
+    HQ: str = "POLE" # Heavy quark mass scheme,  POLE for pole masses (default), MSBAR for running masses (used currently only in eko).
+    IterEv: int = None  # iterations for the evolution of the PDF. Defaults to 40 when ModEv = EXA
+    ModSV: str = None  # Scale variations method in EKO (expanded or exponentiated)
+    DAMPPOWERc: int = None  # Power of the damping factor in FONLL for the c
+    DAMPPOWERb: int = None  # Power of the damping factor in FONLL for the b
+    # N3LO parameters
+    n3lo_ad_variation: list = dataclasses.field(
+        default_factory=lambda: 7 * [0]
+    )  # N3LO anomalous dimension variations
+    n3lo_cf_variation: int = (
+        0 # N3LO coefficient functions variation: -1 = lower bound, 0 = central, 1 = upper bound
+    )
+    use_fhmruvv: bool = (
+        False # N3LO splitting functions approximation: if True use the FHMRUVV parametrization, otherwise use EKO parametrization.
+    )
+    ###### Keys for compatibility with old NNPDF theories, their values will be dropped immediately after reading to avoid problems
+    ###### they will be set to ``None`` immediately after loading the theory
+    MaxNfAs: int = None
+    SxRes: int = None
+    SxOrd: str = None
+    EScaleVar: int = None
+    Qedref: float = None
+    global_nx: int = None
+
+    def __post_init__(self):
+        """Drop deprecated keys and apply some checks"""
+        if self.Qedref is not None and self.QED != 0:
+            # Check that nobody is trying to use this with a wrong Qedref!
+            if self.Qedref != self.Qref:
+                self._raise_or_warn(
+                    f"Trying to use {self.ID} with {self.Qedref} != {self.Qref}. This is not supported!"
+                )
+
+        for key in DEPRECATED_KEYS:
+            object.__setattr__(self, key, None)
+
+    def _raise_or_warn(self, msg):
+        """Raise an error for new theories and a warning for old ones"""
+        if self.ID >= 600:
+            raise TheoryCardError(msg)
+        log.warning(msg)
+
+    def asdict(self):
+        return dataclasses.asdict(self)

nnpdf_data/nnpdf_data/theorydbutils.py

@@ -5,58 +5,15 @@
 low level utilities for querying the theory database file and representing the
 data as a python object.
 """
-from dataclasses import asdict, dataclass
 from functools import lru_cache
 from pathlib import Path
 
 import pandas as pd
 
+from .theory import TheoryCard
 from .utils import parse_yaml_inp
 
 
-@dataclass(frozen=True)
-class _TheoryCard:
-    ID: int
-    PTO: int
-    FNS: str
-    DAMP: int
-    IC: int
-    ModEv: str
-    XIR: float
-    XIF: float
-    NfFF: int
-    MaxNfAs: int
-    MaxNfPdf: int
-    Q0: float
-    alphas: float
-    Qref: float
-    QED: int
-    alphaqed: float
-    Qedref: float
-    SxRes: int
-    SxOrd: str
-    HQ: str
-    mc: float
-    Qmc: float
-    kcThr: float
-    mb: float
-    Qmb: float
-    kbThr: float
-    mt: float
-    Qmt: float
-    ktThr: float
-    CKM: list[float]
-    MZ: float
-    MW: float
-    GF: float
-    SIN2TW: float
-    TMC: int
-    MP: float
-    Comments: str
-    global_nx: int
-    EScaleVar: int
-
-
 class TheoryNotFoundInDatabase(Exception):
     pass
 
@@ -67,7 +24,8 @@ def parse_theory_card(theory_card):
     Returns the theory as a dictionary
     """
     if theory_card.exists():
-        return asdict(parse_yaml_inp(theory_card, _TheoryCard))
+        tcard = parse_yaml_inp(theory_card, TheoryCard)
+        return tcard.asdict()
     raise TheoryNotFoundInDatabase(f"Theory card {theory_card} not found")
 
 

nnpdf_data/nnpdf_data/utils.py

@@ -3,6 +3,11 @@
 import ruamel.yaml as yaml
 from validobj import ValidationError, parse_input
 
+try:
+    Loader = yaml.CLoader
+except AttributeError:
+    Loader = yaml.Loader
+
 
 def parse_yaml_inp(input_yaml, spec):
     """
@@ -12,7 +17,7 @@ def parse_yaml_inp(input_yaml, spec):
     https://validobj.readthedocs.io/en/latest/examples.html#yaml-line-numbers
     """
     input_yaml = pathlib.Path(input_yaml)
-    inp = yaml.round_trip_load(input_yaml.open("r", encoding="utf-8"))
+    inp = yaml.load(input_yaml.read_text(encoding="utf-8"), Loader=Loader)
     try:
         return parse_input(inp, spec)
     except ValidationError as e:

validphys2/src/validphys/photon/alpha.py

@@ -1,4 +1,5 @@
 """Module that implements the alpha running"""
+
 import numpy as np
 from scipy.integrate import solve_ivp
 
@@ -10,6 +11,7 @@
 
 class Alpha:
     """Class that solves the RGE for alpha_qed"""
+
     def __init__(self, theory, q2max):
         self.theory = theory
         self.alpha_em_ref = theory["alphaqed"]
@@ -20,11 +22,11 @@ def __init__(self, theory, q2max):
         self.thresh_c = self.theory["kcThr"] * self.theory["mc"]
         self.thresh_b = self.theory["kbThr"] * self.theory["mb"]
         self.thresh_t = self.theory["ktThr"] * self.theory["mt"]
-        if self.theory["MaxNfAs"] <= 5:
+        if self.theory["MaxNfPdf"] <= 5:
             self.thresh_t = np.inf
-        if self.theory["MaxNfAs"] <= 4:
+        if self.theory["MaxNfPdf"] <= 4:
             self.thresh_b = np.inf
-        if self.theory["MaxNfAs"] <= 3:
+        if self.theory["MaxNfPdf"] <= 3:
             self.thresh_c = np.inf
 
         if self.theory["ModEv"] == "TRN":
@@ -145,7 +147,7 @@ def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl):
         # at LO in QED the TRN solution is exact
         if self.qed_order == 1:
             return self.alphaem_fixed_flavor_trn(q, alphaem_ref, qref, nf, nl)
-        
+
         u = 2 * np.log(q / qref)
 
         # solve RGE
@@ -214,7 +216,7 @@ def compute_betas(self, regions):
         for nf, nl in regions:
             betas_qed[(nf, nl)] = [
                 beta.beta_qed_aem2(nf, nl) / (4 * np.pi),
-                beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2 if self.theory['QED'] == 2 else 0.,
+                beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2 if self.theory['QED'] == 2 else 0.0,
             ]
         return betas_qed
 
@@ -249,6 +251,7 @@ def find_region(self, q):
 def _rge(_t, alpha_t, beta_qed_vec):
     """RGEs for the running of alphaem"""
     rge_qed = -(alpha_t**2) * (
-        beta_qed_vec[0] + np.sum([alpha_t ** (k + 1) * beta for k, beta in enumerate(beta_qed_vec[1:])])
+        beta_qed_vec[0]
+        + np.sum([alpha_t ** (k + 1) * beta for k, beta in enumerate(beta_qed_vec[1:])])
     )
     return rge_qed

validphys2/src/validphys/photon/compute.py

@@ -1,4 +1,5 @@
 """Module that calls fiatlux to add the photon PDF."""
+
 import logging
 import tempfile
 
@@ -55,11 +56,9 @@ def __init__(self, theoryid, lux_params, replicas):
 
         theory = theoryid.get_description()
         fiatlux_runcard = FIATLUX_DEFAULT
-        fiatlux_runcard["qed_running"] = bool(np.isclose(theory["Qedref"], theory["Qref"]))
-        # cast explicitly from np.bool_ to bool otherwise problems in dumping it
-        # TODO: for the time being, we trigger alphaem running if Qedref=Qref.
-        # This is going to be changed in favor of a bool em_running
-        # in the runcard
+        # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
+        # This may be changed in the future in favor of a bool em_running in the runcard
+        fiatlux_runcard["qed_running"] = True
         fiatlux_runcard["mproton"] = float(theory["MP"])
 
         # precision on final integration of double integral

validphys2/src/validphys/tests/photon/test_alpha.py

@@ -34,7 +34,7 @@ def test_set_thresholds_alpha_em():
 
         # test all the thresholds
         np.testing.assert_almost_equal(alpha.alpha_em_ref, theory["alphaqed"])
-        np.testing.assert_almost_equal(alpha.thresh[(5, 3)], theory["Qedref"])
+        np.testing.assert_almost_equal(alpha.thresh[(5, 3)], theory["Qref"])
         np.testing.assert_almost_equal(alpha.thresh[(4, 3)], theory["mb"])
         np.testing.assert_almost_equal(alpha.thresh[(4, 2)], constants.MTAU)
         np.testing.assert_almost_equal(alpha.thresh[(3, 2)], theory["mc"])
@@ -46,7 +46,7 @@ def test_set_thresholds_alpha_em():
         np.testing.assert_almost_equal(alpha.alphaem_thresh[(5, 3)], theory["alphaqed"])
         np.testing.assert_almost_equal(
             alpha.alphaem_thresh[(4, 3)],
-            alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5, 3),
+            alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qref"], 5, 3),
         )
         np.testing.assert_almost_equal(
             alpha.alphaem_thresh[(4, 2)],
@@ -102,7 +102,7 @@ def test_couplings_exa():
                 alphaem=theory["alphaqed"],
                 scale=theory["Qref"],
                 num_flavs_ref=None,
-                max_num_flavs=theory["MaxNfAs"],
+                max_num_flavs=theory["MaxNfPdf"],
                 em_running=True,
             )
         )
@@ -138,7 +138,7 @@ def test_couplings_exa():
             np.testing.assert_allclose(
                 alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6
             )
-        for nf in range(3, theory["MaxNfAs"]):
+        for nf in range(3, theory["MaxNfPdf"]):
             np.testing.assert_allclose(
                 alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)],
                 eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi,

validphys2/src/validphys/tests/photon/test_compute.py

@@ -10,8 +10,8 @@
 from validphys.core import PDF as PDFset
 from validphys.loader import FallbackLoader
 from validphys.photon import structure_functions as sf
-from validphys.photon.compute import FIATLUX_DEFAULT, Photon
 from validphys.photon.alpha import Alpha
+from validphys.photon.compute import FIATLUX_DEFAULT, Photon
 
 from ..conftest import PDF, THEORY_QED
 
@@ -70,7 +70,7 @@ def test_photon():
         # runcard
         fiatlux_default = FIATLUX_DEFAULT.copy()
         fiatlux_default['mproton'] = theory['MP']
-        fiatlux_default["qed_running"] = bool(np.isclose(theory["Qedref"], theory["Qref"]))
+        fiatlux_default["qed_running"] = False
         fiatlux_default["q2_max"] = float(f2.q2_max)
         fiatlux_default["eps_base"] = fiatlux_runcard["eps_base"]