diff --git a/.gitignore b/.gitignore index 14960be3b1..0fbe0bcd81 100644 --- a/.gitignore +++ b/.gitignore @@ -378,6 +378,9 @@ ENV/ env.bak/ venv.bak/ +conda-build +conda_bld/ + # Spyder project settings .spyderproject .spyproject diff --git a/conda-recipe/conda_build_config.yaml b/conda-recipe/conda_build_config.yaml index 71cab38400..a15ac5fdef 100644 --- a/conda-recipe/conda_build_config.yaml +++ b/conda-recipe/conda_build_config.yaml @@ -1,7 +1,7 @@ #For some reason the resolver decides to use an old version of numpy #without this numpy: - - 1.22 + - 1.21 pin_run_as_build: lhapdf: x.x.x diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index db81a1b4ef..510adb0861 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -50,7 +50,11 @@ requirements: - sphinx_rtd_theme >0.5 - sphinxcontrib-bibtex - curio >=1.0 - - pineappl >=0.5.2 + - pineappl >=0.5.8 + - eko ==0.10.2 + - banana-hep ==0.6.6 + - lz4 >=3.1.10 # see https://github.com/conda-forge/eko-feedstock/issues/10 + test: requires: diff --git a/n3fit/setup.py b/n3fit/setup.py index e68a67f4fe..e81c54809b 100644 --- a/n3fit/setup.py +++ b/n3fit/setup.py @@ -13,6 +13,7 @@ entry_points = {'console_scripts': ['n3fit = n3fit.scripts.n3fit_exec:main', + 'evolven3fit_new = n3fit.scripts.evolven3fit_new:main', 'vp-setupfit = n3fit.scripts.vp_setupfit:main', 'varflavors = n3fit.scripts.varflavors:main', ] diff --git a/n3fit/src/evolven3fit_new/__init__.py b/n3fit/src/evolven3fit_new/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/n3fit/src/evolven3fit_new/cli.py b/n3fit/src/evolven3fit_new/cli.py new file mode 100644 index 0000000000..92ecbdf231 --- /dev/null +++ b/n3fit/src/evolven3fit_new/cli.py @@ -0,0 +1,244 @@ +from . import evolve, utils + +import numpy as np + +XGRID = np.array( + [ + 1.00000000000000e-09, + 1.29708482343957e-09, + 1.68242903474257e-09, + 2.18225315420583e-09, + 2.83056741739819e-09, + 3.67148597892941e-09, + 4.76222862935315e-09, + 6.17701427376180e-09, + 8.01211109898438e-09, + 1.03923870607245e-08, + 1.34798064073805e-08, + 1.74844503691778e-08, + 2.26788118881103e-08, + 2.94163370300835e-08, + 3.81554746595878e-08, + 4.94908707232129e-08, + 6.41938295708371e-08, + 8.32647951986859e-08, + 1.08001422993829e-07, + 1.40086873081130e-07, + 1.81704331793772e-07, + 2.35685551545377e-07, + 3.05703512595323e-07, + 3.96522309841747e-07, + 5.14321257236570e-07, + 6.67115245136676e-07, + 8.65299922973143e-07, + 1.12235875241487e-06, + 1.45577995547683e-06, + 1.88824560514613e-06, + 2.44917352454946e-06, + 3.17671650028717e-06, + 4.12035415232797e-06, + 5.34425265752090e-06, + 6.93161897806315e-06, + 8.99034258238145e-06, + 1.16603030112258e-05, + 1.51228312288769e-05, + 1.96129529349212e-05, + 2.54352207134502e-05, + 3.29841683435992e-05, + 4.27707053972016e-05, + 5.54561248105849e-05, + 7.18958313632514e-05, + 9.31954227979614e-05, + 1.20782367731330e-04, + 1.56497209466554e-04, + 2.02708936328495e-04, + 2.62459799331951e-04, + 3.39645244168985e-04, + 4.39234443000422e-04, + 5.67535660104533e-04, + 7.32507615725537e-04, + 9.44112105452451e-04, + 1.21469317686978e-03, + 1.55935306118224e-03, + 1.99627451141338e-03, + 2.54691493736552e-03, + 3.23597510213126e-03, + 4.09103436509565e-03, + 5.14175977083962e-03, + 6.41865096062317e-03, + 7.95137940306351e-03, + 9.76689999624100e-03, + 1.18876139251364e-02, + 1.43298947643919e-02, + 1.71032279460271e-02, + 2.02100733925079e-02, + 2.36463971369542e-02, + 2.74026915728357e-02, + 3.14652506132444e-02, + 3.58174829282429e-02, + 4.04411060163317e-02, + 4.53171343973807e-02, + 5.04266347950069e-02, + 5.57512610084339e-02, + 6.12736019390519e-02, + 6.69773829498255e-02, + 7.28475589986517e-02, + 7.88703322292727e-02, + 8.50331197801452e-02, + 9.13244910278679e-02, + 9.77340879783772e-02, + 1.04252538208639e-01, + 1.10871366547237e-01, + 1.17582909372878e-01, + 1.24380233801599e-01, + 1.31257062945031e-01, + 1.38207707707289e-01, + 1.45227005135651e-01, + 1.52310263065985e-01, + 1.59453210652156e-01, + 1.66651954293987e-01, + 1.73902938455578e-01, + 1.81202910873333e-01, + 1.88548891679097e-01, + 1.95938145999193e-01, + 2.03368159629765e-01, + 2.10836617429103e-01, + 2.18341384106561e-01, + 2.25880487124065e-01, + 2.33452101459503e-01, + 2.41054536011681e-01, + 2.48686221452762e-01, + 2.56345699358723e-01, + 2.64031612468684e-01, + 2.71742695942783e-01, + 2.79477769504149e-01, + 2.87235730364833e-01, + 2.95015546847664e-01, + 3.02816252626866e-01, + 3.10636941519503e-01, + 3.18476762768082e-01, + 3.26334916761672e-01, + 3.34210651149156e-01, + 3.42103257303627e-01, + 3.50012067101685e-01, + 3.57936449985571e-01, + 3.65875810279643e-01, + 3.73829584735962e-01, + 3.81797240286494e-01, + 3.89778271981947e-01, + 3.97772201099286e-01, + 4.05778573402340e-01, + 4.13796957540671e-01, + 4.21826943574548e-01, + 4.29868141614175e-01, + 4.37920180563205e-01, + 4.45982706956990e-01, + 4.54055383887562e-01, + 4.62137890007651e-01, + 4.70229918607142e-01, + 4.78331176755675e-01, + 4.86441384506059e-01, + 4.94560274153348e-01, + 5.02687589545177e-01, + 5.10823085439086e-01, + 5.18966526903235e-01, + 5.27117688756998e-01, + 5.35276355048428e-01, + 5.43442318565661e-01, + 5.51615380379768e-01, + 5.59795349416641e-01, + 5.67982042055800e-01, + 5.76175281754088e-01, + 5.84374898692498e-01, + 5.92580729444440e-01, + 6.00792616663950e-01, + 6.09010408792398e-01, + 6.17233959782450e-01, + 6.25463128838069e-01, + 6.33697780169485e-01, + 6.41937782762089e-01, + 6.50183010158361e-01, + 6.58433340251944e-01, + 6.66688655093089e-01, + 6.74948840704708e-01, + 6.83213786908386e-01, + 6.91483387159697e-01, + 6.99757538392251e-01, + 7.08036140869916e-01, + 7.16319098046733e-01, + 7.24606316434025e-01, + 7.32897705474271e-01, + 7.41193177421404e-01, + 7.49492647227008e-01, + 7.57796032432224e-01, + 7.66103253064927e-01, + 7.74414231541921e-01, + 7.82728892575836e-01, + 7.91047163086478e-01, + 7.99368972116378e-01, + 8.07694250750291e-01, + 8.16022932038457e-01, + 8.24354950923382e-01, + 8.32690244169987e-01, + 8.41028750298844e-01, + 8.49370409522600e-01, + 8.57715163684985e-01, + 8.66062956202683e-01, + 8.74413732009721e-01, + 8.82767437504206e-01, + 8.91124020497459e-01, + 8.99483430165226e-01, + 9.07845617001021e-01, + 9.16210532771399e-01, + 9.24578130473112e-01, + 9.32948364292029e-01, + 9.41321189563734e-01, + 9.49696562735755e-01, + 9.58074441331298e-01, + 9.66454783914439e-01, + 9.74837550056705e-01, + 9.83222700304978e-01, + 9.91610196150662e-01, + 1.00000000000000e00, + ] +) + + +def cli_evolven3fit_new( + configuration_folder, + q_fin, + q_points, + op_card_info, + theory_card_info, + dump, + load, + force, +): + """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 just one of the + two is not specified by the user, the default grid + will be used. + + 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. + """ + 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, + ) diff --git a/n3fit/src/evolven3fit_new/eko_utils.py b/n3fit/src/evolven3fit_new/eko_utils.py new file mode 100644 index 0000000000..076e38144d --- /dev/null +++ b/n3fit/src/evolven3fit_new/eko_utils.py @@ -0,0 +1,83 @@ +from ekobox import gen_theory, gen_op +from eko import run_dglap + +from validphys.loader import Loader +from . import utils + +from typing import Any, Dict, Optional +import logging + +_logger = logging.getLogger(__name__) + + +def construct_eko_cards( + theoryID, + q_fin, + q_points, + x_grid, + op_card_dict: Optional[Dict[str, Any]] = None, + theory_card_dict: Optional[Dict[str, Any]] = None, +): + """ + 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. + 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. + """ + if theory_card_dict is None: + theory_card_dict = {} + if op_card_dict is None: + op_card_dict = {} + # theory_card construction + theory = Loader().check_theoryID(theoryID).get_description() + theory.pop("FNS") + theory.update(theory_card_dict) + theory_card = gen_theory.gen_theory_card(theory["PTO"], theory["Q0"], update=theory) + # construct operator card + q2_grid = utils.generate_q2grid( + theory["Q0"], + q_fin, + q_points, + {theory["mb"]: theory["kbThr"], theory["mt"]: theory["ktThr"]}, + ) + op_card = gen_op.gen_op_card(q2_grid, update={"interpolation_xgrid": x_grid}) + op_card.update(op_card_dict) + return theory_card, op_card + + +def construct_eko_for_fit(theory_card, op_card, save_path=None): + """ + Construct the eko operator needed for evolution of fitted pdfs + + Parameters + ---------- + theory_card: dict + theory card to use for the eko + op_card: dict + operator card to use for the eko + save_path: pathlib.Path + path where the eko will be saved (the eko + won't be saved if save_path is None) + Returns + ------- + : eko.output.Output + eko operator + """ + # generate eko operator + if save_path is not None: + if not save_path.parent.exists(): + raise FileNotFoundError( + f"Path where eko should be dumped does not exist: {save_path}" + ) + eko_op = run_dglap(theory_card, op_card) + if save_path is not None: + # Here we want to catch all possible exceptions in order to avoid losing the computed eko + try: + _logger.info(f"Saving computed eko to : {save_path}") + eko_op.dump_tar(save_path) + except: + _logger.error(f"Error saving the eko to : {save_path}") + pass + return eko_op diff --git a/n3fit/src/evolven3fit_new/evolve.py b/n3fit/src/evolven3fit_new/evolve.py new file mode 100644 index 0000000000..2cf0180040 --- /dev/null +++ b/n3fit/src/evolven3fit_new/evolve.py @@ -0,0 +1,220 @@ +import pathlib +import logging +import sys + +import numpy as np +from reportengine.compat import yaml + +from ekobox import genpdf, gen_info +from ekomark import apply +from eko import basis_rotation as br +from eko import output +from validphys.loader import Loader + +from . import utils, eko_utils + + +_logger = logging.getLogger(__name__) + +LOG_FILE = "evolven3fit_new.log" + +LOGGING_SETTINGS = {"formatter" : "%(asctime)s %(name)s/%(levelname)s: %(message)s", + "level" : logging.INFO +} + +def evolve_fit( + fit_folder, + q_fin, + q_points, + op_card_dict, + theory_card_dict, + force, + eko_path=None, + dump_eko=None, +): + """ + Evolves all the fitted replica in fit_folder/nnfit + + Parameters + ---------- + + fit_folder: str or pathlib.Path + path to the folder containing the fit + q_fin: float + final point of the q_grid + q_points: int + number of points in the q_grid + op_card_dict: dict + user settings for the op_card + theory_card_dict: dict + user settings for the t_card + force: bool + whether to force the evolution to be done again + eko_path: str or pathlib.Path + path where the eko is stored (if None the eko will be + recomputed) + dump_eko: str or pathlib.Path + path where the eko is dumped (if None the eko won't be + stored) + """ + log_file = pathlib.Path(fit_folder) / LOG_FILE + if log_file.exists(): + if force: + log_file.unlink() + else: + raise FileExistsError( + f"Log file already exists: {log_file} evolven3fit_new has already been run?" + ) + log_file = logging.FileHandler(log_file) + stdout_log = logging.StreamHandler(sys.stdout) + for log in [log_file, stdout_log]: + log.setLevel(LOGGING_SETTINGS["level"]) + log.setFormatter( + logging.Formatter(LOGGING_SETTINGS["formatter"]) + ) + for logger in (_logger, *[logging.getLogger("eko")]): + logger.handlers = [] + logger.setLevel(LOGGING_SETTINGS["level"]) + logger.addHandler(log_file) + logger.addHandler(stdout_log) + + usr_path = pathlib.Path(fit_folder) + initial_PDFs_dict = load_fit(usr_path) + x_grid = np.array( + initial_PDFs_dict[list(initial_PDFs_dict.keys())[0]]["xgrid"] + ).astype(np.float) + theoryID = utils.get_theoryID_from_runcard(usr_path) + theory, op = eko_utils.construct_eko_cards( + theoryID, q_fin, q_points, x_grid, op_card_dict, theory_card_dict + ) + if eko_path is not None: + eko_path = pathlib.Path(eko_path) + _logger.info(f"Loading eko from : {eko_path}") + eko_op = output.Output.load_tar(eko_path) + else: + try: + _logger.info(f"Loading eko from theory {theoryID}") + theory_eko_path = (Loader().check_theoryID(theoryID).path)/'eko.tar' + eko_op = output.Output.load_tar(theory_eko_path) + except FileNotFoundError: + _logger.info(f"eko not found in theory {theoryID}, we will construct it") + eko_op = eko_utils.construct_eko_for_fit(theory, op, dump_eko) + eko_op.xgrid_reshape(targetgrid=x_grid, inputgrid=x_grid) + info = gen_info.create_info_file(theory, op, 1, info_update={}) + info["NumMembers"] = "REPLACE_NREP" + info["ErrorType"] = "replicas" + info["AlphaS_Qs"] = list(eko_op["Q2grid"]) + info["XMin"] = float(x_grid[0]) + info["XMax"] = float(x_grid[-1]) + dump_info_file(usr_path, info) + for replica in initial_PDFs_dict.keys(): + evolved_block = evolve_exportgrid(initial_PDFs_dict[replica], eko_op, x_grid) + dump_evolved_replica( + evolved_block, usr_path, int(replica.removeprefix("replica_")) + ) + # remove folder: + # The function dump_evolved_replica dumps the replica files in a temporary folder + # We need then to remove it after fixing the position of those replica files + (usr_path / "nnfit" / usr_path.stem).rmdir() + + +def load_fit(usr_path): + """ + Loads all the replica pdfs at fitting scale in usr_path and return the exportgrids + + Parameters + ---------- + + usr_path: pathlib.Path + path to the folder containing the fit + + Returns + ------- + + : dict + exportgrids info + """ + nnfitpath = usr_path / "nnfit" + pdf_dict = {} + for yaml_file in nnfitpath.glob("replica_*/*.exportgrid"): + data = yaml.safe_load(yaml_file.read_text(encoding="UTF-8")) + pdf_dict[yaml_file.parent.stem] = data + 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 + ------- + : np.array + evolved block + """ + # 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 = list(eko["targetgrid"]) + + def ev_pdf(pid, x, Q2): + return x * evolved_pdf[Q2]["pdfs"][pid][targetgrid.index(x)] + + block = genpdf.generate_block( + ev_pdf, + xgrid=targetgrid, + Q2grid=list(eko["Q2grid"]), + pids=br.flavor_basis_pids, + ) + return block + + +def dump_evolved_replica(evolved_block, usr_path, replica_num): + """ + Dump the evolved replica given by evolved_block as the replica num "replica_num" in + the folder usr_path/nnfit/replica_/usr_path.stem.dat + + Parameters + ---------- + evolved_block: numpy.array + block of an evolved PDF + usr_path: pathlib.Path + path of the fit folder + replica_num: int + replica number + """ + path_where_dump = usr_path / "nnfit" / usr_path.stem + # create folder to dump the evolved replica if it does not exist + path_where_dump.mkdir(exist_ok=True) + to_write_in_head = f"PdfType: replica\nFromMCReplica: {replica_num}\n" + genpdf.export.dump_blocks( + path_where_dump, replica_num, [evolved_block], pdf_type=to_write_in_head + ) + # fixing_replica_path + utils.fix_replica_path(usr_path, replica_num) + + +def dump_info_file(usr_path, info): + """ + Dump the info file given by info in the folder usr_path/nnfit/usr_path.stem.info. + + Parameters + ---------- + usr_path: pathlib.Path + path of the fit folder + info: dict + info of the fit + """ + path_where_dump = usr_path / "nnfit" / usr_path.stem + genpdf.export.dump_info(path_where_dump, info) + utils.fix_info_path(usr_path) diff --git a/n3fit/src/evolven3fit_new/utils.py b/n3fit/src/evolven3fit_new/utils.py new file mode 100644 index 0000000000..afbfa6f445 --- /dev/null +++ b/n3fit/src/evolven3fit_new/utils.py @@ -0,0 +1,201 @@ +import shutil +import pathlib +from scipy.interpolate import interp1d +import numpy as np +from reportengine.compat import yaml +from validphys.pdfbases import PIDS_DICT + +DEFAULT_Q2GRID = np.array( + [ + 1.6500000e00, + 1.7874388e00, + 1.9429053e00, + 2.1193749e00, + 2.3204100e00, + 2.5502944e00, + 2.8142025e00, + 3.1184122e00, + 3.4705775e00, + 3.8800751e00, + 4.3584516e00, + 4.9200000e00, + 4.9200000e00, + 5.5493622e00, + 6.2897452e00, + 7.1650687e00, + 8.2052867e00, + 9.4481248e00, + 1.0941378e01, + 1.2745972e01, + 1.4940062e01, + 1.7624572e01, + 2.0930715e01, + 2.5030298e01, + 3.0149928e01, + 3.6590777e01, + 4.4756282e01, + 5.5191298e01, + 6.8637940e01, + 8.6115921e01, + 1.0903923e02, + 1.3938725e02, + 1.7995815e02, + 2.3474820e02, + 3.0952544e02, + 4.1270732e02, + 5.5671861e02, + 7.6011795e02, + 1.0509694e03, + 1.4722574e03, + 2.0906996e03, + 3.0112909e03, + 4.4016501e03, + 6.5333918e03, + 9.8535186e03, + 1.5109614e04, + 2.3573066e04, + 3.7444017e04, + 6.0599320e04, + 1.0000000e05, + ] + ) **2 + +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_new""" + return yaml.safe_load((usr_path / "filter.yml").read_text(encoding="UTF-8")) + + +def get_theoryID_from_runcard(usr_path): + """Return the theoryID from the runcard""" + # read the runcard + my_runcard = read_runcard(usr_path) + return my_runcard["theory"]["theoryid"] + + +def generate_q2grid(Q0, Qfin, Q_points, match_dict): + """Generate the q2grid used in the final evolved pdfs or use the default grid if Qfin or Q_points is + not provided. + + match_dict contains the couples (mass : factor) where factor is the number to be multiplied to mass + in order to obtain the relative matching scale. + """ + if Qfin is None and Q_points is None: + return DEFAULT_Q2GRID + elif Qfin is None or Q_points is None: + raise ValueError( + "q_fin and q_points must be specified either both or none of them" + ) + else: + grids = [] + Q_ini = Q0 + num_points_list = [] + for masses in match_dict: + match_scale = masses * match_dict[masses] + # Fraction of the total points to be included in this batch is proportional + # to the log of the ratio between the initial scale and final scale of the + # batch itself (normalized to the same log of the global initial and final + # scales) + if match_scale < Qfin: + frac_of_point = np.log(match_scale / Q_ini) / np.log(Qfin / Q0) + num_points = int( + Q_points * frac_of_point + ) + num_points_list.append(num_points) + grids.append(np.geomspace(Q_ini**2, match_scale**2, num=num_points)) + Q_ini = match_scale + num_points = Q_points - sum(num_points_list) + grids.append(np.geomspace(Q_ini**2, Qfin**2, num=num_points)) + return np.concatenate(grids).tolist() + + +def fix_info_path(usr_path): + """Fix the location of the info file from the folder nnfit/usr_path to + just nnfit + + Examples + -------- + Starting from the info path + initial_info_file_path = "/myfolder/myfit/nnfit/myfit/myfit.info" + and using this function with usr_path = "/myfolder/myfit", one gets + final_info_file_path = "/myfolder/myfit/nnfit/myfit.info" + """ + nnfit = usr_path / "nnfit" + info_file = usr_path.stem + ".info" + info_file_path = nnfit / usr_path.stem / info_file + dest_path_info = nnfit / info_file + shutil.move(info_file_path, dest_path_info) + + +def fix_replica_path(usr_path, replica_num): + """Fix the location of the dat file of the replica from the folder nnfit/usr_path to + just nnfit/replica_ + + Examples + -------- + Starting from the replica 5 path + initial_replica_file_path = "/myfolder/myfit/nnfit/myfit/myfit_5.dat" + and using this function with usr_path = "/myfolder/myfit", one gets + final_replica_file_path = "/myfolder/myfit/nnfit/replica_5/myfit.dat" + """ + nnfit = usr_path / "nnfit" + replica_file_path = nnfit / usr_path.stem / f"{usr_path.stem}_{replica_num:04d}.dat" + dest_path_replica = nnfit / f"replica_{replica_num}" / f"{usr_path.stem}.dat" + shutil.move(replica_file_path, dest_path_replica) + + +def check_is_a_fit(config_folder): + """Check if config_folder is a fit folder, i.e. if it contains the filter.yml file + and the nnfit folder.""" + usr_path = pathlib.Path(config_folder) + filter_path = usr_path / "filter.yml" + if not filter_path.is_file(): + raise ValueError("filter.yaml file not found: the path" + str(filter_path.absolute()) + " is not valid") + nnfitpath = usr_path / "nnfit" + if not nnfitpath.is_dir(): + raise ValueError("nnfit folder not found: provided path is not valid") diff --git a/n3fit/src/n3fit/io/writer.py b/n3fit/src/n3fit/io/writer.py index 9eb81068ce..a670a9e365 100644 --- a/n3fit/src/n3fit/io/writer.py +++ b/n3fit/src/n3fit/io/writer.py @@ -11,7 +11,7 @@ import validphys import n3fit from n3fit import vpinterface - +from evolven3fit_new.cli import XGRID class WriterWrapper: def __init__(self, replica_number, pdf_object, stopping_object, q2, timings): @@ -249,8 +249,8 @@ def storefit( q_0^2 """ # build exportgrid - xgrid = np.array([1.00000000000000e-09, 1.29708482343957e-09, 1.68242903474257e-09, 2.18225315420583e-09, 2.83056741739819e-09, 3.67148597892941e-09, 4.76222862935315e-09, 6.17701427376180e-09, 8.01211109898438e-09, 1.03923870607245e-08, 1.34798064073805e-08, 1.74844503691778e-08, 2.26788118881103e-08, 2.94163370300835e-08, 3.81554746595878e-08, 4.94908707232129e-08, 6.41938295708371e-08, 8.32647951986859e-08, 1.08001422993829e-07, 1.40086873081130e-07, 1.81704331793772e-07, 2.35685551545377e-07, 3.05703512595323e-07, 3.96522309841747e-07, 5.14321257236570e-07, 6.67115245136676e-07, 8.65299922973143e-07, 1.12235875241487e-06, 1.45577995547683e-06, 1.88824560514613e-06, 2.44917352454946e-06, 3.17671650028717e-06, 4.12035415232797e-06, 5.34425265752090e-06, 6.93161897806315e-06, 8.99034258238145e-06, 1.16603030112258e-05, 1.51228312288769e-05, 1.96129529349212e-05, 2.54352207134502e-05, 3.29841683435992e-05, 4.27707053972016e-05, 5.54561248105849e-05, 7.18958313632514e-05, 9.31954227979614e-05, 1.20782367731330e-04, 1.56497209466554e-04, 2.02708936328495e-04, 2.62459799331951e-04, 3.39645244168985e-04, 4.39234443000422e-04, 5.67535660104533e-04, 7.32507615725537e-04, 9.44112105452451e-04, 1.21469317686978e-03, 1.55935306118224e-03, 1.99627451141338e-03, 2.54691493736552e-03, 3.23597510213126e-03, 4.09103436509565e-03, 5.14175977083962e-03, 6.41865096062317e-03, 7.95137940306351e-03, 9.76689999624100e-03, 1.18876139251364e-02, 1.43298947643919e-02, 1.71032279460271e-02, 2.02100733925079e-02, 2.36463971369542e-02, 2.74026915728357e-02, 3.14652506132444e-02, 3.58174829282429e-02, 4.04411060163317e-02, 4.53171343973807e-02, 5.04266347950069e-02, 5.57512610084339e-02, 6.12736019390519e-02, 6.69773829498255e-02, 7.28475589986517e-02, 7.88703322292727e-02, 8.50331197801452e-02, 9.13244910278679e-02, 9.77340879783772e-02, 1.04252538208639e-01, 1.10871366547237e-01, 1.17582909372878e-01, 1.24380233801599e-01, 1.31257062945031e-01, 1.38207707707289e-01, 1.45227005135651e-01, 1.52310263065985e-01, 1.59453210652156e-01, 1.66651954293987e-01, 1.73902938455578e-01, 1.81202910873333e-01, 1.88548891679097e-01, 1.95938145999193e-01, 2.03368159629765e-01, 2.10836617429103e-01, 2.18341384106561e-01, 2.25880487124065e-01, 2.33452101459503e-01, 2.41054536011681e-01, 2.48686221452762e-01, 2.56345699358723e-01, 2.64031612468684e-01, 2.71742695942783e-01, 2.79477769504149e-01, 2.87235730364833e-01, 2.95015546847664e-01, 3.02816252626866e-01, 3.10636941519503e-01, 3.18476762768082e-01, 3.26334916761672e-01, 3.34210651149156e-01, 3.42103257303627e-01, 3.50012067101685e-01, 3.57936449985571e-01, 3.65875810279643e-01, 3.73829584735962e-01, 3.81797240286494e-01, 3.89778271981947e-01, 3.97772201099286e-01, 4.05778573402340e-01, 4.13796957540671e-01, 4.21826943574548e-01, 4.29868141614175e-01, 4.37920180563205e-01, 4.45982706956990e-01, 4.54055383887562e-01, 4.62137890007651e-01, 4.70229918607142e-01, 4.78331176755675e-01, 4.86441384506059e-01, 4.94560274153348e-01, 5.02687589545177e-01, 5.10823085439086e-01, 5.18966526903235e-01, 5.27117688756998e-01, 5.35276355048428e-01, 5.43442318565661e-01, 5.51615380379768e-01, 5.59795349416641e-01, 5.67982042055800e-01, 5.76175281754088e-01, 5.84374898692498e-01, 5.92580729444440e-01, 6.00792616663950e-01, 6.09010408792398e-01, 6.17233959782450e-01, 6.25463128838069e-01, 6.33697780169485e-01, 6.41937782762089e-01, 6.50183010158361e-01, 6.58433340251944e-01, 6.66688655093089e-01, 6.74948840704708e-01, 6.83213786908386e-01, 6.91483387159697e-01, 6.99757538392251e-01, 7.08036140869916e-01, 7.16319098046733e-01, 7.24606316434025e-01, 7.32897705474271e-01, 7.41193177421404e-01, 7.49492647227008e-01, 7.57796032432224e-01, 7.66103253064927e-01, 7.74414231541921e-01, 7.82728892575836e-01, 7.91047163086478e-01, 7.99368972116378e-01, 8.07694250750291e-01, 8.16022932038457e-01, 8.24354950923382e-01, 8.32690244169987e-01, 8.41028750298844e-01, 8.49370409522600e-01, 8.57715163684985e-01, 8.66062956202683e-01, 8.74413732009721e-01, 8.82767437504206e-01, 8.91124020497459e-01, 8.99483430165226e-01, 9.07845617001021e-01, 9.16210532771399e-01, 9.24578130473112e-01, 9.32948364292029e-01, 9.41321189563734e-01, 9.49696562735755e-01, 9.58074441331298e-01, 9.66454783914439e-01, 9.74837550056705e-01, 9.83222700304978e-01, 9.91610196150662e-01, 1.00000000000000e+00]).reshape(-1, 1) - + xgrid = XGRID.reshape(-1, 1) + result = pdf_object(xgrid, flavours="n3fit").squeeze() lha = evln2lha(result.T).T diff --git a/n3fit/src/n3fit/scripts/evolven3fit_new.py b/n3fit/src/n3fit/scripts/evolven3fit_new.py new file mode 100644 index 0000000000..a1bfe29370 --- /dev/null +++ b/n3fit/src/n3fit/scripts/evolven3fit_new.py @@ -0,0 +1,159 @@ +""" +This module contains the CLI for evolven3fit_new +""" +import logging +import pathlib +from argparse import ArgumentParser +import sys +import numpy as np + +from evolven3fit_new import eko_utils, evolve, cli + +_logger = logging.getLogger(__name__) + + +def construct_eko_parser(subparsers): + parser = subparsers.add_parser( + "produce_eko", + help=( + """Produce the eko for the specified theory. + 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. + The x_grid starts at x_grid_ini and ends at 1.0 and contains the + 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", + ) + return parser + + +def construct_evolven3fit_parser(subparsers): + parser = subparsers.add_parser( + "evolve", + 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", + ) + parser.add_argument( + "-l", + "--load", + type=pathlib.Path, + default=None, + help="Path of the EKO to be loaded", + ) + parser.add_argument( + "-d", + "--dump", + type=pathlib.Path, + default=None, + help="Path where the EKO is dumped", + ) + parser.add_argument( + "-f", + "--force", + action="store_true", + help="Force the evolution to be done even if it has already been done", + ) + return parser + + +def main(): + parser = ArgumentParser( + description="evolven3fit_new - 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", + ) + subparsers = parser.add_subparsers(title="actions", dest="actions") + eko_parser = construct_eko_parser(subparsers) + evolven3fit_parser = construct_evolven3fit_parser(subparsers) + + args = parser.parse_args() + op_card_info = { + "ev_op_max_order": 10, + "ev_op_iterations": 1, + "n_integration_cores": args.n_cores, + "backward_inversion": "expanded", + } + theory_card_info = {} + if args.actions == "evolve": + cli.cli_evolven3fit_new( + args.configuration_folder, + args.q_fin, + args.q_points, + op_card_info, + theory_card_info, + args.dump, + args.load, + args.force, + ) + elif args.actions == "produce_eko": + stdout_log = logging.StreamHandler(sys.stdout) + stdout_log.setLevel(evolve.LOGGING_SETTINGS["level"]) + stdout_log.setFormatter(evolve.LOGGING_SETTINGS["formatter"]) + for logger_ in (_logger, *[logging.getLogger("eko")]): + logger_.handlers = [] + logger_.setLevel(evolve.LOGGING_SETTINGS["level"]) + logger_.addHandler(stdout_log) + if args.x_grid_ini is None: + if args.x_grid_points is None: + x_grid = cli.XGRID + else: + raise ValueError( + "x_grid_ini and x_grid_points must be specified either both or none of them" + ) + elif args.x_grid_points is None: + raise ValueError( + "x_grid_ini and x_grid_points must be specified either both or none of them" + ) + else: + x_grid = np.geomspace(args.x_grid_ini, 1.0, args.x_grid_points) + tcard, opcard = eko_utils.construct_eko_cards( + args.theoryID, args.q_fin, args.q_points, x_grid, op_card_info, theory_card_info + ) + _ = eko_utils.construct_eko_for_fit(tcard, opcard, args.dump) + + +if __name__ == "__main__": + main() diff --git a/n3fit/src/n3fit/tests/regressions/filter.yml b/n3fit/src/n3fit/tests/regressions/filter.yml new file mode 100644 index 0000000000..efb09d94e3 --- /dev/null +++ b/n3fit/src/n3fit/tests/regressions/filter.yml @@ -0,0 +1,88 @@ +# +# Configuration file for n3fit regression tests +# This runcard includes two DIS datasets, one Hadronic dataset +# and two positivity datasets +# + +############################################################ +description: n3fit regression test + +############################################################ +# frac: training fraction +# ewk: apply ewk k-factors +# sys: systematics treatment (see systypes) +dataset_inputs: +- { dataset: NMC, frac: 0.5 } +- { dataset: SLACP, frac: 0.5} +- { dataset: CMSZDIFF12, frac: 0.5, cfac: ['QCD'], sys: 10 } +- { dataset: CMSTTBARTOT, frac: 1.0, cfac: ['QCD'] } + +############################################################ +datacuts: + t0pdfset: NNPDF40_nnlo_as_01180 # PDF set to generate t0 covmat + q2min : 3.49 # Q2 minimum + w2min : 12.5 # W2 minimum + combocuts : NNPDF31 # NNPDF3.0 final kin. cuts + jetptcut_tev : 0 # jet pt cut for tevatron + jetptcut_lhc : 0 # jet pt cut for lhc + wptcut_lhc : 30.0 # Minimum pT for W pT diff distributions + jetycut_tev : 1e30 # jet rap. cut for tevatron + jetycut_lhc : 1e30 # jet rap. cut for lhc + dymasscut_min: 0 # dy inv.mass. min cut + dymasscut_max: 1e30 # dy inv.mass. max cut + jetcfactcut : 1e30 # jet cfact. cut + +############################################################ +theory: + theoryid: 162 # database id + +############################################################ +genrep: True # True = generate MC replicas, False = use real data +trvlseed: 3 +nnseed: 2 +mcseed: 1 + +load: "weights.h5" + +parameters: # This defines the parameter dictionary that is passed to the Model Trainer + nodes_per_layer: [15, 10, 8] + activation_per_layer: ['sigmoid', 'sigmoid', 'linear'] + initializer: 'glorot_normal' + optimizer: + optimizer_name: 'RMSprop' + learning_rate: 0.00001 + clipnorm: 1.0 + epochs: 1000 + positivity: + multiplier: 1.05 + initial: 1.5 + stopping_patience: 0.30 # percentage of the number of epochs + layer_type: 'dense' + dropout: 0.0 + threshold_chi2: 10.0 + +fitting: + fitbasis: NN31IC # EVOL (7), EVOLQED (8), etc. + basis: + - { fl: sng, smallx: [1.05,1.19], largex: [1.47,2.70] } + - { fl: g, smallx: [0.94,1.25], largex: [0.11,5.87] } + - { fl: v, smallx: [0.54,0.75], largex: [1.15,2.76] } + - { fl: v3, smallx: [0.21,0.57], largex: [1.35,3.08] } + - { fl: v8, smallx: [0.52,0.76], largex: [0.77,3.56] } + - { fl: t3, smallx: [-0.37,1.52], largex: [1.74,3.39] } + - { fl: t8, smallx: [0.56,1.29], largex: [1.45,3.03] } + - { fl: cp, smallx: [0.12,1.19], largex: [1.83,6.70] } + +############################################################ +positivity: + posdatasets: + - { dataset: POSF2U, maxlambda: 1e6 } # Positivity Lagrange Multiplier + - { dataset: POSDYC, maxlambda: 1e5 } + +integrability: + integdatasets: + - {dataset: INTEGXT8, maxlambda: 1e2} + - {dataset: INTEGXT3, maxlambda: 1e2} + +############################################################ +debug: true diff --git a/n3fit/src/n3fit/tests/test_evolven3fit.py b/n3fit/src/n3fit/tests/test_evolven3fit.py new file mode 100644 index 0000000000..278d7212e3 --- /dev/null +++ b/n3fit/src/n3fit/tests/test_evolven3fit.py @@ -0,0 +1,70 @@ +import pathlib +import logging +from numpy.testing import assert_allclose +import numpy as np +from validphys.pdfbases import PIDS_DICT +from evolven3fit_new import utils, eko_utils + +REGRESSION_FOLDER = pathlib.Path(__file__).with_name("regressions") +log = logging.getLogger(__name__) + +def check_consecutive_members(grid, value): + """Check if the first occurrence of value in grid is followed by value again""" + return np.allclose(grid[list(grid).index(value)+1], value) + +def test_utils(): + #Testing the default grid + grid = utils.generate_q2grid(1.65, None, None, {}) + assert_allclose(1.65**2, grid[0]) + assert len(grid) == 50 + # We expect the bottom mass to be repeated twice because it is intended once in 4 flavors and once in 5 flavors. + assert check_consecutive_members(grid, 4.92**2) + #Testing if the points of the matching are correctly repeated twice + matched_grid = utils.generate_q2grid(1.65, 1.0e5, 100, {4.92:2.0, 100:1.0} ) + assert len(matched_grid) == 100 + assert_allclose((1.0e5)**2, matched_grid[-1]) + assert check_consecutive_members(matched_grid, (4.92*2.0)**2) + assert check_consecutive_members(matched_grid, (100.*1.0)**2) + #Testing the fake LHAPDF class + q20 = 1.65**2 + x_grid = np.geomspace(1.0e-7, 1.0, 30) + fake_grids = [[x*(1.-x) for x in x_grid] for pid in PIDS_DICT.keys()] + pdf_grid = dict([(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: + assert_allclose(my_PDF.xfxQ2(pid, x, q20), x*(1.-x) ) + #Testing read_runcard + runcard = utils.read_runcard(REGRESSION_FOLDER) + assert runcard["description"] == "n3fit regression test" + assert runcard["datacuts"]["t0pdfset"] == "NNPDF40_nnlo_as_01180" + #Testing get_theoryID_from_runcard + ID = utils.get_theoryID_from_runcard(REGRESSION_FOLDER) + assert ID == 162 + +def test_eko_utils(): + #Testing construct eko cards + theoryID = 162 + q_fin = 100 + q_points = 5 + x_grid = [1.e-3, 0.1, 1.0] + pto = 2 + qref = 91.2 + comments = "Test" + n_cores = 6 + t_card, op_card = eko_utils.construct_eko_cards(theoryID, q_fin, q_points, x_grid, {'n_integration_cores' : n_cores, 'interpolation_polynomial_degree' : 2}, {'Comments' : comments}) + assert t_card['Qref'] == qref + assert t_card['PTO'] == pto + assert t_card['Comments'] == comments + assert op_card['n_integration_cores'] == n_cores + assert_allclose(op_card["interpolation_xgrid"], x_grid) + assert_allclose(op_card["Q2grid"][0], 1.65**2) + assert_allclose(op_card["Q2grid"][-1], q_fin**2) + #In this case there are not enough points to have twice the bottom matching scale + assert_allclose(op_card["Q2grid"][1], 4.92**2) + #Testing construct_eko_for_fit + eko_op = eko_utils.construct_eko_for_fit(t_card ,op_card) + assert_allclose(eko_op['interpolation_xgrid'], x_grid) + assert_allclose(list(eko_op['Q2grid']), op_card["Q2grid"]) diff --git a/validphys2/src/validphys/pdfbases.py b/validphys2/src/validphys/pdfbases.py index 6f61c138d0..aac0f01ef1 100644 --- a/validphys2/src/validphys/pdfbases.py +++ b/validphys2/src/validphys/pdfbases.py @@ -35,6 +35,23 @@ (21 , r"g"), )) +PIDS_DICT = { + -6: "TBAR", + -5: "BBAR", + -4: "CBAR", + -3: "SBAR", + -2: "UBAR", + -1: "DBAR", + 21: "GLUON", + 1: "D", + 2: "U", + 3: "S", + 4: "C", + 5: "B", + 6: "T", + 22: "PHT", + } + # Canonical ordering of PDG codes (so flavour basis) ALL_FLAVOURS = (-6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6, 22) DEFAULT_FLARR = (-3,-2,-1,0,1,2,3,4)