diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 988a824120..03c74376d5 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -39,7 +39,7 @@ requirements: - prompt_toolkit - validobj - pineappl >=0.8.2 - - eko >=0.14.2,<0.15 + - eko >=0.15.1,<0.16 - fiatlux - sphinx >=5.0.2 - joblib diff --git a/doc/sphinx/.gitignore b/doc/sphinx/.gitignore index 7f1d36e99b..86f516e98f 100644 --- a/doc/sphinx/.gitignore +++ b/doc/sphinx/.gitignore @@ -1,3 +1,4 @@ # Ignore auto generated module references source/modules source/theories.csv +source/theories_central.csv diff --git a/doc/sphinx/Makefile b/doc/sphinx/Makefile index 63bfbfb512..352a64ea2e 100644 --- a/doc/sphinx/Makefile +++ b/doc/sphinx/Makefile @@ -19,6 +19,7 @@ help: @if test $@ != "clean"; then \ sphinx-apidoc -f -o ./$(SOURCEDIR)/modules/validphys ../../validphys2/src/validphys/ ; \ sphinx-apidoc -f -o ./$(SOURCEDIR)/modules/n3fit-code ../../n3fit/src/n3fit/ ; \ + sphinx-apidoc -f -o ./$(SOURCEDIR)/modules/n3fit-code ../../n3fit/src/evolven3fit/ ; \ fi @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/n3fit/src/evolven3fit/cli.py b/n3fit/src/evolven3fit/cli.py index 84f07e4ce9..3943ba6319 100644 --- a/n3fit/src/evolven3fit/cli.py +++ b/n3fit/src/evolven3fit/cli.py @@ -2,7 +2,7 @@ def cli_evolven3fit( - configuration_folder, q_fin, q_points, op_card_info, theory_card_info, force, load, dump, ncores + configuration_folder, q_fin, q_points, op_card_info, theory_card_info, force, load, dump ): """Evolves the fitted PDFs. @@ -23,13 +23,5 @@ def cli_evolven3fit( """ utils.check_is_a_fit(configuration_folder) return evolve.evolve_fit( - configuration_folder, - q_fin, - q_points, - op_card_info, - theory_card_info, - force, - load, - dump, - ncores, + configuration_folder, q_fin, q_points, op_card_info, theory_card_info, force, load, dump ) diff --git a/n3fit/src/evolven3fit/eko_utils.py b/n3fit/src/evolven3fit/eko_utils.py index b3cea6373f..29602138c7 100644 --- a/n3fit/src/evolven3fit/eko_utils.py +++ b/n3fit/src/evolven3fit/eko_utils.py @@ -1,7 +1,6 @@ import logging from typing import Any, Dict, Optional -from ekobox.cards import _operator as default_op_card import numpy as np from eko.io import runcards @@ -30,6 +29,11 @@ NFREF_DEFAULT = 5 +def _eko_theory_from_nnpdf_theory(nnpdf_theory): + """Takes an NNPDF theory and returns an eko theory""" + return runcards.Legacy(nnpdf_theory, {}).new_theory + + def construct_eko_cards( nnpdf_theory, q_fin, @@ -54,10 +58,7 @@ def construct_eko_cards( # 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 a NNPDF theory - legacy_class = runcards.Legacy(theory, {}) - theory_card = legacy_class.new_theory + theory_card = _eko_theory_from_nnpdf_theory(theory) # construct mugrid @@ -127,10 +128,7 @@ def construct_eko_photon_cards( # Now make sure the Legacy class still gets a Qedref, which is equal to Qref theory["Qedref"] = theory["Qref"] theory["MaxNfAs"] = theory["MaxNfPdf"] - - # 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 + theory_card = _eko_theory_from_nnpdf_theory(theory) # The photon needs to be evolved down to Q0 q_fin = theory["Q0"] @@ -177,10 +175,31 @@ def build_opcard(op_card_dict, theory, x_grid, mu0, mugrid): if op_card_dict is None: op_card_dict = {} - op_card = default_op_card - - op_card.update({"mu0": mu0, "mugrid": mugrid}) + # Taken from cards.py https://github.com/NNPDF/eko/blob/master/src/ekobox/cards.py + # 7735fdb + op_card = dict( + init=(1.65, 4), + mugrid=[(100.0, 5)], + xgrid=np.geomspace(1e-7, 1.0, 50).tolist(), + configs=dict( + # These three values might be set by op_card_dict + ev_op_iterations=10, + n_integration_cores=1, + polarized=False, + # ... but these we fix here + ev_op_max_order=[10, 0], + interpolation_polynomial_degree=4, + interpolation_is_log=True, + scvar_method=None, + inversion_method=None, + evolution_method="iterate-exact", + time_like=False, + ), + debug=dict(skip_singlet=False, skip_non_singlet=False), + ) + op_card["init"] = (mu0, theory["nf0"]) + op_card["mugrid"] = mugrid op_card["xgrid"] = x_grid # Specify the evolution options and defaults differently from TRN / EXA diff --git a/n3fit/src/evolven3fit/evolve.py b/n3fit/src/evolven3fit/evolve.py index fb10a4db84..d50b7444ac 100644 --- a/n3fit/src/evolven3fit/evolve.py +++ b/n3fit/src/evolven3fit/evolve.py @@ -1,16 +1,17 @@ from collections import defaultdict +import dataclasses import json import logging import pathlib import sys from ekobox import apply, genpdf, info_file -from joblib import Parallel, delayed import numpy as np -import psutil import eko from eko import basis_rotation, runner +from eko.interpolation import XGrid +from eko.io import manipulate from validphys.utils import yaml_safe from . import eko_utils, utils @@ -24,19 +25,9 @@ "level": logging.DEBUG, } -NUM_CORES = psutil.cpu_count(logical=False) - def evolve_fit( - fit_folder, - q_fin, - q_points, - op_card_dict, - theory_card_dict, - force, - eko_path, - dump_eko=None, - ncores=1, + fit_folder, q_fin, q_points, op_card_dict, theory_card_dict, force, eko_path, dump_eko=None ): """ Evolves all the fitted replica in fit_folder/nnfit @@ -106,10 +97,7 @@ def evolve_fit( else: raise ValueError(f"dump_eko not provided and {eko_path=} not found") - with eko.EKO.edit(eko_path) as eko_op: - x_grid_obj = eko.interpolation.XGrid(x_grid) - eko.io.manipulate.xgrid_reshape(eko_op, targetgrid=x_grid_obj, inputgrid=x_grid_obj) - + # Open the EKO in read-only mode, if it needs to be manipulated keep it in memory with eko.EKO.read(eko_path) as eko_op: # Read the cards directly from the eko to make sure they are consistent theory = eko_op.theory_card @@ -118,6 +106,42 @@ def evolve_fit( _logger.debug(f"Theory card: {json.dumps(theory.raw)}") _logger.debug(f"Operator card: {json.dumps(op.raw)}") + eko_original_xgrid = eko_op.xgrid + if XGrid(x_grid) != eko_original_xgrid: + # If the xgrid of the eko is not directly usable, construct a copy in memory + # by replacing the internal inventory of operators in a readonly copy + new_xgrid = XGrid(x_grid) + new_metadata = dataclasses.replace(eko_op.metadata, xgrid=new_xgrid) + + new_operators = {} + for target_key in eko_op.operators: + elem = eko_op[target_key.ep] + + if eko_op.metadata.version == "0.13.4": + # For eko 0.13.4 xgrid is the internal interpolation so we need to check + # whether the rotation is truly needed + # this means checking whether the operator shape matches the grid + oplen = elem.operator.shape[-1] + if oplen != len(eko_original_xgrid): + # The operator and its xgrid have different shape + # either prepare an identity, or this EKO is not supported + if oplen != len(x_grid): + raise ValueError( + f"The operator at {eko_path} is not usable, version not supported" + ) + eko_original_xgrid = XGrid(x_grid) + + new_operators[target_key] = manipulate.xgrid_reshape( + elem, + eko_original_xgrid, + op.configs.interpolation_polynomial_degree, + targetgrid=XGrid(x_grid), + inputgrid=XGrid(x_grid), + ) + + new_inventory = dataclasses.replace(eko_op.operators, cache=new_operators) + eko_op = dataclasses.replace(eko_op, metadata=new_metadata, operators=new_inventory) + # Modify the info file with the fit-specific info info = info_file.build(theory, op, 1, info_update={}) info["NumMembers"] = "REPLACE_NREP" @@ -126,21 +150,49 @@ def evolve_fit( info["XMax"] = float(x_grid[-1]) # Save the PIDs in the info file in the same order as in the evolution info["Flavors"] = basis_rotation.flavor_basis_pids - info["NumFlavors"] = theory.heavy.num_flavs_max_pdf + info.setdefault("NumFlavors", 5) dump_info_file(usr_path, info) - def _wrap_evolve(pdf, replica): - evolved_blocks = evolve_exportgrid(pdf, eko_op, x_grid) - dump_evolved_replica(evolved_blocks, usr_path, int(replica.removeprefix("replica_"))) - - # Choose the number of cores to be the Minimal value - nb_cores = min(NUM_CORES, abs(ncores)) - Parallel(n_jobs=nb_cores)( - delayed(_wrap_evolve)(pdf, r) for r, pdf in initial_PDFs_dict.items() - ) + # Read the information from all the sorted replicas into what eko wants + n_replicas = len(initial_PDFs_dict) + all_replicas = [] + for rep_idx in range(1, n_replicas + 1): + # swap photon position to match eko.basis_rotation.flavor_basis_pids + pdfgrid = np.array(initial_PDFs_dict[f"replica_{rep_idx}"]["pdfgrid"]) + pdfgrid = np.append(pdfgrid[:, -1].reshape(x_grid.size, 1), pdfgrid[:, :-1], axis=1) + # and divide by x + all_replicas.append(pdfgrid.T / x_grid) + + # output is {(Q2, nf): (replica, flavour, x)} + all_evolved, _ = apply.apply_grids(eko_op, np.array(all_replicas)) + + # Now, replica by replica, break into nf blocks + targetgrid = eko_op.xgrid.tolist() + by_nf = defaultdict(list) + for q2, nf in sorted(eko_op.evolgrid, key=lambda ep: ep[1]): + by_nf[nf].append(q2) + q2block_per_nf = {nf: sorted(q2s) for nf, q2s in by_nf.items()} + + for replica in range(n_replicas): + blocks = [] + for nf, q2grid in q2block_per_nf.items(): + + def pdf_xq2(pid, x, Q2): + x_idx = targetgrid.index(x) + pid_idx = info["Flavors"].index(pid) + return x * all_evolved[(Q2, nf)][replica][pid_idx][x_idx] + + block = genpdf.generate_block( + pdf_xq2, + xgrid=targetgrid, + sorted_q2grid=q2grid, + pids=basis_rotation.flavor_basis_pids, + ) + blocks.append(block) + dump_evolved_replica(blocks, usr_path, replica + 1) # remove folder: - # The function dump_evolved_replica dumps the replica files in a temporary folder + # The function dump_evolved_replica uses a temporary folder # We need then to remove it after fixing the position of those replica files (usr_path / "nnfit" / usr_path.stem).rmdir() @@ -169,52 +221,6 @@ def load_fit(usr_path): return pdf_dict -def evolve_exportgrid(exportgrid, eko, x_grid): - """ - Evolves the provided exportgrid for the desired replica with the eko and returns the evolved block - - Parameters - ---------- - exportgrid: dict - exportgrid of pdf at fitting scale - eko: eko object - eko operator for evolution - xgrid: list - xgrid to be used as the targetgrid - Returns - ------- - : list(np.array) - list of evolved blocks - """ - # construct LhapdfLike object - pdf_grid = np.array(exportgrid["pdfgrid"]).transpose() - pdf_to_evolve = utils.LhapdfLike(pdf_grid, exportgrid["q20"], x_grid) - # evolve pdf - evolved_pdf = apply.apply_pdf(eko, pdf_to_evolve) - # generate block to dump - targetgrid = eko.bases.targetgrid.tolist() - - # Finally separate by nf block (and order per nf/q) - by_nf = defaultdict(list) - for q, nf in sorted(eko.evolgrid, key=lambda ep: ep[1]): - by_nf[nf].append(q) - q2block_per_nf = {nf: sorted(qs) for nf, qs in by_nf.items()} - - blocks = [] - for nf, q2grid in q2block_per_nf.items(): - - def pdf_xq2(pid, x, Q2): - x_idx = targetgrid.index(x) - return x * evolved_pdf[(Q2, nf)]["pdfs"][pid][x_idx] - - block = genpdf.generate_block( - pdf_xq2, xgrid=targetgrid, sorted_q2grid=q2grid, pids=basis_rotation.flavor_basis_pids - ) - blocks.append(block) - - return blocks - - def dump_evolved_replica(evolved_blocks, usr_path, replica_num): """ Dump the evolved replica given by evolved_block as the replica num "replica_num" in diff --git a/n3fit/src/evolven3fit/q2grids.py b/n3fit/src/evolven3fit/q2grids.py index a10dced774..ef4385b273 100644 --- a/n3fit/src/evolven3fit/q2grids.py +++ b/n3fit/src/evolven3fit/q2grids.py @@ -1,12 +1,13 @@ """ - Definition of default Q2 grids +Definition of default Q2 grids - This file includes: +This file includes: - ``Q2GRID_DEFAULT``: default NNPDF Q2 grid for evolution (55 points, starts at Q=1GeV) - ``Q2GRID_NNPDF40``: q2 grid used in the fits for the NNPDF4.0 release (49 points, starts at Q=1.65 GeV) - ``Q2GRID_Nf03``: q2 grid used in the perturvative charm fits for the NNPDF4.0 release (48 points, starts at Q=1GeV) + - ``Q2GRID_DEFAULT``: default NNPDF Q2 grid for evolution (55 points, starts at Q=1GeV) + - ``Q2GRID_NNPDF40``: q2 grid used in the fits for the NNPDF4.0 release (49 points, starts at Q=1.65 GeV) + - ``Q2GRID_Nf03``: q2 grid used in the perturvative charm fits for the NNPDF4.0 release (48 points, starts at Q=1GeV) """ + import numpy as np Q2GRID_NNPDF40 = ( diff --git a/n3fit/src/evolven3fit/utils.py b/n3fit/src/evolven3fit/utils.py index cdd649d972..7d31e3c4ef 100644 --- a/n3fit/src/evolven3fit/utils.py +++ b/n3fit/src/evolven3fit/utils.py @@ -10,51 +10,6 @@ from .q2grids import Q2GRID_DEFAULT, Q2GRID_NNPDF40 -class LhapdfLike: - """ - Class which emulates lhapdf but only for an initial condition PDF (i.e. with only one q2 value). - - Q20 is the fitting scale fo the pdf and it is the only available scale for the objects of this class. - - X_GRID is the grid of x values on top of which the pdf is interpolated. - - PDF_GRID is a dictionary containing the pdf grids at fitting scale for each pid. - """ - - def __init__(self, pdf_grid, q20, x_grid): - self.pdf_grid = pdf_grid - self.q20 = q20 - self.x_grid = x_grid - self.funcs = [ - interp1d(self.x_grid, self.pdf_grid[pid], kind="cubic") for pid in range(len(PIDS_DICT)) - ] - - def xfxQ2(self, pid, x, q2): - """Return the value of the PDF for the requested pid, x value and, whatever the requested - q2 value, for the fitting q2. - - Parameters - ---------- - - pid: int - pid index of particle - x: float - x-value - q2: float - Q square value - - Returns - ------- - : float - x * PDF value - """ - return self.funcs[list(PIDS_DICT.values()).index(PIDS_DICT[pid])](x) - - def hasFlavor(self, pid): - """Check if the requested pid is in the PDF.""" - return pid in PIDS_DICT - - def read_runcard(usr_path): """Read the runcard and return the relevant information for evolven3fit""" return yaml_safe.load((usr_path / "filter.yml").read_text(encoding="UTF-8")) diff --git a/n3fit/src/n3fit/scripts/evolven3fit.py b/n3fit/src/n3fit/scripts/evolven3fit.py index 277de17858..03ddd42f24 100644 --- a/n3fit/src/n3fit/scripts/evolven3fit.py +++ b/n3fit/src/n3fit/scripts/evolven3fit.py @@ -19,6 +19,27 @@ _logger = logging.getLogger(__name__) +def _add_production_options(parser): + """Adds the production options that are shared between the normal and the QED EKOs to their parsers""" + parser.add_argument("theoryID", type=int, help="ID of the theory used to produce the eko") + parser.add_argument("dump", type=pathlib.Path, help="Path where the EKO is dumped") + parser.add_argument( + "-i", "--x-grid-ini", default=None, type=float, help="Starting point of the x-grid" + ) + parser.add_argument( + "-p", "--x-grid-points", default=None, type=int, help="Number of points of the x-grid" + ) + parser.add_argument("-n", "--n-cores", type=int, default=1, help="Number of cores to be used") + parser.add_argument('--use_polarized', action='store_true', help="Use polarized evolution") + parser.add_argument( + "-e", + "--ev-op-iterations", + type=int, + default=None, + help="ev_op_iterations for the EXA theory. Overrides the settings given in the theory card.", + ) + + def construct_eko_parser(subparsers): parser = subparsers.add_parser( "produce_eko", @@ -30,14 +51,7 @@ def construct_eko_parser(subparsers): provided number of points. The eko will be dumped in the provided path.""" ), ) - parser.add_argument("theoryID", type=int, help="ID of the theory used to produce the eko") - parser.add_argument("dump", type=pathlib.Path, help="Path where the EKO is dumped") - parser.add_argument( - "-i", "--x-grid-ini", default=None, type=float, help="Starting point of the x-grid" - ) - parser.add_argument( - "-p", "--x-grid-points", default=None, type=int, help="Number of points of the x-grid" - ) + _add_production_options(parser) parser.add_argument( "--legacy40", action="store_true", @@ -56,14 +70,7 @@ def construct_eko_photon_parser(subparsers): provided number of points. The eko will be dumped in the provided path.""" ), ) - parser.add_argument("theoryID", type=int, help="ID of the theory used to produce the eko") - parser.add_argument("dump", type=pathlib.Path, help="Path where the EKO is dumped") - parser.add_argument( - "-i", "--x-grid-ini", default=None, type=float, help="Starting point of the x-grid" - ) - parser.add_argument( - "-p", "--x-grid-points", default=None, type=int, help="Number of points of the x-grid" - ) + _add_production_options(parser) parser.add_argument( "-g", "--q-gamma", default=100, type=float, help="Scale at which the photon is generated" ) @@ -76,7 +83,7 @@ def construct_evolven3fit_parser(subparsers): help="Evolves the fitted PDFs. The q_grid starts at the Q0 given by the theory but the last point is q_fin and its number of points can be specified by q_points. If a path is given for the dump option, the eko will be dumped in that path after the computation. If a path is given for the load option, the eko to be used for the evolution will be loaded from that path. The two options are mutually exclusive.", ) parser.add_argument( - "configuration_folder", help="Path to the folder containing the (pre-DGLAP) fit result" + "fit_folder", help="Path to the folder containing the (pre-DGLAP) fit result" ) parser.add_argument( "-l", "--load", type=pathlib.Path, default=None, help="Path of the EKO to be loaded" @@ -99,37 +106,15 @@ def evolven3fit_new(): def main(): - parser = ArgumentParser( - description="evolven3fit - a script with tools to evolve PDF fits", - usage="""evolven3fit [-h] [-q Q_FIN] [-p Q_POINTS] [-n N_CORES] [-e EV_OP_ITERATIONS] [--use-fhmruvv] - {produce_eko,produce_eko_photon,evolve} [fit folder] - - Note that with the now removed apfel-based version of `evolven3fit` the syntax was - `evolven3fit [fit folder] [number of replicas]`. This syntax is no longer supported in the - eko-based version of evolven3fit. - """, - ) - parser.add_argument('--use_polarized', action='store_true', help="Use polarized evolution") + parser = ArgumentParser(description="evolven3fit - a script with tools to evolve PDF fits") parser.add_argument( "-q", "--q-fin", type=float, default=None, help="Final q-value of the evolution" ) parser.add_argument( "-p", "--q-points", type=int, default=None, help="Number of q points for the evolution" ) - parser.add_argument("-n", "--n-cores", type=int, default=1, help="Number of cores to be used") parser.add_argument("--no-net", action="store_true", help="Emulates validphys' --no-net") - parser.add_argument( - "-e", - "--ev-op-iterations", - type=int, - default=None, - help="ev_op_iterations for the EXA theory. Overrides the settings given in the theory card.", - ) - parser.add_argument( - "--use-fhmruvv", - action="store_true", - help="Use the FHMRUVV N3LO splitting splitting functions", - ) + subparsers = parser.add_subparsers(title="actions", dest="actions") construct_eko_parser(subparsers) construct_eko_photon_parser(subparsers) @@ -137,15 +122,9 @@ def main(): args = parser.parse_args() - op_card_info = { - "configs": {"n_integration_cores": args.n_cores, "polarized": args.use_polarized} - } - if args.ev_op_iterations is not None: - op_card_info["configs"]["ev_op_iterations"] = args.ev_op_iterations - + op_card_info = {} + # Here we do not allow any modification of the theory card, for the moment. theory_card_info = {} - if args.use_fhmruvv: - theory_card_info["use_fhmruvv"] = args.use_fhmruvv if args.no_net: loader = Loader() @@ -154,7 +133,7 @@ def main(): if args.actions == "evolve": - fit_folder = pathlib.Path(args.configuration_folder) + fit_folder = pathlib.Path(args.fit_folder) if args.load is None: theoryID = utils.get_theoryID_from_runcard(fit_folder) _logger.info(f"Loading eko from theory {theoryID}") @@ -171,7 +150,6 @@ def main(): args.force, eko_path, None, - args.n_cores, ) else: # If we are in the business of producing an eko, do some checks before starting: @@ -200,6 +178,14 @@ def main(): ) else: x_grid = np.geomspace(args.x_grid_ini, 1.0, args.x_grid_points) + + # Prepare the op card config + op_card_info = { + "configs": {"n_integration_cores": args.n_cores, "polarized": args.use_polarized} + } + if args.ev_op_iterations is not None: + op_card_info["configs"]["ev_op_iterations"] = args.ev_op_iterations + if args.actions == "produce_eko": tcard, opcard = eko_utils.construct_eko_cards( nnpdf_theory, diff --git a/n3fit/src/n3fit/tests/test_evolven3fit.py b/n3fit/src/n3fit/tests/test_evolven3fit.py index 34039cb8ba..106682a882 100644 --- a/n3fit/src/n3fit/tests/test_evolven3fit.py +++ b/n3fit/src/n3fit/tests/test_evolven3fit.py @@ -98,17 +98,6 @@ def test_generate_q2grid(): def test_utils(): - # Testing the fake LHAPDF class - q20 = 1.65**2 - x_grid = np.geomspace(1.0e-7, 1.0, 30) - fake_grids = [[x * (1.0 - x) for x in x_grid] for _ in PIDS_DICT.keys()] - pdf_grid = {pid: v for pid, v in zip(range(len(PIDS_DICT)), fake_grids)} - my_PDF = utils.LhapdfLike(pdf_grid, q20, x_grid) - assert my_PDF.hasFlavor(6) - assert not my_PDF.hasFlavor(0) - for pid in PIDS_DICT: - for x in x_grid: - np.testing.assert_allclose(my_PDF.xfxQ2(pid, x, q20), x * (1.0 - x)) # Testing read_runcard runcard = utils.read_runcard(REGRESSION_FOLDER) assert isinstance(runcard["description"], str) diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py index 18f45a6d3b..49c6fa7b01 100644 --- a/n3fit/src/n3fit/vpinterface.py +++ b/n3fit/src/n3fit/vpinterface.py @@ -26,7 +26,6 @@ import numpy as np import numpy.linalg as la -from n3fit.backends import PREPROCESSING_LAYER_ALL_REPLICAS from validphys.arclength import arc_lengths, integrability_number from validphys.core import PDF, MCStats from validphys.covmats import covmat_from_systematics, sqrt_covmat @@ -239,9 +238,13 @@ def get_preprocessing_factors(self, replica=None): otherwise outputs a Nx2 array where [:,0] are alphas and [:,1] betas """ # If no replica is explicitly requested, get the preprocessing layer for the first model + # remember replicas start counting at one if replica is None: replica = 1 - # Replicas start counting in 1 so: + + # Keep the import here to avoid loading the backend when it is not necessary + from n3fit.backends import PREPROCESSING_LAYER_ALL_REPLICAS + preprocessing_layer = self._models[replica - 1].get_layer(PREPROCESSING_LAYER_ALL_REPLICAS) alphas_and_betas = None diff --git a/pyproject.toml b/pyproject.toml index 8b41b6b48d..cc05b18be8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -73,7 +73,7 @@ reportengine = ">=0.32" psutil = "*" tensorflow = "*" keras = "^3.1" -eko = "^0.14.1" +eko = "^0.15.1" joblib = "*" # Hyperopt hyperopt = "*" diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 6466c4ddb7..4b5cf29622 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -8,6 +8,7 @@ from scipy.interpolate import interp1d import yaml +from eko import basis_rotation from eko.io import EKO from n3fit.io.writer import XGRID from validphys.n3fit_data import replica_luxseed @@ -165,26 +166,27 @@ def compute_photon_array(self, replica): # TODO : the different x points could be even computed in parallel # Load eko and reshape it - with EKO.read(self.path_to_eko_photon) as eko: + with EKO.read(self.path_to_eko_photon) as eko_photon: # TODO : if the eko has not the correct grid we have to reshape it # it has to be done inside vp-setupfit - # construct PDFs - pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID))) - for j, pid in enumerate(eko.bases.inputpids): - if pid == 22: - pdfs_init[j] = photon_qin - ph_id = j - else: - if pid not in self.luxpdfset.flavors: + # NB: the eko should contain a single operator + for _, elem in eko_photon.items(): + eko_op = elem.operator + + pdfs_init = np.zeros_like(eko_op[0, 0]) + for j, pid in enumerate(basis_rotation.flavor_basis_pids): + if pid == 22: + pdfs_init[j] = photon_qin + ph_id = j + elif pid not in self.luxpdfset.flavors: continue - pdfs_init[j] = np.array( - [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID] - ) + else: + pdfs_init[j] = np.array( + [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID] + ) - # Apply EKO to PDFs - for _, elem in eko.items(): - pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init) + pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init) photon_Q0 = pdfs_final[ph_id] diff --git a/validphys2/src/validphys/tests/photon/test_alpha.py b/validphys2/src/validphys/tests/photon/test_alpha.py index 118ee2499b..07654b6748 100644 --- a/validphys2/src/validphys/tests/photon/test_alpha.py +++ b/validphys2/src/validphys/tests/photon/test_alpha.py @@ -100,9 +100,7 @@ def test_couplings_exa(): dict( alphas=theory["alphas"], alphaem=theory["alphaqed"], - scale=theory["Qref"], - num_flavs_ref=None, - max_num_flavs=theory["MaxNfPdf"], + ref=(theory["Qref"], theory["MaxNfPdf"]), em_running=True, ) ) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index ff0b64710b..274f6673fa 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -4,6 +4,7 @@ import numpy as np import yaml +from eko import basis_rotation from eko.io import EKO from n3fit.io.writer import XGRID from validphys.api import API @@ -94,8 +95,8 @@ def test_photon(): photon_fiatlux_qin = np.array([lux.EvaluatePhoton(x, eko.mu20).total for x in XGRID]) photon_fiatlux_qin /= XGRID # construct PDFs - pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID))) - for j, pid in enumerate(eko.bases.inputpids): + pdfs_init = np.zeros((len(basis_rotation.flavor_basis_pids), len(XGRID))) + for j, pid in enumerate(basis_rotation.flavor_basis_pids): if pid == 22: pdfs_init[j] = photon_fiatlux_qin ph_id = j