Merge pull request #99 from NNPDF/num_fonll

Numerical FONLL

benchmarks/bench_evolve.py

@@ -6,13 +6,47 @@
 import pineappl
 import pytest
 import yaml
-from eko import couplings as sc
 
 import pineko
 import pineko.evolve
 import pineko.theory_card
 
 
+@pytest.mark.parametrize("theoryid,is_mixed", [(400, False), (208, True)])
+def benchmark_write_operator_card_from_file_num_fonll(
+    tmp_path, test_files, test_configs, theoryid, is_mixed
+):
+    tcard = pineko.theory_card.load(theoryid)
+    tcards_path_list = pineko.fonll.dump_tcards(tcard, tmp_path, theoryid)
+    pine_path = (
+        test_files
+        / "data"
+        / "grids"
+        / "400"
+        / "HERA_NC_225GEV_EP_SIGMARED.pineappl.lz4"
+    )
+    default_path = test_files / "data" / "operator_cards" / "400" / "_template.yaml"
+    num_opcard = 7 if is_mixed else 5
+    targets_path_list = [
+        tmp_path / f"test_opcard_{num}.yaml" for num in range(num_opcard)
+    ]
+    for target_path, tcard_path in zip(targets_path_list, tcards_path_list):
+        with open(tcard_path, encoding="utf-8") as f:
+            tcard = yaml.safe_load(f)
+        _x_grid, _q2_grid = pineko.evolve.write_operator_card_from_file(
+            pine_path, default_path, target_path, tcard
+        )
+    # Check if the opcards are ok
+    for opcard_path, cfg in zip(
+        targets_path_list,
+        pineko.fonll.MIXED_FNS_CONFIG if is_mixed else pineko.fonll.FNS_CONFIG,
+    ):
+        with open(opcard_path, encoding="utf-8") as f:
+            ocard = yaml.safe_load(f)
+        for entry in ocard["mugrid"]:
+            assert entry[1] == int(cfg.nf)
+
+
 def benchmark_write_operator_card_from_file(tmp_path, test_files, test_configs):
     pine_path = test_files / "data/grids/400/HERA_NC_225GEV_EP_SIGMARED.pineappl.lz4"
     default_path = test_files / "data/operator_cards/400/_template.yaml"

benchmarks/bench_fonll.py

@@ -0,0 +1,28 @@
+import numpy as np
+import pytest
+import yaml
+
+import pineko
+
+
+@pytest.mark.parametrize("theoryid,is_mixed", [(400, False), (208, True)])
+def benchmark_produce_fonll_tcards(
+    tmp_path, test_files, test_configs, theoryid, is_mixed
+):
+    tcard = pineko.theory_card.load(theoryid)
+    tcard_paths_list = pineko.fonll.dump_tcards(tcard, tmp_path, theoryid)
+    # Check they are correct
+    theorycards = []
+    for path in tcard_paths_list:
+        with open(path, encoding="UTF-8") as f:
+            theorycards.append(yaml.safe_load(f))
+    base_pto = pineko.fonll.FNS_BASE_PTO[tcard["FNS"]]
+    for num_fonll_tcard, cfg in zip(
+        theorycards,
+        pineko.fonll.MIXED_FNS_CONFIG if is_mixed else pineko.fonll.FNS_CONFIG,
+    ):
+        po = int(base_pto) + (cfg.delta_pto if is_mixed else 0)
+        assert num_fonll_tcard["FNS"] == cfg.scheme
+        assert num_fonll_tcard["NfFF"] == int(cfg.nf)
+        assert num_fonll_tcard["PTO"] == po - 1 if is_mixed and cfg.asy else po
+        assert num_fonll_tcard["FONLLParts"] == cfg.parts

src/pineko/check.py

@@ -116,29 +116,49 @@ def check_grid_and_eko_compatible(pineappl_grid, operators, xif, max_as, max_al)
         raise ValueError("x grid in pineappl grid and eko operator are NOT compatible!")
 
 
-def is_fonll_b(fns, lumi):
-    """Check if the fktable we are computing is a DIS FONLL-B fktable.
+def is_dis(lumi):
+    """Check if the fktable we are computing is a DIS fktable.
 
     Parameters
     ----------
-    fns : str
-          flavor number scheme (from the theory card)
     lumi : list(list(tuple))
            luminosity info
 
     Returns
     -------
     bool
-        true if the fktable is a FONLL-B DIS fktable
+        true if the fktable is a DIS fktable
     """
-    for lists in lumi:
-        for el in lists:
+    for entry in lumi:
+        for el in entry:
             if (not islepton(el[0])) and (not islepton(el[1])):
-                # in this case we are sure it is not DIS so for sure it is not FONLL-B
+                # If neither of the incoming particles is a lepton we are sure
+                # it is not DIS
                 return False
-    if fns == "FONLL-B":
-        return True
-    return False
+    return True
+
+
+def is_fonll_mixed(fns, lumi):
+    """Check if the fktable we are computing is FONLL-B, FONLL-D or, in general, a mixed FONLL fktable.
+
+    Parameters
+    ----------
+    fns : str
+          flavor number scheme (from the theory card)
+    lumi : list(list(tuple))
+           luminosity info
+
+    Returns
+    -------
+    bool
+        true if the fktable is a mixed FONLL DIS fktable
+    """
+    return is_dis(lumi) and fns in ["FONLL-B", "FONLL-D"]
+
+
+def is_num_fonll(fns):
+    """Check if the FNS is a nFONLL FNS."""
+    return fns in ["FONLL-FFNS", "FONLL-FFN0"]
 
 
 def orders(grid, max_as, max_al) -> list:

src/pineko/cli/__init__.py

@@ -1,3 +1,13 @@
 """CLI entry point."""
-from . import check, compare, convolute, gen_sv, kfactor, opcard, scaffold, theory_
+from . import (
+    check,
+    compare,
+    convolute,
+    fonll,
+    gen_sv,
+    kfactor,
+    opcard,
+    scaffold,
+    theory_,
+)
 from ._base import command

src/pineko/cli/fonll.py

@@ -0,0 +1,174 @@
+"""CLI entry point to FONLL."""
+import pathlib
+
+import click
+import rich
+
+from .. import configs, fonll, parser, theory_card
+from ._base import command
+
+config_setting = click.option(
+    "-c",
+    "--configs",
+    "cfg",
+    default=None,
+    type=click.Path(resolve_path=True, path_type=pathlib.Path),
+    help="Explicitly specify config file (it has to be a valid TOML file).",
+)
+
+
+class TheoryCardError(Exception):
+    """Raised when asked for FONLL theory cards with an original tcard as input that is not asking for FONLL."""
+
+
+class InconsistentInputsError(Exception):
+    """Raised if the inputs are not consistent with FONLL."""
+
+
+def cfgpath(name, grid):
+    """Path of the fktable in 'name' called 'grid' if it exists, else None."""
+    path = configs.configs["paths"]["fktables"] / name / grid
+    return path if path.exists() else None
+
+
+def grids_names(yaml_file):
+    """Return the list of the grids in the yaml file."""
+    yaml_content = parser._load_yaml(yaml_file)
+    # Turn the operands and the members into paths (and check all of them exist)
+    ret = []
+    for operand in yaml_content["operands"]:
+        for member in operand:
+            ret.append(f"{member}.{parser.EXT}")
+    return ret
+
+
+@command.command("combine_fonll")
+@click.argument("theoryID", type=int)
+@click.argument("dataset", type=str)
+@click.option("--FFNS3", type=int, help="theoryID containing the ffns3 fktable")
+@click.option("--FFN03", type=int, help="theoryID containing the ffn03 fktable")
+@click.option("--FFNS4", type=int, help="theoryID containing the ffns4 fktable")
+@click.option("--FFNS4til", type=int, help="theoryID containing the ffns4til fktable")
+@click.option("--FFNS4bar", type=int, help="theoryID containing the ffns4bar fktable")
+@click.option("--FFN04", type=int, help="theoryID containing the ffn04 fktable")
+@click.option("--FFNS5", type=int, help="theoryID containing the ffns5 fktable")
+@click.option("--FFNS5til", type=int, help="theoryID containing the ffns5til fktable")
+@click.option("--FFNS5bar", type=int, help="theoryID containing the ffns5bar fktable")
+@click.option("--overwrite", is_flag=True, help="Allow files to be overwritten")
+@config_setting
+def subcommand(
+    theoryid,
+    dataset,
+    ffns3,
+    ffn03,
+    ffns4,
+    ffns4til,
+    ffns4bar,
+    ffn04,
+    ffns5,
+    ffns5til,
+    ffns5bar,
+    overwrite,
+    cfg,
+):
+    """Combine the different FKs needed to produce the FONLL prescription."""
+    path = configs.detect(cfg)
+    base_configs = configs.load(path)
+    configs.configs = configs.defaults(base_configs)
+    if cfg is not None:
+        print(f"Configurations loaded from '{path}'")
+
+    # Checks
+
+    if not ffns3 or not ffn03:
+        raise InconsistentInputsError("ffns3 and/or ffn03 is not provided.")
+
+    if any([ffns4, ffns4til, ffns4bar]):
+        if ffns4:
+            if any([ffns4til, ffns4bar]):
+                raise InconsistentInputsError(
+                    "If ffns4 is provided no ffnstil or ffnsbar should be provided."
+                )
+        else:
+            if ffns4til is None or ffns4bar is None:
+                raise InconsistentInputsError(
+                    "if ffnstil is provided also ffnsbar should be provided, and vice versa."
+                )
+    else:
+        raise InconsistentInputsError("ffns4 is not provided.")
+
+    # Do we consider two masses, i.e. mc and mb
+    two_masses = False
+    if any([ffns5, ffns5til, ffns5bar]):
+        two_masses = True
+        if ffns5:
+            if any([ffns5til, ffns5bar]):
+                raise InconsistentInputsError(
+                    "If ffns5 is provided no ffnstil or ffnsbar should be provided."
+                )
+        else:
+            if ffns5til is None or ffns5bar is None:
+                raise InconsistentInputsError(
+                    "if ffnstil is provided also ffnsbar should be provided, and vice versa."
+                )
+
+    if (ffn04 is None and two_masses) or (ffn04 is not None and not two_masses):
+        raise InconsistentInputsError(
+            "If two masses are to be considered, both ffn04 and the nf=5 coefficient should be provided"
+        )
+
+    # Get theory info
+    tcard = theory_card.load(theoryid)
+    if not "DAMPPOWER" in tcard:
+        if tcard["DAMP"] != 0:
+            raise InconsistentInputsError("If DAMP is set, set also DAMPPOWER")
+        tcard["DAMPPOWER"] = None
+    # Getting the paths to the grids
+    grids_name = grids_names(configs.configs["paths"]["ymldb"] / f"{dataset}.yaml")
+    for grid in grids_name:
+        # Checking if it already exists
+        new_fk_path = configs.configs["paths"]["fktables"] / str(theoryid) / grid
+        if new_fk_path.exists():
+            if not overwrite:
+                rich.print(
+                    f"[green]Success:[/] skipping existing FK Table {new_fk_path}"
+                )
+                return
+        fonll.produce_combined_fk(
+            *(
+                cfgpath(str(name), grid)
+                for name in (
+                    ffns3,
+                    ffn03,
+                    ffns4,
+                    ffns4til,
+                    ffns4bar,
+                    ffn04,
+                    ffns5,
+                    ffns5til,
+                    ffns5bar,
+                )
+            ),
+            theoryid,
+            damp=(tcard["DAMP"], tcard["DAMPPOWER"]),
+            cfg=cfg,
+        )
+        if new_fk_path.exists():
+            rich.print(f"[green]Success:[/] Wrote FK table to {new_fk_path}")
+        else:
+            rich.print(f"[red]Failure:[/]")
+
+
+@command.command("fonll_tcards")
+@click.argument("theoryID", type=int)
+@config_setting
+def fonll_tcards(theoryid, cfg):
+    """Produce the FONLL tcards starting from the original tcard given by the theoryID."""
+    path = configs.detect(cfg)
+    base_configs = configs.load(path)
+    configs.configs = configs.defaults(base_configs)
+    tcard = theory_card.load(theoryid)
+    tcard_parent_path = theory_card.path(theoryid).parent
+    if "FONLL" not in tcard["FNS"]:
+        raise TheoryCardError("The theorycard does not correspond to an FONLL scheme.")
+    fonll.dump_tcards(tcard, tcard_parent_path, theoryid)

src/pineko/evolve.py

@@ -105,7 +105,8 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):
 
     """
     # Add a +1 to the orders for the difference in convention between nnpdf and pineappl
-    is_fns = int(check.is_fonll_b(tcard["FNS"], pineappl_grid.lumi()))
+    # NB: This would not happen for nFONLL
+    is_fns = int(check.is_fonll_mixed(tcard["FNS"], pineappl_grid.lumi()))
     max_as = 1 + tcard["PTO"] + is_fns
     max_al = 1 + tcard["QED"]
     # ... in order to create a mask ...
@@ -128,9 +129,14 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):
         matching_scales=heavy_quarks.MatchingScales(masses * thresholds_ratios),
         origin=(tcard["Q0"] ** 2, tcard["nf0"]),
     )
-    operators_card["mugrid"] = [
-        (float(np.sqrt(q2)), int(nf_default(q2, atlas))) for q2 in q2_grid
-    ]
+    # If we are producing nFONLL FKs we need to look to NfFF...
+    if check.is_num_fonll(tcard["FNS"]):
+        nf = tcard["NfFF"]
+        operators_card["mugrid"] = [(float(np.sqrt(q2)), int(nf)) for q2 in q2_grid]
+    else:
+        operators_card["mugrid"] = [
+            (float(np.sqrt(q2)), nf_default(q2, atlas)) for q2 in q2_grid
+        ]
     if "integrability_version" in pineappl_grid.key_values():
         x_grid = evol_info.x1
         x_grid = np.append(x_grid, 1.0)