diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 5e8b44970c..c2c3028999 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -51,9 +51,9 @@ requirements: - sphinxcontrib-bibtex - curio >=1.0 - pineappl >=0.5.8 - - eko >=0.12.0 + - eko >=0.12.2 - banana-hep >=0.6.8 - + - fiatlux test: requires: diff --git a/doc/sphinx/source/n3fit/runcard_detailed.rst b/doc/sphinx/source/n3fit/runcard_detailed.rst index 180de72667..58aa9b94b6 100644 --- a/doc/sphinx/source/n3fit/runcard_detailed.rst +++ b/doc/sphinx/source/n3fit/runcard_detailed.rst @@ -419,4 +419,5 @@ It is however possible to disable them by setting to false the ``sum_rules`` fla It is also possible to impose just the valence or the momentum sum rules by using the -``VSR`` or ``MSR`` flags, respectively (``True`` is equal to ``All``). \ No newline at end of file +``VSR`` or ``MSR`` flags, respectively (``True`` is equal to ``All``). + diff --git a/doc/sphinx/source/references.bib b/doc/sphinx/source/references.bib index 3dc9bb8782..bb646dcae4 100644 --- a/doc/sphinx/source/references.bib +++ b/doc/sphinx/source/references.bib @@ -604,3 +604,32 @@ @article{Carrazza:2016htc pages = "205", year = "2016" } + +@article{Manohar:2017eqh, + author = "Manohar, Aneesh V. and Nason, Paolo and Salam, Gavin P. and Zanderighi, Giulia", + title = "{The Photon Content of the Proton}", + eprint = "1708.01256", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "CERN-TH-2017-141", + doi = "10.1007/JHEP12(2017)046", + journal = "JHEP", + volume = "12", + pages = "046", + year = "2017" +} + +@article{Manohar:2016nzj, + author = "Manohar, Aneesh and Nason, Paolo and Salam, Gavin P. and Zanderighi, Giulia", + title = "{How bright is the proton? A precise determination of the photon parton distribution function}", + eprint = "1607.04266", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "CERN-TH-2016-155", + doi = "10.1103/PhysRevLett.117.242002", + journal = "Phys. Rev. Lett.", + volume = "117", + number = "24", + pages = "242002", + year = "2016" +} diff --git a/doc/sphinx/source/tutorials/run-fit.md b/doc/sphinx/source/tutorials/run-fit.md index dc26ec8d94..f0a1d10201 100644 --- a/doc/sphinx/source/tutorials/run-fit.md +++ b/doc/sphinx/source/tutorials/run-fit.md @@ -256,3 +256,8 @@ Iterate the fit It may be desirable to iterate a fit to achieve a higher degree of convergence/stability in the fit. To read more about this, see [How to run an iterated fit](run-iterated-fit). + +QED fit +------- + +In order to run a QED fit see [How to run a QED fit](run-qed-fit) diff --git a/doc/sphinx/source/tutorials/run-qed-fit.rst b/doc/sphinx/source/tutorials/run-qed-fit.rst new file mode 100644 index 0000000000..051694cc8d --- /dev/null +++ b/doc/sphinx/source/tutorials/run-qed-fit.rst @@ -0,0 +1,30 @@ +.. _run-qed-fit: + +========================== +How to run a QED fit +========================== + +It is possible to perform a QED fit adding the key `fiatlux` to the runcard. In this way +a photon PDF will be generated using the `FiatLux` public library that implements the `LuxQED` +(see :cite:p:`Manohar:2016nzj` and :cite:p:`Manohar:2017eqh`) approach. +The parameters to be added are the following: + +.. code-block:: yaml + + fiatlux: + luxset: NNPDF40_nnlo_as_01180 + additional_errors: true + luxseed: 1234567890 + +`luxset` is the PDF set used to generate the photon PDF with `FiatLux `. +The code will generate as many photon replicas as the number of replicas contained in the `luxset`. Therefore, if the user +tries to generate a replica with ID higher than the maximum ID of the `luxset`, the code will +raise an error. Moreover, being the `LuxQED` approach an iterated prcedure, and given that some replicas +do not pass the `postfit` selection criteria, the user should make sure that the number of replicas in +the `luxset` is high enough such that in the final iteration there will be a number of replicas +higher than the final replicas desired. +`additional_errors` is the parameter that switches on and off the additional errors of the `LuxQED` approach, +while `luxseed` is the seed used to generate such errors. +This errors should be switched on only in the very last iteration of the procedure. + +Whenever the photon PDF is generated, it will remain constant during the fit and will enter in the `MSR`. diff --git a/n3fit/runcards/examples/Basic_runcard_qed.yml b/n3fit/runcards/examples/Basic_runcard_qed.yml new file mode 100644 index 0000000000..309f5901d5 --- /dev/null +++ b/n3fit/runcards/examples/Basic_runcard_qed.yml @@ -0,0 +1,171 @@ +# +# Configuration file for n3fit +# +############################################################ +description: Basic runcard qed + +############################################################ +# frac: training fraction +# ewk: apply ewk k-factors +# sys: systematics treatment (see systypes) +dataset_inputs: +# - {dataset: NMCPD_dw_ite, frac: 0.75} +- {dataset: NMC, frac: 0.75} +- {dataset: SLACP_dwsh, frac: 0.75} +# - {dataset: SLACD_dw_ite, frac: 0.75} +# - {dataset: BCDMSP_dwsh, frac: 0.75} +# - {dataset: BCDMSD_dw_ite, frac: 0.75} +# - {dataset: CHORUSNUPb_dw_ite, frac: 0.75} +# - {dataset: CHORUSNBPb_dw_ite, frac: 0.75} +# - {dataset: NTVNUDMNFe_dw_ite, frac: 0.75, cfac: [MAS]} +# - {dataset: NTVNBDMNFe_dw_ite, frac: 0.75, cfac: [MAS]} +# - {dataset: HERACOMBNCEM, frac: 0.75} +# - {dataset: HERACOMBNCEP460, frac: 0.75} +# - {dataset: HERACOMBNCEP575, frac: 0.75} +# - {dataset: HERACOMBNCEP820, frac: 0.75} +# - {dataset: HERACOMBNCEP920, frac: 0.75} +# - {dataset: HERACOMBCCEM, frac: 0.75} +# - {dataset: HERACOMBCCEP, frac: 0.75} +# - {dataset: HERACOMB_SIGMARED_C, frac: 0.75} +# - {dataset: HERACOMB_SIGMARED_B, frac: 0.75} +# - {dataset: DYE886R_dw_ite, frac: 0.75, cfac: [QCD]} +# - {dataset: DYE886P, frac: 0.75, cfac: [QCD]} +# - {dataset: DYE605_dw_ite, frac: 0.75, cfac: [QCD]} +# - {dataset: DYE906R_dw_ite, frac: 0.75, cfac: [ACC, QCD]} +# - {dataset: CDFZRAP_NEW, frac: 0.75, cfac: [QCD]} +# - {dataset: D0ZRAP_40, frac: 0.75, cfac: [QCD]} +# - {dataset: D0WMASY, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASWZRAP36PB, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASZHIGHMASS49FB, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASLOMASSDY11EXT, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASWZRAP11CC, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASWZRAP11CF, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASDY2D8TEV, frac: 0.75, cfac: [QCDEWK]} +# - {dataset: ATLAS_DY_2D_8TEV_LOWMASS, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_WZ_TOT_13TEV, frac: 0.75, cfac: [NRM, QCD]} +# - {dataset: ATLAS_WP_JET_8TEV_PT, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_WM_JET_8TEV_PT, frac: 0.75, cfac: [QCD]} +- {dataset: ATLASZPT8TEVMDIST, frac: 0.75, cfac: [QCD], sys: 10} +# - {dataset: ATLASZPT8TEVYDIST, frac: 0.75, cfac: [QCD], sys: 10} +# - {dataset: ATLASTTBARTOT7TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASTTBARTOT8TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_TTBARTOT_13TEV_FULLLUMI, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_TTB_DIFF_8TEV_LJ_TRAPNORM, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_TTB_DIFF_8TEV_LJ_TTRAPNORM, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_TOPDIFF_DILEPT_8TEV_TTRAPNORM, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_1JET_8TEV_R06_DEC, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_2JET_7TEV_R06, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLASPHT15_SF, frac: 0.75, cfac: [QCD, EWK]} +# - {dataset: ATLAS_SINGLETOP_TCH_R_7TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_SINGLETOP_TCH_R_13TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_7TEV_T_RAP_NORM, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_7TEV_TBAR_RAP_NORM, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_8TEV_T_RAP_NORM, frac: 0.75, cfac: [QCD]} +# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_8TEV_TBAR_RAP_NORM, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSWEASY840PB, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSWMASY47FB, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSDY2D11, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSWMU8TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSZDIFF12, frac: 0.75, cfac: [QCD, NRM], sys: 10} +# - {dataset: CMS_2JET_7TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_1JET_8TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSTTBARTOT7TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSTTBARTOT8TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSTTBARTOT13TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSTOPDIFF8TEVTTRAPNORM, frac: 0.75, cfac: [QCD]} +# - {dataset: CMSTTBARTOT5TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_TTBAR_2D_DIFF_MTT_TRAP_NORM, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_TTB_DIFF_13TEV_2016_2L_TRAP, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_TTB_DIFF_13TEV_2016_LJ_TRAP, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_SINGLETOP_TCH_TOT_7TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_SINGLETOP_TCH_R_8TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: CMS_SINGLETOP_TCH_R_13TEV, frac: 0.75, cfac: [QCD]} +# - {dataset: LHCBZ940PB, frac: 0.75, cfac: [QCD]} +# - {dataset: LHCBZEE2FB_40, frac: 0.75, cfac: [QCD]} +# - {dataset: LHCBWZMU7TEV, frac: 0.75, cfac: [NRM, QCD]} +# - {dataset: LHCBWZMU8TEV, frac: 0.75, cfac: [NRM, QCD]} +# - {dataset: LHCB_Z_13TEV_DIMUON, frac: 0.75, cfac: [QCD]} +# - {dataset: LHCB_Z_13TEV_DIELECTRON, frac: 0.75, cfac: [QCD]} + +############################################################ +datacuts: + t0pdfset: NNPDF40_nnlo_as_01180 # PDF set to generate t0 covmat + q2min: 3.49 # Q2 minimum + w2min: 12.5 # W2 minimum + +############################################################ +theory: + theoryid: 522 # database id + +############################################################ +trvlseed: 1551864071 +nnseed: 676150632 +mcseed: 619859729 +save: false +genrep: true # true = generate MC replicas, false = use real data + +parameters: # This defines the parameter dictionary that is passed to the Model Trainer + nodes_per_layer: [25, 20, 8] + activation_per_layer: [tanh, tanh, linear] + initializer: glorot_normal + optimizer: + clipnorm: 6.073e-6 + learning_rate: 2.621e-3 + optimizer_name: Nadam + epochs: 900 + positivity: + initial: 184.8 + multiplier: + integrability: + initial: 10 + multiplier: + stopping_patience: 0.1 + layer_type: dense + dropout: 0.0 + threshold_chi2: 3.5 + +fitting: + fitbasis: EVOL # EVOL (7), EVOLQED (8), etc. + basis: + - {fl: sng, trainable: false, smallx: [1.089, 1.117], largex: [1.462, 3.008]} + - {fl: g, trainable: false, smallx: [0.7542, 1.105], largex: [2.826, 5.407]} + - {fl: v, trainable: false, smallx: [0.4715, 0.7253], largex: [1.564, 3.48]} + - {fl: v3, trainable: false, smallx: [0.1372, 0.4205], largex: [1.755, 3.451]} + - {fl: v8, trainable: false, smallx: [0.5641, 0.7702], largex: [1.513, 3.433]} + - {fl: t3, trainable: false, smallx: [-0.4942, 0.9992], largex: [1.751, 3.383]} + - {fl: t8, trainable: false, smallx: [0.532, 0.8572], largex: [1.541, 3.349]} + - {fl: t15, trainable: false, smallx: [1.052, 1.14], largex: [1.487, 3.09]} + +############################################################ +positivity: + posdatasets: + - {dataset: POSF2U, maxlambda: 1e6} # Positivity Lagrange Multiplier + # - {dataset: POSF2DW, maxlambda: 1e6} + # - {dataset: POSF2S, maxlambda: 1e6} + - {dataset: POSFLL, maxlambda: 1e6} + # - {dataset: POSDYU, maxlambda: 1e10} + # - {dataset: POSDYD, maxlambda: 1e10} + # - {dataset: POSDYS, maxlambda: 1e10} + # - {dataset: POSF2C, maxlambda: 1e6} + # - {dataset: POSXUQ, maxlambda: 1e6} # Positivity of MSbar PDFs + # - {dataset: POSXUB, maxlambda: 1e6} + # - {dataset: POSXDQ, maxlambda: 1e6} + # - {dataset: POSXDB, maxlambda: 1e6} + # - {dataset: POSXSQ, maxlambda: 1e6} + # - {dataset: POSXSB, maxlambda: 1e6} + # - {dataset: POSXGL, maxlambda: 1e6} + +############################################################ +integrability: + integdatasets: + # - {dataset: INTEGXT8, maxlambda: 1e2} + - {dataset: INTEGXT3, maxlambda: 1e2} + +############################################################ +debug: True +maxcores: 8 + +fiatlux: + luxset: NNPDF40_nnlo_as_01180 + additional_errors: true # should be set to true only for the last iteration + luxseed: 1234567890 diff --git a/n3fit/src/evolven3fit_new/eko_utils.py b/n3fit/src/evolven3fit_new/eko_utils.py index 17467dc1d9..0dc4872979 100644 --- a/n3fit/src/evolven3fit_new/eko_utils.py +++ b/n3fit/src/evolven3fit_new/eko_utils.py @@ -10,7 +10,7 @@ _logger = logging.getLogger(__name__) -EVOLVEN3FIT_CONFIGS_DEFAULTS = { +EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN = { "ev_op_iterations": 1, "ev_op_max_order": (1, 0), "evolution_method": "truncated", @@ -18,6 +18,14 @@ "n_integration_cores": 1, } +EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = { + "ev_op_iterations": 10, + "ev_op_max_order": (1, 0), + "evolution_method": "iterate-exact", + "inversion_method": "exact", + "n_integration_cores": 1, +} + NFREF_DEFAULT = 5 NF0_DEFAULT = 4 @@ -68,7 +76,10 @@ def construct_eko_cards( ) op_card["rotations"]["xgrid"] = x_grid # Specific defaults for evolven3fit evolution - op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS) + 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 # Note that every entry that is not a dictionary should not be diff --git a/n3fit/src/evolven3fit_new/evolve.py b/n3fit/src/evolven3fit_new/evolve.py index f9141ddb2b..08f139e8e7 100644 --- a/n3fit/src/evolven3fit_new/evolve.py +++ b/n3fit/src/evolven3fit_new/evolve.py @@ -6,8 +6,7 @@ import eko from eko import basis_rotation from eko import runner -from ekobox import genpdf, info_file -from ekomark import apply +from ekobox import genpdf, info_file, apply from reportengine.compat import yaml from validphys.loader import Loader @@ -69,7 +68,7 @@ def evolve_fit( 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"])) + log.setFormatter(LOGGING_SETTINGS["formatter"]) for logger in (_logger, *[logging.getLogger("eko")]): logger.handlers = [] logger.setLevel(LOGGING_SETTINGS["level"]) @@ -85,6 +84,7 @@ def evolve_fit( theory, op = eko_utils.construct_eko_cards( theoryID, q_fin, q_points, x_grid, op_card_dict, theory_card_dict ) + qed = theory.order[1] > 0 if eko_path is not None: eko_path = pathlib.Path(eko_path) _logger.info(f"Loading eko from : {eko_path}") @@ -105,10 +105,12 @@ def evolve_fit( info["XMin"] = float(x_grid[0]) info["XMax"] = float(x_grid[-1]) with eko.EKO.read(eko_path) as eko_op: - info["AlphaS_Qs"] = eko_op.mu2grid.tolist() + if eko.__version__ >= "0.13": + raise ModuleNotFoundError("Please, fix evolven3fit np.sqrt(Q) hack") + info["AlphaS_Qs"] = np.sqrt(info["AlphaS_Qs"]).tolist() 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) + evolved_block = evolve_exportgrid(initial_PDFs_dict[replica], eko_op, x_grid, qed) dump_evolved_replica( evolved_block, usr_path, int(replica.removeprefix("replica_")) ) @@ -142,7 +144,7 @@ def load_fit(usr_path): return pdf_dict -def evolve_exportgrid(exportgrid, eko, x_grid): +def evolve_exportgrid(exportgrid, eko, x_grid, qed): """ Evolves the provided exportgrid for the desired replica with the eko and returns the evolved block @@ -154,6 +156,8 @@ def evolve_exportgrid(exportgrid, eko, x_grid): eko operator for evolution xgrid: list xgrid to be used as the targetgrid + qed: bool + whether qed is activated or not Returns ------- : np.array @@ -163,7 +167,7 @@ def evolve_exportgrid(exportgrid, eko, x_grid): 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) + evolved_pdf = apply.apply_pdf(eko, pdf_to_evolve, qed=qed) # generate block to dump targetgrid = eko.rotations.targetgrid.tolist() diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py index b651f512c5..dc9c6a9d22 100644 --- a/n3fit/src/n3fit/checks.py +++ b/n3fit/src/n3fit/checks.py @@ -1,14 +1,17 @@ """ This module contains checks to be perform by n3fit on the input """ -import os import logging import numbers +import os + import numpy as np -from reportengine.checks import make_argcheck, CheckError -from validphys.pdfbases import check_basis + from n3fit.hyper_optimization import penalties as penalties_module from n3fit.hyper_optimization import rewards as rewards_module +from reportengine.checks import CheckError, make_argcheck +from validphys.core import PDF +from validphys.pdfbases import check_basis log = logging.getLogger(__name__) @@ -393,3 +396,15 @@ def check_deprecated_options(fitting): for option in nnfit_options: if option in fitting: log.warning("'fitting::%s' is an nnfit-only key, it will be ignored", option) + + +@make_argcheck +def check_fiatlux_pdfs_id(replicas, fiatlux, replica_path): + if fiatlux is not None: + luxset = fiatlux["luxset"] + pdfs_ids = luxset.get_members() - 1 # get_members counts also replica0 + max_id = max(replicas) + if max_id > pdfs_ids: + raise CheckError( + f"Cannot generate a photon replica with id larger than the number of replicas of the PDFs set {luxset.name}:\nreplica id={max_id}, replicas of {luxset.name} = {pdfs_ids}" + ) diff --git a/n3fit/src/n3fit/layers/__init__.py b/n3fit/src/n3fit/layers/__init__.py index e5839046cc..a9168bda13 100644 --- a/n3fit/src/n3fit/layers/__init__.py +++ b/n3fit/src/n3fit/layers/__init__.py @@ -1,8 +1,9 @@ -from .preprocessing import Preprocessing -from .rotations import FkRotation, FlavourToEvolution, ObsRotation -from .x_operations import xIntegrator, xDivide -from .msr_normalization import MSR_Normalization +from n3fit.backends import MetaLayer + from .DIS import DIS from .DY import DY from .mask import Mask -from n3fit.backends import MetaLayer +from .msr_normalization import MSR_Normalization +from .preprocessing import Preprocessing +from .rotations import AddPhoton, FkRotation, FlavourToEvolution, ObsRotation +from .x_operations import xDivide, xIntegrator diff --git a/n3fit/src/n3fit/layers/msr_normalization.py b/n3fit/src/n3fit/layers/msr_normalization.py index c31bf8fb41..51219ad73f 100644 --- a/n3fit/src/n3fit/layers/msr_normalization.py +++ b/n3fit/src/n3fit/layers/msr_normalization.py @@ -16,7 +16,8 @@ class MSR_Normalization(MetaLayer): _msr_enabled = False _vsr_enabled = False - def __init__(self, output_dim=14, mode="ALL", **kwargs): + def __init__(self, output_dim=14, mode="ALL", photons_contribution=None, **kwargs): + self._photons = photons_contribution if mode == True or mode.upper() == "ALL": self._msr_enabled = True self._vsr_enabled = True @@ -39,9 +40,9 @@ def __init__(self, output_dim=14, mode="ALL", **kwargs): super().__init__(**kwargs, name="normalizer") - def call(self, pdf_integrated): + def call(self, pdf_integrated, ph_replica): """Imposes the valence and momentum sum rules: - A_g = (1-sigma)/g + A_g = (1-sigma-photon)/g A_v = A_v24 = A_v35 = 3/V A_v3 = 1/V_3 A_v8 = 3/V_8 @@ -52,8 +53,15 @@ def call(self, pdf_integrated): y = op.flatten(pdf_integrated) norm_constants = [] + if self._photons: + photon_integral = self._photons[ph_replica] + else: + photon_integral = 0.0 + if self._msr_enabled: - n_ag = [(1.0 - y[GLUON_IDX[0][0]-1]) / y[GLUON_IDX[0][0]]] * len(GLUON_IDX) + n_ag = [(1.0 - y[GLUON_IDX[0][0] - 1] - photon_integral) / y[GLUON_IDX[0][0]]] * len( + GLUON_IDX + ) norm_constants += n_ag if self._vsr_enabled: diff --git a/n3fit/src/n3fit/layers/rotations.py b/n3fit/src/n3fit/layers/rotations.py index 347571094d..686d3c5bbb 100644 --- a/n3fit/src/n3fit/layers/rotations.py +++ b/n3fit/src/n3fit/layers/rotations.py @@ -2,6 +2,7 @@ This module includes rotation layers """ import numpy as np + from n3fit.backends import MetaLayer from n3fit.backends import operations as op from validphys import pdfbases @@ -28,7 +29,7 @@ def __init__(self, rotation_matrix, axes=1, **kwargs): super().__init__(**kwargs) def is_identity(self): - """ Returns true if the rotation is an identity """ + """Returns true if the rotation is an identity""" # check whether it is a mxm matrix if self.rotation_matrix.shape[0] == self.rotation_matrix.shape[1]: # check whether it is the identity @@ -41,12 +42,15 @@ def call(self, x_raw): class FlavourToEvolution(Rotation): """ - Rotates from the flavour basis to - the evolution basis. + Rotates from the flavour basis to + the evolution basis. """ def __init__( - self, flav_info, fitbasis, **kwargs, + self, + flav_info, + fitbasis, + **kwargs, ): rotation_matrix = pdfbases.fitbasis_to_NN31IC(flav_info, fitbasis) super().__init__(rotation_matrix, axes=1, **kwargs) @@ -92,6 +96,33 @@ def call(self, pdf_raw): return op.batchit(ret) +class AddPhoton(MetaLayer): + """ + Changes the value of the photon component of the PDF to non-zero. + The photon idx in the dimension-14 PDF basis of the FKTables is always index 0. + + In order to avoid bottlenecks, this layer can only compute the photon + for a given fixed shape. + In order to change the shape it is necessary to rebuild the photon. + """ + + def __init__(self, photons, **kwargs): + self._photons_generator = photons + self._pdf_ph = None + super().__init__(**kwargs) + + def register_photon(self, xgrid): + """Compute the photon array and set the layer to be rebuilt""" + if self._photons_generator: + self._pdf_ph = self._photons_generator(xgrid) + self.built = False + + def call(self, pdfs, ph_replica): + if self._pdf_ph is None: + return pdfs + return op.concatenate([self._pdf_ph[ph_replica], pdfs[:, :, 1:]], axis=-1) + + class ObsRotation(MetaLayer): """ Rotation is a layer used to apply a rotation transformation diff --git a/n3fit/src/n3fit/model_gen.py b/n3fit/src/n3fit/model_gen.py index 2814c2f656..ca7844d12a 100644 --- a/n3fit/src/n3fit/model_gen.py +++ b/n3fit/src/n3fit/model_gen.py @@ -10,16 +10,24 @@ """ from dataclasses import dataclass + import numpy as np -from n3fit.msr import msr_impose -from n3fit.layers import DIS, DY, ObsRotation, losses -from n3fit.layers import Preprocessing, FkRotation, FlavourToEvolution -from n3fit.layers.observable import is_unique -from n3fit.backends import MetaModel, Input +from n3fit.backends import Input, Lambda, MetaLayer, MetaModel, base_layer_selector from n3fit.backends import operations as op -from n3fit.backends import MetaLayer, Lambda -from n3fit.backends import base_layer_selector, regularizer_selector +from n3fit.backends import regularizer_selector +from n3fit.layers import ( + DIS, + DY, + AddPhoton, + FkRotation, + FlavourToEvolution, + ObsRotation, + Preprocessing, + losses, +) +from n3fit.layers.observable import is_unique +from n3fit.msr import msr_impose @dataclass @@ -403,6 +411,7 @@ def pdfNN_layer_generator( impose_sumrule=None, scaler=None, parallel_models=1, + photons=None, ): # pylint: disable=too-many-locals """ Generates the PDF model which takes as input a point in x (from 0 to 1) @@ -495,6 +504,10 @@ def pdfNN_layer_generator( will be a (1, None, 2) tensor where dim [:,:,0] is scaled parallel_models: int How many models should be trained in parallel + photon: :py:class:`validphys.photon.compute.Photon` + If given, gives the AddPhoton layer a function to compute a photon which will be added at the + index 0 of the 14-size FK basis + This same function will also be used to compute the MSR component for the photon Returns ------- @@ -560,12 +573,17 @@ def pdfNN_layer_generator( # Evolution layer layer_evln = FkRotation(input_shape=(last_layer_nodes,), output_dim=out) + # Photon layer + layer_photon = AddPhoton(photons=photons) + # Basis rotation basis_rotation = FlavourToEvolution(flav_info=flav_info, fitbasis=fitbasis) # Normalization and sum rules if impose_sumrule: - sumrule_layer, integrator_input = msr_impose(mode=impose_sumrule, scaler=scaler) + sumrule_layer, integrator_input = msr_impose( + mode=impose_sumrule, scaler=scaler, photons=photons + ) model_input["integrator_input"] = integrator_input else: sumrule_layer = lambda x: x @@ -640,7 +658,17 @@ def layer_pdf(x): return layer_evln(layer_fitbasis(x)) # Final PDF (apply normalization) - final_pdf = sumrule_layer(layer_pdf) + normalized_pdf = sumrule_layer(layer_pdf, i) + + # Photon layer, changes the photon from zero to non-zero + def apply_photon(x): + # if photon is None then the photon layer is not applied + if photons: + return layer_photon(normalized_pdf(x), i) + else: + return normalized_pdf(x) + + final_pdf = apply_photon # Create the model pdf_model = MetaModel( diff --git a/n3fit/src/n3fit/model_trainer.py b/n3fit/src/n3fit/model_trainer.py index b756c0e26b..1f1a5a4011 100644 --- a/n3fit/src/n3fit/model_trainer.py +++ b/n3fit/src/n3fit/model_trainer.py @@ -8,18 +8,21 @@ This allows to use hyperscanning libraries, that need to change the parameters of the network between iterations while at the same time keeping the amount of redundant calls to a minimum """ -import logging from collections import namedtuple from itertools import zip_longest +import logging + import numpy as np from scipy.interpolate import PchipInterpolator + from n3fit import model_gen -from n3fit.backends import MetaModel, clear_backend_state, callbacks +from n3fit.backends import MetaModel, callbacks, clear_backend_state from n3fit.backends import operations as op -from n3fit.stopping import Stopping -from n3fit.vpinterface import N3PDF import n3fit.hyper_optimization.penalties import n3fit.hyper_optimization.rewards +from n3fit.stopping import Stopping +from n3fit.vpinterface import N3PDF +from validphys.photon.compute import Photon log = logging.getLogger(__name__) @@ -36,6 +39,7 @@ # See ModelTrainer::_xgrid_generation for the definition of each field and how they are generated InputInfo = namedtuple("InputInfo", ["input", "split", "idx"]) + def _pdf_injection(pdf_layers, observables, masks): """ Takes as input a list of PDF layers each corresponding to one observable (also given as a list) @@ -99,6 +103,9 @@ def __init__( model_file=None, sum_rules=None, parallel_models=1, + theoryid=None, + lux_params=None, + replicas=None, ): """ Parameters @@ -131,6 +138,12 @@ def __init__( whether sum rules should be enabled (All, MSR, VSR, False) parallel_models: int number of models to fit in parallel + theoryid: validphys.core.TheoryIDSpec + object contining info for generating the photon + lux_params: dict + dictionary containing the params needed from LuxQED + replicas_id: list + list with the replicas ids to be fitted """ # Save all input information self.exp_info = exp_info @@ -151,6 +164,9 @@ def __init__( self.all_datasets = [] self._scaler = None self._parallel_models = parallel_models + self.theoryid = theoryid + self.lux_params = lux_params + self.replicas = replicas # Initialise internal variables which define behaviour if debug: @@ -352,9 +368,7 @@ def _xgrid_generation(self): # now the output needs to be splitted so that each experiment takes its corresponding input sp_ar = [[i.shape[1] for i in inputs_unique]] sp_kw = {"axis": 1} - sp_layer = op.as_layer( - op.split, op_args=sp_ar, op_kwargs=sp_kw, name="pdf_split" - ) + sp_layer = op.as_layer(op.split, op_args=sp_ar, op_kwargs=sp_kw, name="pdf_split") return InputInfo(input_layer, sp_layer, inputs_idx) @@ -410,9 +424,7 @@ def _model_generation(self, xinput, pdf_models, partition, partition_idx): for pdf_model in pdf_models: # The input to the full model also works as the input to the PDF model # We apply the Model as Layers and save for later the model (full_pdf) - full_model_input_dict, full_pdf = pdf_model.apply_as_layer( - {"pdf_input": xinput.input} - ) + full_model_input_dict, full_pdf = pdf_model.apply_as_layer({"pdf_input": xinput.input}) all_replicas_pdf.append(full_pdf) # Note that all models share the same symbolic input so we take as input the last @@ -646,6 +658,7 @@ def _generate_pdf( regularizer, regularizer_args, seed, + photons, ): """ Defines the internal variable layer_pdf @@ -672,6 +685,8 @@ def _generate_pdf( dictionary of arguments for the regularizer seed: int seed for the NN + photons: :py:class:`validphys.photon.compute.Photon` + function to compute the photon PDF see model_gen.pdfNN_layer_generator for more information Returns @@ -697,6 +712,7 @@ def _generate_pdf( impose_sumrule=self.impose_sumrule, scaler=self._scaler, parallel_models=self._parallel_models, + photons=photons, ) return pdf_models @@ -867,7 +883,15 @@ def hyperparametrizable(self, params): # Generate the grid in x, note this is the same for all partitions xinput = self._xgrid_generation() - + # Initialize all photon classes for the different replicas: + if self.lux_params: + photons = Photon( + theoryid=self.theoryid, + lux_params=self.lux_params, + replicas=self.replicas, + ) + else: + photons = None ### Training loop for k, partition in enumerate(self.kpartitions): # Each partition of the kfolding needs to have its own separate model @@ -886,8 +910,14 @@ def hyperparametrizable(self, params): params.get("regularizer", None), # regularizer optional params.get("regularizer_args", None), seeds, + photons, ) + if photons: + for m in pdf_models: + pl = m.get_layer("add_photon") + pl.register_photon(xinput.input.tensor_content) + # Model generation joins all the different observable layers # together with pdf model generated above models = self._model_generation(xinput, pdf_models, partition, k) diff --git a/n3fit/src/n3fit/msr.py b/n3fit/src/n3fit/msr.py index e99006de47..9db718f285 100644 --- a/n3fit/src/n3fit/msr.py +++ b/n3fit/src/n3fit/msr.py @@ -2,11 +2,11 @@ The constraint module include functions to impose the momentum sum rules on the PDFs """ import logging + import numpy as np -from n3fit.layers import xDivide, MSR_Normalization, xIntegrator from n3fit.backends import operations as op - +from n3fit.layers import MSR_Normalization, xDivide, xIntegrator log = logging.getLogger(__name__) @@ -36,26 +36,28 @@ def gen_integration_input(nx): return xgrid, weights_array -def msr_impose(nx=int(2e3), mode='All', scaler=None): +def msr_impose(nx=int(2e3), mode='All', scaler=None, photons=None): """ - This function receives: - Generates a function that applies a normalization layer to the fit. - - fit_layer: the 8-basis layer of PDF which we fit - The normalization is computed from the direct output of the NN (so the 7,8-flavours basis) - - final_layer: the 14-basis which is fed to the fktable - and it is applied to the input of the fktable (i.e., to the 14-flavours fk-basis). - It uses pdf_fit to compute the sum rule and returns a modified version of - the final_pdf layer with a normalisation by which the sum rule is imposed - - Parameters - ---------- - nx: int - number of points for the integration grid, default: 2000 - mode: str - what sum rules to compute (MSR, VSR or All), default: All - scaler: scaler - Function to apply to the input. If given the input to the model - will be a (1, None, 2) tensor where dim [:,:,0] is scaled + This function receives: + Generates a function that applies a normalization layer to the fit. + - fit_layer: the 8-basis layer of PDF which we fit + The normalization is computed from the direct output of the NN (so the 7,8-flavours basis) + - final_layer: the 14-basis which is fed to the fktable + and it is applied to the input of the fktable (i.e., to the 14-flavours fk-basis). + It uses pdf_fit to compute the sum rule and returns a modified version of + the final_pdf layer with a normalisation by which the sum rule is imposed + + Parameters + ---------- + nx: int + number of points for the integration grid, default: 2000 + mode: str + what sum rules to compute (MSR, VSR or All), default: All + scaler: scaler + Function to apply to the input. If given the input to the model + will be a (1, None, 2) tensor where dim [:,:,0] is scaled + photon: :py:class:`validphys.photon.compute.Photon` + If given, gives the AddPhoton layer a function to compute the MSR component for the photon """ # 1. Generate the fake input which will be used to integrate @@ -70,23 +72,28 @@ def msr_impose(nx=int(2e3), mode='All', scaler=None): # 3. Now create the integration layer (the layer that will simply integrate, given some weight integrator = xIntegrator(weights_array, input_shape=(nx,)) + # 3.1 If a photon is given, compute the photon component of the MSR + photons_c = None + if photons: + photons_c = photons.integral + # 4. Now create the normalization by selecting the right integrations - normalizer = MSR_Normalization(mode=mode) + normalizer = MSR_Normalization(mode=mode, photons_contribution=photons_c) # 5. Make the xgrid array into a backend input layer so it can be given to the normalization xgrid_input = op.numpy_to_input(xgrid, name="integration_grid") # Finally prepare a function which will take as input the output of the PDF model # and will return it appropiately normalized. - def apply_normalization(layer_pdf): + def apply_normalization(layer_pdf, ph_replica): """ - layer_pdf: output of the PDF, unnormalized, ready for the fktable + layer_pdf: output of the PDF, unnormalized, ready for the fktable """ x_original = op.op_gather_keep_dims(xgrid_input, -1, axis=-1) pdf_integrand = op.op_multiply([division_by_x(x_original), layer_pdf(xgrid_input)]) - normalization = normalizer(integrator(pdf_integrand)) + normalization = normalizer(integrator(pdf_integrand), ph_replica) def ultimate_pdf(x): - return layer_pdf(x)*normalization + return layer_pdf(x) * normalization return ultimate_pdf diff --git a/n3fit/src/n3fit/performfit.py b/n3fit/src/n3fit/performfit.py index a98040dd96..431bb1f5a6 100644 --- a/n3fit/src/n3fit/performfit.py +++ b/n3fit/src/n3fit/performfit.py @@ -5,7 +5,9 @@ # Backend-independent imports import copy import logging + import numpy as np + import n3fit.checks from n3fit.vpinterface import N3PDF @@ -15,14 +17,16 @@ # Action to be called by validphys # All information defining the NN should come here in the "parameters" dict @n3fit.checks.can_run_multiple_replicas +@n3fit.checks.check_fiatlux_pdfs_id def performfit( *, - n3fit_checks_action, # wrapper for all checks - replicas, # checks specific to performfit + n3fit_checks_action, # wrapper for all checks + replicas, # checks specific to performfit replicas_nnseed_fitting_data_dict, posdatasets_fitting_pos_dict, integdatasets_fitting_integ_dict, theoryid, + fiatlux, basis, fitbasis, sum_rules=True, @@ -40,84 +44,86 @@ def performfit( parallel_models=False, ): """ - This action will (upon having read a validcard) process a full PDF fit - for a set of replicas. - - The input to this function is provided by validphys - and/or defined in the runcards or commandline arguments. - - This controller is provided with: - 1. Seeds generated using the replica number and the seeds defined in the runcard. - 2. Loaded datasets with replicas generated. - 2.1 Loaded positivity/integrability sets. - - The workflow of this controller is as follows: - 1. Generate a ModelTrainer object holding information to create the NN and perform a fit - (at this point no NN object has been generated) - 1.1 (if hyperopt) generates the hyperopt scanning dictionary - taking as a base the fitting dictionary and the runcard's hyperscanner dictionary - 2. Pass the dictionary of parameters to ModelTrainer - for the NN to be generated and the fit performed - 2.1 (if hyperopt) Loop over point 4 for `hyperopt` number of times - 3. Once the fit is finished, output the PDF grid and accompanying files - - Parameters - ---------- - genrep: bool - Whether or not to generate MC replicas. (Only used for checks) - data: validphys.core.DataGroupSpec - containing the datasets to be included in the fit. (Only used - for checks) - replicas_nnseed_fitting_data_dict: list[tuple] - list with element for each replica (typically just one) to be - fitted. Each element - is a tuple containing the replica number, nnseed and - ``fitted_data_dict`` containing all of the data, metadata - for each group of datasets which is to be fitted. - posdatasets_fitting_pos_dict: list[dict] - list of dictionaries containing all data and metadata for each - positivity dataset - integdatasets_fitting_integ_dict: list[dict] - list of dictionaries containing all data and metadata for each - integrability dataset - theoryid: validphys.core.TheoryIDSpec - Theory which is used to generate theory predictions from model - during fit. Object also contains some metadata on the theory - settings. - basis: list[dict] - preprocessing information for each flavour to be fitted. - fitbasis: str - Valid basis which the fit is to be ran in. Available bases can - be found in :py:mod:`validphys.pdfbases`. - sum_rules: bool - Whether to impose sum rules in fit. By default set to True - parameters: dict - Mapping containing parameters which define the network - architecture/fitting methodology. - replica_path: pathlib.Path - path to the output of this run - output_path: str - name of the fit - save: None, str - model file where weights will be saved, used in conjunction with - ``load``. - load: None, str - model file from which to load weights from. - hyperscanner: dict - dictionary containing the details of the hyperscanner - hyperopt: int - if given, number of hyperopt iterations to run - kfold_parameters: None, dict - dictionary with kfold settings used in hyperopt. - tensorboard: None, dict - mapping containing tensorboard settings if it is to be used. By - default it is None and tensorboard is not enabled. - debug: bool - activate some debug options - maxcores: int - maximum number of (logical) cores that the backend should be aware of - parallel_models: bool - whether to run models in parallel + This action will (upon having read a validcard) process a full PDF fit + for a set of replicas. + + The input to this function is provided by validphys + and/or defined in the runcards or commandline arguments. + + This controller is provided with: + 1. Seeds generated using the replica number and the seeds defined in the runcard. + 2. Loaded datasets with replicas generated. + 2.1 Loaded positivity/integrability sets. + + The workflow of this controller is as follows: + 1. Generate a ModelTrainer object holding information to create the NN and perform a fit + (at this point no NN object has been generated) + 1.1 (if hyperopt) generates the hyperopt scanning dictionary + taking as a base the fitting dictionary and the runcard's hyperscanner dictionary + 2. Pass the dictionary of parameters to ModelTrainer + for the NN to be generated and the fit performed + 2.1 (if hyperopt) Loop over point 4 for `hyperopt` number of times + 3. Once the fit is finished, output the PDF grid and accompanying files + + Parameters + ---------- + genrep: bool + Whether or not to generate MC replicas. (Only used for checks) + data: validphys.core.DataGroupSpec + containing the datasets to be included in the fit. (Only used + for checks) + replicas_nnseed_fitting_data_dict: list[tuple] + list with element for each replica (typically just one) to be + fitted. Each element + is a tuple containing the replica number, nnseed and + ``fitted_data_dict`` containing all of the data, metadata + for each group of datasets which is to be fitted. + posdatasets_fitting_pos_dict: list[dict] + list of dictionaries containing all data and metadata for each + positivity dataset + integdatasets_fitting_integ_dict: list[dict] + list of dictionaries containing all data and metadata for each + integrability dataset + theoryid: validphys.core.TheoryIDSpec + Theory which is used to generate theory predictions from model + during fit. Object also contains some metadata on the theory + settings. + fiatlux: dict + dictionary containing the params needed from LuxQED + basis: list[dict] + preprocessing information for each flavour to be fitted. + fitbasis: str + Valid basis which the fit is to be ran in. Available bases can + be found in :py:mod:`validphys.pdfbases`. + sum_rules: bool + Whether to impose sum rules in fit. By default set to True + parameters: dict + Mapping containing parameters which define the network + architecture/fitting methodology. + replica_path: pathlib.Path + path to the output of this run + output_path: str + name of the fit + save: None, str + model file where weights will be saved, used in conjunction with + ``load``. + load: None, str + model file from which to load weights from. + hyperscanner: dict + dictionary containing the details of the hyperscanner + hyperopt: int + if given, number of hyperopt iterations to run + kfold_parameters: None, dict + dictionary with kfold settings used in hyperopt. + tensorboard: None, dict + mapping containing tensorboard settings if it is to be used. By + default it is None and tensorboard is not enabled. + debug: bool + activate some debug options + maxcores: int + maximum number of (logical) cores that the backend should be aware of + parallel_models: bool + whether to run models in parallel """ from n3fit.backends import set_initial_state @@ -130,8 +136,8 @@ def performfit( # All potentially backend dependent imports should come inside the fit function # so they can eventually be set from the runcard - from n3fit.model_trainer import ModelTrainer from n3fit.io.writer import WriterWrapper + from n3fit.model_trainer import ModelTrainer # Note: there are three possible scenarios for the loop of replicas: # 1.- Only one replica is being run, in this case the loop is only evaluated once @@ -193,7 +199,10 @@ def performfit( max_cores=maxcores, model_file=load, sum_rules=sum_rules, - parallel_models=n_models + parallel_models=n_models, + theoryid=theoryid, + lux_params=fiatlux, + replicas=replica_idxs, ) # This is just to give a descriptive name to the fit function @@ -285,13 +294,12 @@ def performfit( replica_path_set, output_path.name, training_chi2, val_chi2, exp_chi2 ) log.info( - "Best fit for replica #%d, chi2=%.3f (tr=%.3f, vl=%.3f)", - replica_number, - exp_chi2, - training_chi2, - val_chi2 - ) - + "Best fit for replica #%d, chi2=%.3f (tr=%.3f, vl=%.3f)", + replica_number, + exp_chi2, + training_chi2, + val_chi2, + ) # Save the weights to some file for the given replica if save: diff --git a/n3fit/src/n3fit/scripts/n3fit_exec.py b/n3fit/src/n3fit/scripts/n3fit_exec.py index a8c4265fe8..cf6b276037 100755 --- a/n3fit/src/n3fit/scripts/n3fit_exec.py +++ b/n3fit/src/n3fit/scripts/n3fit_exec.py @@ -3,28 +3,22 @@ n3fit - performs fit using ml external frameworks """ -import sys +import argparse +import logging +import pathlib import re import shutil -import pathlib -import logging +import sys import warnings -import argparse -from validphys.app import App -from validphys.config import Environment, Config -from validphys.config import EnvironmentError_, ConfigError -from validphys.core import FitSpec from reportengine import colors from reportengine.compat import yaml from reportengine.namespaces import NSList +from validphys.app import App +from validphys.config import Config, ConfigError, Environment, EnvironmentError_ +from validphys.core import FitSpec - -N3FIT_FIXED_CONFIG = dict( - use_cuts = 'internal', - use_t0 = True, - actions_ = [] -) +N3FIT_FIXED_CONFIG = dict(use_cuts='internal', use_t0=True, actions_=[]) FIT_NAMESPACE = "datacuts::theory::fitting " CLOSURE_NAMESPACE = "datacuts::theory::closuretest::fitting " @@ -45,6 +39,7 @@ INPUT_FOLDER = "input" TAB_FOLDER = "tables" + class N3FitError(Exception): """Exception raised when n3fit cannot succeed and knows why""" @@ -64,12 +59,12 @@ def init_output(self): # check if results folder exists self.output_path = pathlib.Path(self.output_path).absolute() - if not (self.output_path/"nnfit").is_dir(): + if not (self.output_path / "nnfit").is_dir(): if not re.fullmatch(r"[\w.\-]+", self.output_path.name): raise N3FitError("Invalid output folder name. Must be alphanumeric.") try: self.output_path.mkdir(exist_ok=True) - (self.output_path /"nnfit").mkdir(exist_ok=True) + (self.output_path / "nnfit").mkdir(exist_ok=True) except OSError as e: raise EnvironmentError_(e) from e @@ -124,7 +119,9 @@ def from_yaml(cls, o, *args, **kwargs): except yaml.error.YAMLError as e: raise ConfigError(f"Failed to parse yaml file: {e}") if not isinstance(file_content, dict): - raise ConfigError(f"Expecting input runcard to be a mapping, " f"not '{type(file_content)}'.") + raise ConfigError( + f"Expecting input runcard to be a mapping, not '{type(file_content)}'." + ) if file_content.get('closuretest') is not None: namespace = CLOSURE_NAMESPACE @@ -147,21 +144,39 @@ def from_yaml(cls, o, *args, **kwargs): validation_action = namespace + "validation_pseudodata" N3FIT_FIXED_CONFIG['actions_'].extend((training_action, validation_action)) - #Theorycovmat flags and defaults + + if thconfig := file_content.get('fiatlux'): + N3FIT_FIXED_CONFIG['fiatlux'] = thconfig + else: + N3FIT_FIXED_CONFIG['fiatlux'] = None + # Theorycovmat flags and defaults N3FIT_FIXED_CONFIG['theory_covmat_flag'] = False N3FIT_FIXED_CONFIG['use_thcovmat_in_fitting'] = False N3FIT_FIXED_CONFIG['use_thcovmat_in_sampling'] = False - if (thconfig:=file_content.get('theorycovmatconfig')) is not None: - N3FIT_FIXED_CONFIG['use_thcovmat_in_fitting'] = thconfig.get('use_thcovmat_in_fitting', True) - N3FIT_FIXED_CONFIG['use_thcovmat_in_sampling'] = thconfig.get('use_thcovmat_in_sampling', True) - if N3FIT_FIXED_CONFIG['use_thcovmat_in_sampling'] or N3FIT_FIXED_CONFIG['use_thcovmat_in_fitting']: + if (thconfig := file_content.get('theorycovmatconfig')) is not None: + N3FIT_FIXED_CONFIG['use_thcovmat_in_fitting'] = thconfig.get( + 'use_thcovmat_in_fitting', True + ) + N3FIT_FIXED_CONFIG['use_thcovmat_in_sampling'] = thconfig.get( + 'use_thcovmat_in_sampling', True + ) + if ( + N3FIT_FIXED_CONFIG['use_thcovmat_in_sampling'] + or N3FIT_FIXED_CONFIG['use_thcovmat_in_fitting'] + ): N3FIT_FIXED_CONFIG['theory_covmat_flag'] = True - N3FIT_FIXED_CONFIG['use_user_uncertainties'] = thconfig.get('use_user_uncertainties', False) - N3FIT_FIXED_CONFIG['use_scalevar_uncertainties'] = thconfig.get('use_scalevar_uncertainties', True) - #Sampling flags - if (sam_t0:=file_content.get('sampling')) is not None: - N3FIT_FIXED_CONFIG['separate_multiplicative'] = sam_t0.get('separate_multiplicative', True) - #Fitting flag + N3FIT_FIXED_CONFIG['use_user_uncertainties'] = thconfig.get( + 'use_user_uncertainties', False + ) + N3FIT_FIXED_CONFIG['use_scalevar_uncertainties'] = thconfig.get( + 'use_scalevar_uncertainties', True + ) + # Sampling flags + if (sam_t0 := file_content.get('sampling')) is not None: + N3FIT_FIXED_CONFIG['separate_multiplicative'] = sam_t0.get( + 'separate_multiplicative', True + ) + # Fitting flag file_content.update(N3FIT_FIXED_CONFIG) return cls(file_content, *args, **kwargs) @@ -178,12 +193,13 @@ def parse_fakedata(self, fakedata: bool): """ if fakedata: log.warning("using filtered closure data") - if not (self.environment.output_path/'filter').is_dir(): + if not (self.environment.output_path / 'filter').is_dir(): raise ConfigError( "Could not find filter result at " f"{self.environment.output_path/'filter'} " "to load commondata from. Did you run filter on the " - "runcard?") + "runcard?" + ) return fakedata def produce_use_fitcommondata(self, fakedata): @@ -193,8 +209,7 @@ def produce_use_fitcommondata(self, fakedata): return fakedata def produce_kfold_parameters(self, kfold=None, hyperopt=None): - """Return None even if there are kfolds in the runcard if the hyperopt flag is not active - """ + """Return None even if there are kfolds in the runcard if the hyperopt flag is not active""" if hyperopt is not None: return kfold return None @@ -232,7 +247,9 @@ def __init__(self): @property def argparser(self): parser = super().argparser - parser.add_argument("-o", "--output", help="Output folder and name of the fit", default=None) + parser.add_argument( + "-o", "--output", help="Output folder and name of the fit", default=None + ) def check_positive(value): ivalue = int(value) @@ -243,7 +260,10 @@ def check_positive(value): parser.add_argument("--hyperopt", help="Enable hyperopt scan", default=None, type=int) parser.add_argument("replica", help="MC replica number", type=check_positive) parser.add_argument( - "-r", "--replica_range", help="End of the range of replicas to compute", type=check_positive + "-r", + "--replica_range", + help="End of the range of replicas to compute", + type=check_positive, ) return parser diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py index 19722802af..c2ce728acf 100644 --- a/n3fit/src/n3fit/scripts/vp_setupfit.py +++ b/n3fit/src/n3fit/scripts/vp_setupfit.py @@ -25,39 +25,39 @@ # top. -import sys +import hashlib +import logging +import pathlib import re import shutil -import pathlib -import logging -import hashlib +import sys import warnings -from validphys.config import Environment, Config, EnvironmentError_, ConfigError -from validphys.app import App -from reportengine.compat import yaml from reportengine import colors - +from reportengine.compat import yaml +from validphys.app import App +from validphys.config import Config, ConfigError, Environment, EnvironmentError_ SETUPFIT_FIXED_CONFIG = dict( actions_=[ 'datacuts check_t0pdfset', 'theory check_positivity', - ]) + ] +) -SETUPFIT_PROVIDERS = ['validphys.filters', - 'validphys.theorycovariance.construction', - 'validphys.results', - 'validphys.covmats', - 'n3fit.n3fit_checks_provider' +SETUPFIT_PROVIDERS = [ + 'validphys.filters', + 'validphys.theorycovariance.construction', + 'validphys.results', + 'validphys.covmats', + 'n3fit.n3fit_checks_provider', ] SETUPFIT_DEFAULTS = dict( - use_cuts = 'internal', + use_cuts='internal', ) - log = logging.getLogger(__name__) RUNCARD_COPY_FILENAME = "filter.yml" @@ -69,11 +69,13 @@ class SetupFitError(Exception): """Exception raised when setup-fit cannot succeed and knows why""" + pass class SetupFitEnvironment(Environment): """Container for information to be filled at run time""" + def init_output(self): # check file exists, is a file, has extension. if not self.config_yml.exists(): @@ -122,8 +124,7 @@ def save_md5(self): @classmethod def ns_dump_description(cls): - return {'filter_path': "The filter output folder", - **super().ns_dump_description()} + return {'filter_path': "The filter output folder", **super().ns_dump_description()} class SetupFitConfig(Config): @@ -133,19 +134,19 @@ class SetupFitConfig(Config): def from_yaml(cls, o, *args, **kwargs): try: with warnings.catch_warnings(): - warnings.simplefilter('ignore', - yaml.error.MantissaNoDotYAML1_1Warning) - #We need to specify the older version 1.1 to support the - #older configuration files, which liked to use on/off for - #booleans. - #The floating point parsing yields warnings everywhere, which - #we suppress. + warnings.simplefilter('ignore', yaml.error.MantissaNoDotYAML1_1Warning) + # We need to specify the older version 1.1 to support the + # older configuration files, which liked to use on/off for + # booleans. + # The floating point parsing yields warnings everywhere, which + # we suppress. file_content = yaml.safe_load(o, version='1.1') except yaml.error.YAMLError as e: raise ConfigError(f"Failed to parse yaml file: {e}") if not isinstance(file_content, dict): - raise ConfigError(f"Expecting input runcard to be a mapping, " - f"not '{type(file_content)}'.") + raise ConfigError( + f"Expecting input runcard to be a mapping, " f"not '{type(file_content)}'." + ) if file_content.get('closuretest') is not None: filter_action = 'datacuts::closuretest::theory::fitting filter' @@ -156,8 +157,13 @@ def from_yaml(cls, o, *args, **kwargs): SETUPFIT_FIXED_CONFIG['actions_'] += [check_n3fit_action, filter_action] if file_content.get('theorycovmatconfig') is not None: SETUPFIT_FIXED_CONFIG['actions_'].append( - 'datacuts::theory::theorycovmatconfig nnfit_theory_covmat') - for k,v in SETUPFIT_DEFAULTS.items(): + 'datacuts::theory::theorycovmatconfig nnfit_theory_covmat' + ) + if file_content.get('fiatlux') is not None: + SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset') + if file_content.get('fiatlux')["additional_errors"]: + SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors') + for k, v in SETUPFIT_DEFAULTS.items(): file_content.setdefault(k, v) file_content.update(SETUPFIT_FIXED_CONFIG) return cls(file_content, *args, **kwargs) @@ -165,19 +171,19 @@ def from_yaml(cls, o, *args, **kwargs): class SetupFitApp(App): """The class which parsers and perform the filtering""" + environment_class = SetupFitEnvironment config_class = SetupFitConfig def __init__(self): - super(SetupFitApp, self).__init__(name='setup-fit', - providers=SETUPFIT_PROVIDERS) + super(SetupFitApp, self).__init__(name='setup-fit', providers=SETUPFIT_PROVIDERS) @property def argparser(self): parser = super().argparser - parser.add_argument('-o','--output', - help="Output folder and name of the fit", - default=None) + parser.add_argument( + '-o', '--output', help="Output folder and name of the fit", default=None + ) return parser def get_commandline_arguments(self, cmdline=None): @@ -201,9 +207,7 @@ def run(self): sys.exit(1) except Exception as e: log.critical(f"Bug in setup-fit ocurred. Please report it.") - print( - colors.color_exception(e.__class__, e, e.__traceback__), - file=sys.stderr) + print(colors.color_exception(e.__class__, e, e.__traceback__), file=sys.stderr) sys.exit(1) diff --git a/n3fit/src/n3fit/tests/test_layers.py b/n3fit/src/n3fit/tests/test_layers.py index b1d662450f..b8cbb844a0 100644 --- a/n3fit/src/n3fit/tests/test_layers.py +++ b/n3fit/src/n3fit/tests/test_layers.py @@ -259,3 +259,22 @@ def test_mask(): masker = layers.Mask(bool_mask=np_mask, c=rn_val) ret = masker(fi) np.testing.assert_allclose(ret, masked_fi * rn_val, rtol=1e-5) + +def test_addphoton_init(): + """Test AddPhoton class.""" + addphoton = layers.AddPhoton(photons=None) + np.testing.assert_equal(addphoton._photons_generator, None) + addphoton = layers.AddPhoton(photons=1234) + np.testing.assert_equal(addphoton._photons_generator, 1234) + np.testing.assert_equal(addphoton._pdf_ph, None) + +class FakePhoton(): + def __call__(self, xgrid): + return [np.exp(-xgrid)] + +def test_compute_photon(): + photon = FakePhoton() + addphoton = layers.AddPhoton(photons=photon) + xgrid = np.geomspace(1e-4, 1., 10) + addphoton.register_photon(xgrid) + np.testing.assert_allclose(addphoton._pdf_ph, [np.exp(-xgrid)]) diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py index 777cac91ac..52cf001980 100644 --- a/n3fit/src/n3fit/vpinterface.py +++ b/n3fit/src/n3fit/vpinterface.py @@ -18,14 +18,16 @@ """ -import logging from collections.abc import Iterable +import logging + import numpy as np import numpy.linalg as la + +from validphys.arclength import arc_lengths, integrability_number from validphys.core import PDF, MCStats -from validphys.pdfbases import ALL_FLAVOURS, check_basis from validphys.lhapdfset import LHAPDFSet -from validphys.arclength import integrability_number, arc_lengths +from validphys.pdfbases import ALL_FLAVOURS, check_basis log = logging.getLogger(__name__) # Order of the evolution basis output from n3fit @@ -78,6 +80,18 @@ def xfxQ(self, x, Q, n, fl): ) return self.grid_values([fl], [x]).squeeze()[n] + def _register_photon(self, xgrid): + """If the PDF models contain photons, register the xgrid with them""" + for m in self._lhapdf_set: + pl = m.get_layer_re("add_photon") + # if pl is an empy list there's no photon + if not pl: + continue + pl[0].register_photon(xgrid) + # Recompile the model if necessary + if not pl[0].built: + m.compile() + def __call__(self, xarr, flavours=None, replica=None): """Uses the internal model to produce pdf values for the grid The output is on the evolution basis. @@ -103,9 +117,14 @@ def __call__(self, xarr, flavours=None, replica=None): # as the scaling is done by the model itself mod_xgrid = xarr.reshape(1, -1, 1) + # Register the grid with the photon + self._register_photon(mod_xgrid) + if replica is None or replica == 0: # We need generate output values for all replicas - result = np.concatenate([m.predict({"pdf_input": mod_xgrid}) for m in self._lhapdf_set], axis=0) + result = np.concatenate( + [m.predict({"pdf_input": mod_xgrid}) for m in self._lhapdf_set], axis=0 + ) if replica == 0: # We want _only_ the central value result = np.mean(result, axis=0, keepdims=True) diff --git a/nnpdfcpp/data/theory.db b/nnpdfcpp/data/theory.db index ef5c669b78..1f953e2459 100644 Binary files a/nnpdfcpp/data/theory.db and b/nnpdfcpp/data/theory.db differ diff --git a/validphys2/src/validphys/config.py b/validphys2/src/validphys/config.py index f82781536c..a9ee17a651 100644 --- a/validphys2/src/validphys/config.py +++ b/validphys2/src/validphys/config.py @@ -4,58 +4,47 @@ @author: Zahari Kassabov """ -import logging -import pathlib +from collections import ChainMap, defaultdict +from collections.abc import Mapping, Sequence +import copy import functools +import glob +from importlib.resources import contents, read_text import inspect +import logging import numbers -import copy -from importlib.resources import read_text, contents - -from collections import ChainMap, defaultdict -from collections.abc import Mapping, Sequence +import pathlib import pandas as pd -import glob -from reportengine import configparser +from reportengine import configparser, report +from reportengine.compat import yaml +from reportengine.configparser import ConfigError, _parse_func, element_of, record_from_defaults from reportengine.environment import Environment, EnvironmentError_ -from reportengine.configparser import ( - ConfigError, - element_of, - _parse_func, - record_from_defaults, -) from reportengine.helputils import get_parser_type from reportengine.namespaces import NSList -from reportengine import report -from reportengine.compat import yaml - from validphys.core import ( + CutsPolicy, DataGroupSpec, DataSetInput, ExperimentInput, - CutsPolicy, MatchedCuts, SimilarCuts, ThCovMatSpec, ) from validphys.fitdata import fitted_replica_indexes, num_fitted_replicas +from validphys.gridvalues import LUMI_CHANNELS from validphys.loader import ( + DataNotFoundError, + FallbackLoader, + InconsistentMetaDataError, Loader, LoaderError, LoadFailedError, - DataNotFoundError, PDFNotFound, - FallbackLoader, - InconsistentMetaDataError, ) -from validphys.gridvalues import LUMI_CHANNELS - from validphys.paramfits.config import ParamfitsConfig - from validphys.plotoptions import get_info - import validphys.scalevariations log = logging.getLogger(__name__) @@ -65,7 +54,13 @@ class Environment(Environment): """Container for information to be filled at run time""" def __init__( - self, *, this_folder=None, net=True, upload=False, dry=False, **kwargs, + self, + *, + this_folder=None, + net=True, + upload=False, + dry=False, + **kwargs, ): if this_folder: self.this_folder = pathlib.Path(this_folder) @@ -84,10 +79,7 @@ def __init__( try: self.loader = loader_class() except LoaderError as e: - log.error( - "Failed to find the paths. These are configured " - "in the nnprofile settings" - ) + log.error("Failed to find the paths. These are configured " "in the nnprofile settings") raise EnvironmentError_(e) from e self.deta_path = self.loader.datapath self.results_path = self.loader.resultspath @@ -207,9 +199,7 @@ def parse_use_cuts(self, use_cuts: (bool, str)): elif isinstance(use_cuts, str) and use_cuts in valid_cuts: res = CutsPolicy(use_cuts) else: - raise ConfigError( - f"Invalid use_cuts setting: '{use_cuts}'.", use_cuts, valid_cuts - ) + raise ConfigError(f"Invalid use_cuts setting: '{use_cuts}'.", use_cuts, valid_cuts) return res @@ -223,7 +213,7 @@ def produce_inclusive_use_scalevar_uncertainties( point_prescription: (str, None) = None, ): """Whether to use a scale variation uncertainty theory covmat. - Checks whether a point prescription is included in the runcard and if so + Checks whether a point prescription is included in the runcard and if so assumes scale uncertainties are to be used.""" if (not use_scalevar_uncertainties) and (point_prescription is not None): use_scalevar_uncertainties = True @@ -255,8 +245,7 @@ def produce_pdfreplicas(self, fitpdf): return NSList(replicas, nskey="replica") def produce_fitcontextwithcuts(self, fit, fitinputcontext): - """Like fitinputcontext but setting the cuts policy. - """ + """Like fitinputcontext but setting the cuts policy.""" theoryid = fitinputcontext["theoryid"] data_input = fitinputcontext["data_input"] @@ -331,18 +320,13 @@ def parse_hyperscan(self, hyperscan): try: return self.loader.check_hyperscan(hyperscan) except LoadFailedError as e: - raise ConfigError( - str(e), hyperscan, self.loader.available_hyperscans - ) from e + raise ConfigError(str(e), hyperscan, self.loader.available_hyperscans) from e def parse_hyperscan_config(self, hyperscan_config, hyperopt=None): - """Configuration of the hyperscan - """ + """Configuration of the hyperscan""" if "from_hyperscan" in hyperscan_config: hyperscan = self.parse_hyperscan(hyperscan_config["from_hyperscan"]) - log.info( - "Using previous hyperscan: '%s' to generate the search space", hyperscan - ) + log.info("Using previous hyperscan: '%s' to generate the search space", hyperscan) return hyperscan.as_input().get("hyperscan_config") if "use_tries_from" in hyperscan_config: @@ -378,8 +362,7 @@ def produce_multiclosure_underlyinglaw(self, fits): if len(laws) != 1: raise ConfigError( - "Did not find unique underlying law from fits, " - f"instead found: {laws}" + "Did not find unique underlying law from fits, " f"instead found: {laws}" ) return self.parse_pdf(laws.pop()) @@ -398,7 +381,7 @@ def produce_basisfromfit(self, fit): return {"basis": basis} def produce_fitpdfandbasis(self, fitpdf, basisfromfit): - """ Set the PDF and basis from the fit config. """ + """Set the PDF and basis from the fit config.""" return {**fitpdf, **basisfromfit} @element_of("dataset_inputs") @@ -417,9 +400,7 @@ def parse_dataset_input(self, dataset: Mapping): if name.startswith("POS"): raise ConfigError("Please, use `posdataset` for positivity") except KeyError: - raise ConfigError( - "'dataset' must be a mapping with " "'dataset' and 'sysnum'" - ) + raise ConfigError("'dataset' must be a mapping with " "'dataset' and 'sysnum'") sysnum = dataset.get("sys") cfac = dataset.get("cfac", tuple()) @@ -437,9 +418,7 @@ def parse_dataset_input(self, dataset: Mapping): kdiff = dataset.keys() - known_keys for k in kdiff: # Abuse ConfigError to get the suggestions. - log.warning( - ConfigError(f"Key '{k}' in dataset_input not known.", k, known_keys) - ) + log.warning(ConfigError(f"Key '{k}' in dataset_input not known.", k, known_keys)) return DataSetInput( name=name, sys=sysnum, @@ -502,15 +481,11 @@ def _produce_matched_cuts(self, commondata): self._check_dataspecs_type(nss) if not nss: - raise ConfigError( - "'cuts_intersection_spec' must contain at least one namespace." - ) + raise ConfigError("'cuts_intersection_spec' must contain at least one namespace.") for ns in nss: with self.set_context( - ns=self._curr_ns.new_child(ns).new_child( - {"use_cuts": CutsPolicy.INTERNAL} - ) + ns=self._curr_ns.new_child(ns).new_child({"use_cuts": CutsPolicy.INTERNAL}) ): # Note: Do not call _produce_internal_cuts directly here: # That doesn't correctly set the namespace in a way that `rules` @@ -520,7 +495,7 @@ def _produce_matched_cuts(self, commondata): return MatchedCuts(cut_list, ndata=ndata) def _produce_similarity_cuts(self, commondata): - """ Compute the intersection between two namespaces (similar to + """Compute the intersection between two namespaces (similar to `fromintersection`) but additionally require that the predictions computed for each dataset across the namespaces are *similar*, specifically that the ratio between the absolute difference in the @@ -544,22 +519,24 @@ def _produce_similarity_cuts(self, commondata): ) try: - _, exclusion_list = self.parse_from_( - None, "do_not_require_similarity_for", write=False - ) + _, exclusion_list = self.parse_from_(None, "do_not_require_similarity_for", write=False) except configparser.InputNotFoundError: exclusion_list = [] name = commondata.name # slightly circular here, since matched cuts will re-produce nss if name in exclusion_list: - with self.set_context( - ns=self._curr_ns.new_child({"use_cuts": CutsPolicy.INTERNAL}) - ): + with self.set_context(ns=self._curr_ns.new_child({"use_cuts": CutsPolicy.INTERNAL})): return self.parse_from_(None, "cuts", write=False)[1] matched_cuts = self._produce_matched_cuts(commondata) inps = [] for i, ns in enumerate(nss): - with self.set_context(ns=self._curr_ns.new_child({**ns,})): + with self.set_context( + ns=self._curr_ns.new_child( + { + **ns, + } + ) + ): # TODO: find a way to not duplicate this and use a dict # instead of a linear search _, dins = self.parse_from_(None, "dataset_inputs", write=False) @@ -613,9 +590,9 @@ def produce_dataset( check_plotting: bool = False, ): """Dataset specification from the theory and CommonData. - Use the cuts from the fit, if provided. If check_plotting is set to - True, attempt to lod and check the PLOTTING files - (note this may cause a noticeable slowdown in general).""" + Use the cuts from the fit, if provided. If check_plotting is set to + True, attempt to lod and check the PLOTTING files + (note this may cause a noticeable slowdown in general).""" name = dataset_input.name sysnum = dataset_input.sys cfac = dataset_input.cfac @@ -652,8 +629,8 @@ def produce_dataset( @configparser.element_of("experiments") def parse_experiment(self, experiment: dict): """A set of datasets where correlated systematics are taken - into account. It is a mapping where the keys are the experiment - name 'experiment' and a list of datasets.""" + into account. It is a mapping where the keys are the experiment + name 'experiment' and a list of datasets.""" try: name, datasets = experiment["experiment"], experiment["datasets"] except KeyError as e: @@ -684,9 +661,7 @@ def parse_experiment_input(self, ei: dict): return ExperimentInput(name=name, datasets=datasets) # TODO: Do away with the mapping and make the conversion implicitly - def produce_experiment_from_input( - self, experiment_input, theoryid, use_cuts, fit=None - ): + def produce_experiment_from_input(self, experiment_input, theoryid, use_cuts, fit=None): """Return a mapping containing a single experiment from an experiment input. NOTE: This might be deprecated in the future.""" return { @@ -700,8 +675,8 @@ def produce_experiment_from_input( def produce_sep_mult(self, separate_multiplicative=None): """ - Specifies whether to separate the multiplicative errors in the - experimental covmat construction. The default is True. + Specifies whether to separate the multiplicative errors in the + experimental covmat construction. The default is True. """ if separate_multiplicative is False: return False @@ -761,7 +736,7 @@ def produce_loaded_theory_covmat( ): """ Loads the theory covmat from the correct file according to how it - was generated by vp-setupfit. + was generated by vp-setupfit. """ if theory_covmat_flag is False: return 0.0 @@ -769,34 +744,21 @@ def produce_loaded_theory_covmat( generic_path = "datacuts_theory_theorycovmatconfig_theory_covmat_custom.csv" if use_user_uncertainties is True: if use_scalevar_uncertainties is True: - generic_path = ( - "datacuts_theory_theorycovmatconfig_total_theory_covmat.csv" - ) + generic_path = "datacuts_theory_theorycovmatconfig_total_theory_covmat.csv" else: generic_path = "datacuts_theory_theorycovmatconfig_user_covmat.csv" # check if there are multiple files files = glob.glob(str(output_path / "tables/*theorycovmat*")) paths = [ - str( - output_path - / "tables/datacuts_theory_theorycovmatconfig_theory_covmat_custom.csv" - ), - str( - output_path - / "tables/datacuts_theory_theorycovmatconfig_total_theory_covmat.csv" - ), - str( - output_path - / "tables/datacuts_theory_theorycovmatconfig_user_covmat.csv" - ), + str(output_path / "tables/datacuts_theory_theorycovmatconfig_theory_covmat_custom.csv"), + str(output_path / "tables/datacuts_theory_theorycovmatconfig_total_theory_covmat.csv"), + str(output_path / "tables/datacuts_theory_theorycovmatconfig_user_covmat.csv"), ] paths.remove(str(output_path / "tables" / generic_path)) for f in files: for path in paths: if f == path: - raise ValueError( - "More than one theory_covmat file in folder tables" - ) + raise ValueError("More than one theory_covmat file in folder tables") theorypath = output_path / "tables" / generic_path theory_covmat = pd.read_csv( theorypath, @@ -863,8 +825,7 @@ def produce_dataset_inputs_covariance_matrix(self, use_pdferr: bool = False): def _check_dataspecs_type(dataspecs): if not isinstance(dataspecs, Sequence): raise ConfigError( - "dataspecs should be a sequence of mappings, not " - f"{type(dataspecs).__name__}" + "dataspecs should be a sequence of mappings, not " f"{type(dataspecs).__name__}" ) for spec in dataspecs: @@ -918,7 +879,12 @@ def produce_matched_datasets_from_dataspecs(self, dataspecs): inner_spec_list = inres["dataspecs"] = [] for ispec, spec in enumerate(dataspecs): # Passing spec by referene - d = ChainMap({"dataset_input": all_names[ispec][k],}, spec) + d = ChainMap( + { + "dataset_input": all_names[ispec][k], + }, + spec, + ) inner_spec_list.append(d) res.append(inres) res.sort(key=lambda x: (x["process"], x["dataset_name"])) @@ -942,7 +908,12 @@ def produce_matched_positivity_from_dataspecs(self, dataspecs): l = inres["dataspecs"] = [] for ispec, spec in enumerate(dataspecs): # Passing spec by referene - d = ChainMap({"posdataset": all_names[ispec][k],}, spec) + d = ChainMap( + { + "posdataset": all_names[ispec][k], + }, + spec, + ) l.append(d) res.append(inres) res.sort(key=lambda x: (x["posdataset_name"])) @@ -1011,13 +982,9 @@ def produce_combined_shift_and_theory_dataspecs(self, theoryconfig, shiftconfig) # TODO: Worth it to do some black magic to not pass params explicitly? # Note that `parse_experiments` doesn't exist yet. - def parse_reweighting_experiments( - self, experiments, *, theoryid, use_cuts, fit=None - ): + def parse_reweighting_experiments(self, experiments, *, theoryid, use_cuts, fit=None): """A list of experiments to be used for reweighting.""" - return self.parse_experiments( - experiments, theoryid=theoryid, use_cuts=use_cuts, fit=fit - ) + return self.parse_experiments(experiments, theoryid=theoryid, use_cuts=use_cuts, fit=fit) def parse_t0pdfset(self, name): """PDF set used to generate the t0 covmat.""" @@ -1029,7 +996,9 @@ def parse_use_t0(self, do_use_t0: bool): # TODO: Find a good name for this def produce_t0set( - self, t0pdfset=None, use_t0=False, + self, + t0pdfset=None, + use_t0=False, ): """Return the t0set if use_t0 is True and None otherwise. Raises an error if t0 is requested but no t0set is given. @@ -1040,18 +1009,36 @@ def produce_t0set( return t0pdfset return None + def parse_luxset(self, name): + """PDF set used to generate the photon with fiatlux.""" + return self.parse_pdf(name) + + def parse_additional_errors(self, bool): + """PDF set used to generate the photon additional errors: + they are constructed using the replicas 101-107 of the PDF set + LUXqed17_plus_PDF4LHC15_nnlo_100 (that are obtained varying some + parameters of the LuxQED approach) in the way described + in sec. 2.5 of https://arxiv.org/pdf/1712.07053.pdf + """ + if bool: + return self.parse_pdf("LUXqed17_plus_PDF4LHC15_nnlo_100") + else: + return False + def parse_fakepdf(self, name): """PDF set used to generate the fake data in a closure test.""" return self.parse_pdf(name) def _parse_lagrange_multiplier(self, kind, theoryid, setdict): - """ Lagrange multiplier constraints are mappings + """Lagrange multiplier constraints are mappings containing a `dataset` and a `maxlambda` argument which - defines the maximum value allowed for the multiplier """ - bad_msg = f"{kind} must be a mapping with a name ('dataset') and a float multiplier (maxlambda)" + defines the maximum value allowed for the multiplier""" + bad_msg = ( + f"{kind} must be a mapping with a name ('dataset') and a float multiplier (maxlambda)" + ) theoryno, _ = theoryid lambda_key = "maxlambda" - #BCH allow for old-style runcards with 'poslambda' instead of 'maxlambda' + # BCH allow for old-style runcards with 'poslambda' instead of 'maxlambda' if "poslambda" in setdict and "maxlambda" not in setdict: log.warning("The `poslambda` argument has been deprecated in favour of `maxlambda`") lambda_key = "poslambda" @@ -1078,8 +1065,7 @@ def parse_posdataset(self, posset: dict, *, theoryid): def produce_posdatasets(self, positivity): if not isinstance(positivity, dict) or "posdatasets" not in positivity: raise ConfigError( - "Failed to get 'posdatasets' from positivity. " - "Expected that key to be present." + "Failed to get 'posdatasets' from positivity. " "Expected that key to be present." ) return positivity["posdatasets"] @@ -1102,12 +1088,8 @@ def produce_reweight_all_datasets(self, experiments): ret = [] for experiment in experiments: for dsinput, dataset in zip(experiment, experiment.datasets): - single_exp = DataGroupSpec( - experiment.name, datasets=[dataset], dsinputs=[dsinput] - ) - ret.append( - {"reweighting_experiments": [single_exp], "dataset_input": dsinput} - ) + single_exp = DataGroupSpec(experiment.name, datasets=[dataset], dsinputs=[dsinput]) + ret.append({"reweighting_experiments": [single_exp], "dataset_input": dsinput}) return ret """ @@ -1169,16 +1151,12 @@ def produce_nnfit_theory_covmat( if inclusive_use_scalevar_uncertainties: if use_user_uncertainties: # Both scalevar and user uncertainties - from validphys.theorycovariance.construction import ( - total_theory_covmat_fitting, - ) + from validphys.theorycovariance.construction import total_theory_covmat_fitting f = total_theory_covmat_fitting else: # Only scalevar uncertainties - from validphys.theorycovariance.construction import ( - theory_covmat_custom_fitting, - ) + from validphys.theorycovariance.construction import theory_covmat_custom_fitting f = theory_covmat_custom_fitting elif use_user_uncertainties: @@ -1203,20 +1181,14 @@ def produce_fitthcovmat( theory covariance matrix then returns `False`. """ if not isinstance(use_thcovmat_if_present, bool): - raise ConfigError( - "use_thcovmat_if_present should be a boolean, by default it is False" - ) + raise ConfigError("use_thcovmat_if_present should be a boolean, by default it is False") if use_thcovmat_if_present and not fit: - raise ConfigError( - "`use_thcovmat_if_present` was true but no `fit` was specified." - ) + raise ConfigError("`use_thcovmat_if_present` was true but no `fit` was specified.") if use_thcovmat_if_present and fit: try: - thcovmat_present = fit.as_input()["theorycovmatconfig"][ - "use_thcovmat_in_fitting" - ] + thcovmat_present = fit.as_input()["theorycovmatconfig"]["use_thcovmat_in_fitting"] except KeyError: # assume covmat wasn't used and fill in key accordingly but warn user log.warning( @@ -1264,21 +1236,17 @@ def parse_fitdeclaration(self, label: str): return label def produce_all_commondata(self): - """produces all commondata using the loader function """ + """produces all commondata using the loader function""" ds_names = self.loader.available_datasets ds_inputs = [self.parse_dataset_input({"dataset": ds}) for ds in ds_names] - cd_out = [ - self.produce_commondata(dataset_input=ds_input) for ds_input in ds_inputs - ] + cd_out = [self.produce_commondata(dataset_input=ds_input) for ds_input in ds_inputs] return cd_out def parse_groupby(self, grouping: str): """parses the groupby key and checks it is an allowed grouping""" # TODO: think if better way to do this properly if grouping not in ["experiment", "nnpdf31_process"]: - raise ConfigError( - f"Grouping not available: {grouping}, did you spell it " "correctly?" - ) + raise ConfigError(f"Grouping not available: {grouping}, did you spell it " "correctly?") return grouping def parse_norm_threshold(self, val: (numbers.Number, type(None))): @@ -1312,9 +1280,7 @@ def load_default_default_filter_rules(self, spec): lock_token = "_filters.lock.yaml" try: - return yaml.safe_load( - read_text(validphys.cuts.lockfiles, f"{spec}{lock_token}") - ) + return yaml.safe_load(read_text(validphys.cuts.lockfiles, f"{spec}{lock_token}")) except FileNotFoundError as e: alternatives = [ el.strip(lock_token) @@ -1352,11 +1318,7 @@ def produce_rules( ): """Produce filter rules based on the user defined input and defaults.""" - from validphys.filters import ( - Rule, - RuleProcessingError, - default_filter_rules_input, - ) + from validphys.filters import Rule, RuleProcessingError, default_filter_rules_input theory_parameters = theoryid.get_description() @@ -1393,9 +1355,7 @@ def load_default_default_filter_settings(self, spec): lock_token = "_defaults.lock.yaml" try: - return yaml.safe_load( - read_text(validphys.cuts.lockfiles, f"{spec}{lock_token}") - ) + return yaml.safe_load(read_text(validphys.cuts.lockfiles, f"{spec}{lock_token}")) except FileNotFoundError as e: alternatives = alternatives = [ el.strip(lock_token) @@ -1431,23 +1391,13 @@ def produce_defaults( """ from validphys.filters import default_filter_settings_input - if ( - q2min is not None - and "q2min" in filter_defaults - and q2min != filter_defaults["q2min"] - ): + if q2min is not None and "q2min" in filter_defaults and q2min != filter_defaults["q2min"]: raise ConfigError("q2min defined multiple times with different values") - if ( - w2min is not None - and "w2min" in filter_defaults - and w2min != filter_defaults["w2min"] - ): + if w2min is not None and "w2min" in filter_defaults and w2min != filter_defaults["w2min"]: raise ConfigError("w2min defined multiple times with different values") if default_filter_settings_recorded_spec_ is not None: - filter_defaults = default_filter_settings_recorded_spec_[ - default_filter_settings - ] + filter_defaults = default_filter_settings_recorded_spec_[default_filter_settings] # If we find recorded specs return immediately and don't read q2min and w2min # from runcard return filter_defaults @@ -1468,7 +1418,10 @@ def produce_defaults( return filter_defaults def produce_data( - self, data_input, *, group_name="data", + self, + data_input, + *, + group_name="data", ): """A set of datasets where correlated systematics are taken into account @@ -1525,15 +1478,9 @@ def _parse_data_input_from_( ) # We need to make theoryid available if using experiments try: - _, experiments = self.parse_from_( - parse_from_value, data_key, write=False - ) + _, experiments = self.parse_from_(parse_from_value, data_key, write=False) data_val = NSList( - [ - dsinput - for experiment in experiments - for dsinput in experiment.dsinputs - ], + [dsinput for experiment in experiments for dsinput in experiment.dsinputs], nskey="dataset_input", ) except ConfigError as inner_error: @@ -1570,22 +1517,23 @@ def load_default_data_grouping(self, spec): """Load the default grouping of data""" # slightly superfluous, only one default at present but perhaps # somebody will want to add to this at some point e.g for th. uncertainties - allowed = { - "standard_report": "experiment", - "thcovmat_fit": "ALL" - } + allowed = {"standard_report": "experiment", "thcovmat_fit": "ALL"} return allowed[spec] def produce_processed_data_grouping( - self, use_thcovmat_in_fitting=False, use_thcovmat_in_sampling=False, data_grouping=None, data_grouping_recorded_spec_=None + self, + use_thcovmat_in_fitting=False, + use_thcovmat_in_sampling=False, + data_grouping=None, + data_grouping_recorded_spec_=None, ): """Process the data_grouping key from the runcard, or lockfile. If `data_grouping_recorded_spec_` is present then its value is taken, and the runcard is assumed to be a lockfile. If data_grouping is None, then, if either use_thcovmat_in_fitting or use_thcovmat_in_sampling - (or both) are true (which means that the fit is a thcovmat fit), group all the datasets - together, otherwise fall back to the default behaviour of grouping by + (or both) are true (which means that the fit is a thcovmat fit), group all the datasets + together, otherwise fall back to the default behaviour of grouping by experiment (called standard_report). Else, the user can specfiy their own grouping, for example metadata_process. @@ -1599,9 +1547,7 @@ def produce_processed_data_grouping( return data_grouping_recorded_spec_[data_grouping] return self.load_default_data_grouping(data_grouping) - def produce_processed_metadata_group( - self, processed_data_grouping, metadata_group=None - ): + def produce_processed_metadata_group(self, processed_data_grouping, metadata_group=None): """Expose the final data grouping result. Either metadata_group is specified by user, in which case uses `processed_data_grouping` which is experiment by default. @@ -1611,7 +1557,9 @@ def produce_processed_metadata_group( return metadata_group def produce_group_dataset_inputs_by_metadata( - self, data_input, processed_metadata_group, + self, + data_input, + processed_metadata_group, ): """Take the data and the processed_metadata_group key and attempt to group the data, returns a list where each element specifies the data_input @@ -1665,28 +1613,21 @@ def produce_group_dataset_inputs_by_experiment(self, data_input): return self.produce_group_dataset_inputs_by_metadata(data_input, "experiment") def produce_group_dataset_inputs_by_process(self, data_input): - return self.produce_group_dataset_inputs_by_metadata( - data_input, "nnpdf31_process" - ) + return self.produce_group_dataset_inputs_by_metadata(data_input, "nnpdf31_process") def produce_scale_variation_theories(self, theoryid, point_prescription): """Produces a list of theoryids given a theoryid at central scales and a point - prescription. The options for the latter are '3 point', '5 point', '5bar point', '7 point' - and '9 point'. Note that these are defined in arXiv:1906.10698. This hard codes the - theories needed for each prescription to avoid user error.""" + prescription. The options for the latter are '3 point', '5 point', '5bar point', '7 point' + and '9 point'. Note that these are defined in arXiv:1906.10698. This hard codes the + theories needed for each prescription to avoid user error.""" pp = point_prescription th = theoryid.id - lsv = yaml.safe_load( - read_text(validphys.scalevariations, "scalevariationtheoryids.yaml") - ) + lsv = yaml.safe_load(read_text(validphys.scalevariations, "scalevariationtheoryids.yaml")) scalevarsfor_list = lsv["scale_variations_for"] # Allowed central theoryids - cent_thids = [ - str(scalevarsfor_dict["theoryid"]) - for scalevarsfor_dict in scalevarsfor_list - ] + cent_thids = [str(scalevarsfor_dict["theoryid"]) for scalevarsfor_dict in scalevarsfor_list] if th not in cent_thids: valid_thids = ", ".join(cent_thids) raise ConfigError( diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py index e18b764d06..91c80ce2fb 100644 --- a/validphys2/src/validphys/filters.py +++ b/validphys2/src/validphys/filters.py @@ -2,20 +2,18 @@ Filters for NNPDF fits """ -import logging -import re from collections.abc import Mapping from importlib.resources import read_text +import logging +import re import numpy as np from reportengine.checks import check, make_check from reportengine.compat import yaml +from validphys.commondatawriter import write_commondata_to_file, write_systype_to_file import validphys.cuts -from validphys.commondatawriter import ( - write_commondata_to_file, - write_systype_to_file, - ) + log = logging.getLogger(__name__) KIN_LABEL = { @@ -53,6 +51,7 @@ class BadPerturbativeOrder(ValueError): """Exception raised when the perturbative order string is not recognized.""" + class MissingRuleAttribute(RuleProcessingError, AttributeError): """Exception raised when a rule is missing required attributes.""" @@ -75,13 +74,14 @@ def default_filter_rules_input(): return yaml.safe_load(read_text(validphys.cuts, "filters.yaml")) - def check_nonnegative(var: str): """Ensure that `var` is positive""" + @make_check def run_check(ns, **kwargs): val = ns[var] check(val >= 0, f"'{var}' must be positive or equal zero, but it is {val!r}.") + return run_check @@ -124,9 +124,7 @@ def filter_closure_data_by_experiment( res = [] for exp in experiments_data: - experiment_index = experiments_index[ - experiments_index.isin([exp.name], level=0) - ] + experiment_index = experiments_index[experiments_index.isin([exp.name], level=0)] res.append( _filter_closure_data( filter_path, exp, fakepdf, fakenoise, filterseed, experiment_index, sep_mult @@ -137,8 +135,7 @@ def filter_closure_data_by_experiment( def filter_real_data(filter_path, data): - """Filter real data, cutting any points which do not pass the filter rules. - """ + """Filter real data, cutting any points which do not pass the filter rules.""" log.info('Filtering real data.') return _filter_real_data(filter_path, data) @@ -158,8 +155,9 @@ def _write_ds_cut_data(path, dataset): log.info("All {all_ndata} points in in {dataset.name} passed kinematic cuts.") else: filtered_dsndata = len(datamask) - log.info(f"{len(datamask)}/{all_dsndata} datapoints " - f"in {dataset.name} passed kinematic cuts.") + log.info( + f"{len(datamask)}/{all_dsndata} datapoints " f"in {dataset.name} passed kinematic cuts." + ) # save to disk if datamask is not None: export_mask(path / f'FKMASK_{dataset.name}.dat', datamask) @@ -181,8 +179,7 @@ def _filter_real_data(filter_path, data): def _filter_closure_data( - filter_path, data, fakepdf, fakenoise, filterseed, experiments_index - , sep_mult + filter_path, data, fakepdf, fakenoise, filterseed, experiments_index, sep_mult ): """ This function is accessed within a closure test only, that is, the fakedata @@ -234,25 +231,19 @@ def _filter_closure_data( closure_data = level0_commondata_wc(data, fakepdf) for dataset in data.datasets: - #== print number of points passing cuts, make dataset directory and write FKMASK ==# + # == print number of points passing cuts, make dataset directory and write FKMASK ==# path = filter_path / dataset.name nfull, ncut = _write_ds_cut_data(path, dataset) make_dataset_dir(path / "systypes") total_data_points += nfull total_cut_data_points += ncut - + if fakenoise: - #======= Level 1 closure test =======# - - closure_data = make_level1_data( - data, - closure_data, - filterseed, - experiments_index, - sep_mult - ) + # ======= Level 1 closure test =======# + + closure_data = make_level1_data(data, closure_data, filterseed, experiments_index, sep_mult) - #====== write commondata and systype files ======# + # ====== write commondata and systype files ======# if fakenoise: log.info("Writing Level1 data") else: @@ -260,12 +251,7 @@ def _filter_closure_data( for cd in closure_data: path_cd = filter_path / cd.setname / f"DATA_{cd.setname}.dat" - path_sys = ( - filter_path - / cd.setname - / "systypes" - / f"SYSTYPE_{cd.setname}_DEFAULT.dat" - ) + path_sys = filter_path / cd.setname / "systypes" / f"SYSTYPE_{cd.setname}_DEFAULT.dat" write_commondata_to_file(commondata=cd, path=path_cd) write_systype_to_file(commondata=cd, path=path_sys) @@ -278,6 +264,18 @@ def check_t0pdfset(t0pdfset): log.info(f'{t0pdfset} T0 checked.') +def check_luxset(luxset): + """Lux pdf check""" + luxset.load() + log.info(f'{luxset} Lux pdf checked.') + + +def check_additional_errors(additional_errors): + """Lux additional errors pdf check""" + additional_errors.load() + log.info(f'{additional_errors} Lux additional errors pdf checked.') + + def check_positivity(posdatasets): """Verify positive datasets are ready for the fit.""" log.info('Verifying positivity tables:') @@ -285,6 +283,7 @@ def check_positivity(posdatasets): pos.load_commondata() log.info(f'{pos.name} checked.') + def check_integrability(integdatasets): """Verify positive datasets are ready for the fit.""" log.info('Verifying integrability tables:') @@ -292,6 +291,7 @@ def check_integrability(integdatasets): integ.load_commondata() log.info(f'{integ.name} checked.') + class PerturbativeOrder: """Class that conveniently handles perturbative order declarations for use @@ -362,6 +362,7 @@ def __contains__(self, i): else: return i == self.numeric_pto + class Rule: """Rule object to be used to generate cuts mask. @@ -410,9 +411,7 @@ def __init__( raise MissingRuleAttribute("No rule defined.") if self.dataset is None and self.process_type is None: - raise MissingRuleAttribute( - "Please define either a process type or dataset." - ) + raise MissingRuleAttribute("Please define either a process type or dataset.") if self.process_type is None: from validphys.loader import Loader, LoaderError @@ -422,9 +421,7 @@ def __init__( try: cd = loader.check_commondata(self.dataset) except LoaderError as e: - raise RuleProcessingError( - f"Could not find dataset {self.dataset}" - ) from e + raise RuleProcessingError(f"Could not find dataset {self.dataset}") from e if cd.process_type[:3] == "DIS": self.variables = KIN_LABEL["DIS"] else: @@ -445,9 +442,7 @@ def __init__( if hasattr(self, "PTO"): if not isinstance(self.PTO, str): - raise RuleProcessingError( - f"Expecting PTO to be a string, not {type(self.PTO)}." - ) + raise RuleProcessingError(f"Expecting PTO to be a string, not {type(self.PTO)}.") try: self.PTO = PerturbativeOrder(self.PTO) except BadPerturbativeOrder as e: @@ -482,9 +477,7 @@ def __init__( try: self.rule = compile(self.rule, "rule", "eval") except Exception as e: - raise RuleProcessingError( - f"Could not process rule {self.rule_string!r}: {e}" - ) from e + raise RuleProcessingError(f"Could not process rule {self.rule_string!r}: {e}") from e for name in self.rule.co_names: if name not in ns: raise RuleProcessingError( @@ -527,9 +520,7 @@ def __call__(self, dataset, idat): if k == "PTO" and hasattr(self, "PTO"): if v not in self.PTO: return None - elif hasattr(self, k) and ( - getattr(self, k) != v - ): + elif hasattr(self, k) and (getattr(self, k) != v): return None # Will return True if datapoint passes through the filter @@ -543,12 +534,10 @@ def __call__(self, dataset, idat): **ns, }, ) - except Exception as e: # pragma: no cover - raise FatalRuleError( - f"Error when applying rule {self.rule_string!r}: {e}" - ) from e + except Exception as e: # pragma: no cover + raise FatalRuleError(f"Error when applying rule {self.rule_string!r}: {e}") from e - def __repr__(self): # pragma: no cover + def __repr__(self): # pragma: no cover return self.rule_string def _make_kinematics_dict(self, dataset, idat) -> dict: @@ -565,6 +554,7 @@ def _make_point_namespace(self, dataset, idat) -> dict: ns[key] = eval(value, {**self.numpy_functions, **ns}) return ns + def get_cuts_for_dataset(commondata, rules) -> list: """Function to generate a list containing the index of all experimental points that passed kinematic diff --git a/validphys2/src/validphys/n3fit_data.py b/validphys2/src/validphys/n3fit_data.py index 75e5efdc12..38f946248d 100644 --- a/validphys2/src/validphys/n3fit_data.py +++ b/validphys2/src/validphys/n3fit_data.py @@ -4,8 +4,8 @@ Providers which prepare the data ready for :py:func:`n3fit.performfit.performfit`. """ -import functools from collections import defaultdict +import functools import hashlib import logging @@ -14,11 +14,8 @@ from reportengine import collect from reportengine.table import table - -from validphys.n3fit_data_utils import ( - validphys_group_extractor, -) from validphys.core import IntegrabilitySetSpec, TupleComp +from validphys.n3fit_data_utils import validphys_group_extractor log = logging.getLogger(__name__) @@ -53,6 +50,13 @@ def replica_mcseed(replica, mcseed, genrep): return res +def replica_luxseed(replica, luxseed): + """Generate the ``luxseed`` for a ``replica``. + Identical to replica_nnseed but used for a different purpose. + """ + return replica_nnseed(replica, luxseed) + + class _TrMasks(TupleComp): """Class holding the training validation mask for a group of datasets If the same group of dataset receives the same trvlseed then the mask @@ -91,13 +95,11 @@ def tr_masks(data, replica_trvlseed): ndata = len(cuts.load()) if cuts else dataset.commondata.ndata frac = dataset.frac # We do this so that a given dataset will always have the same number of points masked - trmax = int(ndata*frac) + trmax = int(ndata * frac) if trmax == 0: # If that number is 0, then get 1 point with probability frac trmax = int(rng.random() < frac) - mask = np.concatenate( - [np.ones(trmax, dtype=bool), np.zeros(ndata - trmax, dtype=bool)] - ) + mask = np.concatenate([np.ones(trmax, dtype=bool), np.zeros(ndata - trmax, dtype=bool)]) rng.shuffle(mask) trmask_partial.append(mask) return _TrMasks(str(data), replica_trvlseed, trmask_partial) @@ -236,7 +238,7 @@ def fitting_data_dict( expdata = make_replica tr_masks = tr_masks.masks - covmat = dataset_inputs_fitting_covmat # t0 covmat, or theory covmat or whatever was decided by the runcard + covmat = dataset_inputs_fitting_covmat # t0 covmat, or theory covmat or whatever was decided by the runcard inv_true = np.linalg.inv(covmat) fittable_datasets = fittable_datasets_masked @@ -345,7 +347,7 @@ def pseudodata_table(groups_replicas_indexed_make_replica, replicas): Notes ----- Whilst running ``n3fit``, this action will only be called if - `fitting::savepseudodata` is `true` (as per the default setting) and + `fitting::savepseudodata` is `true` (as per the default setting) and replicas are fitted one at a time. The table can be found in the replica folder i.e. /nnfit/replica_*/ diff --git a/validphys2/src/validphys/pdfbases.py b/validphys2/src/validphys/pdfbases.py index 7c9f4a3664..10ad6543ed 100644 --- a/validphys2/src/validphys/pdfbases.py +++ b/validphys2/src/validphys/pdfbases.py @@ -525,6 +525,9 @@ def f_(transform_func): EVOL = evolution +LUX = copy.deepcopy(evolution) +LUX.default_elements = (r'\Sigma', 'V', 'T3', 'V3', 'T8', 'V8', 'T15', 'gluon', 'photon') + CCBAR_ASYMM = copy.deepcopy(evolution) CCBAR_ASYMM.default_elements = (r'\Sigma', 'V', 'T3', 'V3', 'T8', 'V8', 'T15', 'gluon', 'V15') diff --git a/validphys2/src/validphys/photon/__init__.py b/validphys2/src/validphys/photon/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py new file mode 100644 index 0000000000..bbbdfc650f --- /dev/null +++ b/validphys2/src/validphys/photon/compute.py @@ -0,0 +1,334 @@ +"""Script that calls fiatlux to add the photon PDF.""" +import logging +import tempfile + +import fiatlux +import numpy as np +from scipy.integrate import trapezoid +from scipy.interpolate import interp1d +import yaml + +from eko.io import EKO +from n3fit.io.writer import XGRID +from validphys.n3fit_data import replica_luxseed + +from . import structure_functions as sf +from .constants import ED2, EU2, NC, NL + +log = logging.getLogger(__name__) + +Q_IN = 100 +FIATLUX_DEFAULT = { + "apfel": False, + "eps_base": 1e-5, # precision on final integration of double integral. + "eps_rel": 1e-1, # extra precision on any single integration. + "mum_proton": 2.792847356, # proton magnetic moment, from + # http://pdglive.lbl.gov/DataBlock.action?node=S016MM which itself + # gets it from arXiv:1203.5425 (CODATA) + # the elastic param type, options: + # dipole + # A1_world_spline + # A1_world_pol_spline + "elastic_param": "A1_world_pol_spline", + "elastic_electric_rescale": 1, + "elastic_magnetic_rescale": 1, + # the inelastic param type, options: + "inelastic_param": "LHAPDF_Hermes_ALLM_CLAS", # Hermes_ALLM_CLAS, LHAPDF_Hermes_ALLM_CLAS + "rescale_r_twist4": 0, + "rescale_r": 1, + "allm_limits": 0, + "rescale_non_resonance": 1, + "rescale_resonance": 1, + "use_mu2_as_upper_limit": False, + "q2min_inel_override": 0.0, + "q2max_inel_override": 1e300, + "lhapdf_transition_q2": 9, + # general + "verbose": False, +} + + +class Photon: + """Photon class computing the photon array with the LuxQED approach.""" + + 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 + fiatlux_runcard["mproton"] = theory["MP"] + self.replicas = replicas + + # structure functions + self.luxpdfset = lux_params["luxset"].load() + self.additional_errors = lux_params["additional_errors"] + self.luxseed = lux_params["luxseed"] + + # TODO : maybe find a different name for fiatlux_dis_F2 + path_to_F2 = theoryid.path / "fastkernel/fiatlux_dis_F2.pineappl.lz4" + path_to_FL = theoryid.path / "fastkernel/fiatlux_dis_FL.pineappl.lz4" + self.path_to_eko_photon = theoryid.path / "eko_photon.tar" + + # set fiatlux + self.lux = {} + + alpha = Alpha(theory) + mb_thr = theory["kbThr"] * theory["mb"] + mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100 + + self.interpolator = [] + self.integral = [] + + for replica in replicas: + f2 = sf.InterpStructureFunction(path_to_F2, self.luxpdfset.members[replica]) + fl = sf.InterpStructureFunction(path_to_FL, self.luxpdfset.members[replica]) + if not np.isclose(f2.q2_max, fl.q2_max): + log.error( + "FKtables for fiatlux_dis_F2 and fiatlux_dis_FL have two different q2_max" + ) + + fiatlux_runcard["q2_max"] = float(f2.q2_max) + f2lo = sf.F2LO(self.luxpdfset.members[replica], theory) + with tempfile.NamedTemporaryFile(mode="w") as tmp: + with tmp.file as tmp_file: + tmp_file.write(yaml.dump(fiatlux_runcard)) + self.lux[replica] = fiatlux.FiatLux(tmp_file.name) + # we have a dict but fiatlux wants a yaml file + # TODO : once that fiatlux will allow dictionaries + # pass directly fiatlux_runcard + + self.lux[replica].PlugAlphaQED(alpha.alpha_em, alpha.qref) + self.lux[replica].InsertInelasticSplitQ( + [ + mb_thr, + mt_thr, + ] + ) + self.lux[replica].PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq) + + photon_array = self.compute_photon_array(replica) + self.interpolator.append(interp1d(XGRID, photon_array, fill_value=0.0, kind="cubic")) + self.integral.append(trapezoid(photon_array, XGRID)) + + def compute_photon_array(self, replica): + r""" + Compute the photon PDF for every point in the grid xgrid. + + Parameters + ---------- + replica: int + replica id + + Returns + ------- + compute_photon_array: numpy.array + photon PDF at the fitting scale Q0 + """ + # Compute photon PDF + log.info(f"Computing photon") + photon_qin = np.array([self.lux[replica].EvaluatePhoton(x, Q_IN**2).total for x in XGRID]) + photon_qin += self.generate_errors(replica) + # fiatlux computes x * gamma(x) + photon_qin /= XGRID + # 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: + # 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.rotations.inputpids), len(XGRID))) + for j, pid in enumerate(eko.rotations.inputpids): + if pid == 22: + pdfs_init[j] = photon_qin + ph_id = j + else: + if pid not in self.luxpdfset.flavors: + continue + pdfs_init[j] = np.array( + [self.luxpdfset.xfxQ(x, Q_IN, replica, pid) / x for x in XGRID] + ) + + # Apply EKO to PDFs + q2 = eko.mu2grid[0] + with eko.operator(q2) as elem: + pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init) + + photon_Q0 = pdfs_final[ph_id] + + # we want x * gamma(x) + return XGRID * photon_Q0 + + def __call__(self, xgrid): + """ + Compute the photon interpolating the values of self.photon_array. + + Parameters + ---------- + xgrid : nd.array + array of x values with shape (1,xgrid,1) + + Returns + ------- + photon values : nd.array + array of photon values with shape (1,xgrid,1) + """ + return [ + self.interpolator[id](xgrid[0, :, 0])[np.newaxis, :, np.newaxis] + for id in range(len(self.replicas)) + ] + + @property + def error_matrix(self): + """Generate error matrix to be used in generate_errors.""" + if not self.additional_errors: + return None + extra_set = self.additional_errors.load() + qs = [Q_IN] * len(XGRID) + res_central = np.array(extra_set.central_member.xfxQ(22, XGRID, qs)) + res = [] + for idx_member in range(101, 107 + 1): + tmp = np.array(extra_set.members[idx_member].xfxQ(22, XGRID, qs)) + res.append(tmp - res_central) + # first index must be x, while second one must be replica index + return np.stack(res, axis=1) + + def generate_errors(self, replica): + """ + Generate LUX additional errors according to the procedure + described in sec. 2.5 of https://arxiv.org/pdf/1712.07053.pdf + """ + if self.error_matrix is None: + return np.zeros_like(XGRID) + log.info(f"Generating photon additional errors") + seed = replica_luxseed(replica, self.luxseed) + rng = np.random.default_rng(seed=seed) + u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False) + errors = u @ (s * rng.normal(size=7)) + return errors + + +class Alpha: + def __init__(self, theory): + self.theory = theory + self.alpha_em_ref = theory["alphaqed"] + self.qref = self.theory["Qedref"] + self.beta0, self.b1 = self.set_betas() + self.thresh, self.alpha_thresh = self.set_thresholds_alpha_em() + + def alpha_em(self, q): + r""" + Compute the value of alpha_em. + + Parameters + ---------- + q: float + value in which the coupling is computed + + Returns + ------- + alpha_em: float + electromagnetic coupling + """ + if q < self.thresh_c: + nf = 3 + elif q < self.thresh_b: + nf = 4 + elif q < self.thresh_t: + nf = 5 + else: + nf = 6 + return self.alpha_em_fixed_flavor(q, self.alpha_thresh[nf], self.thresh[nf], nf) + + def alpha_em_fixed_flavor(self, q, alpha_ref, qref, nf): + """ + Compute the alpha_em running for nf fixed at NLO. + + Parameters + ---------- + q : float + target scale + alph_aem_ref : float + reference value of alpha_em + qref: float + reference scale + nf: int + number of flavors + + Returns + ------- + alpha_em at NLO : float + target value of a + """ + lmu = 2 * np.log(q / qref) + den = 1.0 + self.beta0[nf] * alpha_ref * lmu + alpha_LO = alpha_ref / den + alpha_NLO = alpha_LO * (1 - self.b1[nf] * alpha_LO * np.log(den)) + return alpha_NLO + + def set_thresholds_alpha_em(self): + """ + Compute and store the couplings at thresholds to speed up the calling + to alpha_em inside fiatlux: + when q is in a certain range (e.g. thresh_c < q < thresh_b) and qref in a different one + (e.g. thresh_b < q < thresh_t) we need to evolve from qref to thresh_b with nf=5 and then + from thresh_b to q with nf=4. Given that the value of alpha at thresh_b is always the same + we can avoid computing the first step storing the values of alpha in the threshold points. + It is done for qref in a generic range (not necessarly qref=91.2). + + """ + 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: + self.thresh_t = np.inf + if self.theory["MaxNfAs"] <= 4: + self.thresh_b = np.inf + if self.theory["MaxNfAs"] <= 3: + self.thresh_c = np.inf + + thresh_list = [self.thresh_c, self.thresh_b, self.thresh_t] + # determine nfref + if self.qref < self.thresh_c: + nfref = 3 + elif self.qref < self.thresh_b: + nfref = 4 + elif self.qref < self.thresh_t: + nfref = 5 + else: + nfref = 6 + thresh_list.insert(nfref - 3, self.qref) + + thresh = {nf: thresh_list[nf - 3] for nf in range(3, self.theory["MaxNfAs"] + 1)} + + alpha_thresh = {nfref: self.alpha_em_ref} + + # determine the values of alpha in the threshold points, depending on the value of qref + for nf in range(nfref + 1, self.theory["MaxNfAs"] + 1): + alpha_thresh[nf] = self.alpha_em_fixed_flavor( + thresh[nf], alpha_thresh[nf - 1], thresh[nf - 1], nf - 1 + ) + + for nf in reversed(range(3, nfref)): + alpha_thresh[nf] = self.alpha_em_fixed_flavor( + thresh[nf], alpha_thresh[nf + 1], thresh[nf + 1], nf + 1 + ) + + return thresh, alpha_thresh + + def set_betas(self): + """Compute and store beta0 / 4pi and b1 = (beta1/beta0)/4pi as a function of nf.""" + beta0 = {} + b1 = {} + for nf in range(3, 6 + 1): + nu = nf // 2 + nd = nf - nu + beta0[nf] = (-4.0 / 3 * (NL + NC * (nu * EU2 + nd * ED2))) / (4 * np.pi) + b1[nf] = ( + -4.0 * (NL + NC * (nu * EU2**2 + nd * ED2**2)) / beta0[nf] / (4 * np.pi) ** 2 + ) + return beta0, b1 diff --git a/validphys2/src/validphys/photon/constants.py b/validphys2/src/validphys/photon/constants.py new file mode 100644 index 0000000000..d820fb12a9 --- /dev/null +++ b/validphys2/src/validphys/photon/constants.py @@ -0,0 +1,4 @@ +EU2 = 4.0 / 9 +ED2 = 1.0 / 9 +NL = 3 +NC = 3 diff --git a/validphys2/src/validphys/photon/structure_functions.py b/validphys2/src/validphys/photon/structure_functions.py new file mode 100644 index 0000000000..fe17e32ad8 --- /dev/null +++ b/validphys2/src/validphys/photon/structure_functions.py @@ -0,0 +1,115 @@ +from abc import ABC, abstractmethod + +import numpy as np +import pineappl +from scipy.interpolate import RectBivariateSpline + +from . import constants + + +class StructureFunction(ABC): + """Abstract class for the DIS structure functions""" + + @abstractmethod + def fxq(self, x, Q): + r""" + Abstract method returning the value of a DIS structure functions. + + Parameters + ---------- + x : float + Bjorken x + Q : float + DIS hard scale + + Returns + ------- + F_{2,L}: float + Structure function F2 or FL + """ + pass + + +class InterpStructureFunction(StructureFunction): + """ + Compute an interpolated structure function convoluting an FKtable + with a PDF. + """ + + def __init__(self, path_to_fktable, pdfs): + self.fktable = pineappl.fk_table.FkTable.read(path_to_fktable) + x = np.unique(self.fktable.bin_left(1)) + q2 = np.unique(self.fktable.bin_left(0)) + + self.q2_max = max(q2) + + predictions = self.fktable.convolute_with_one(2212, pdfs.xfxQ2) + + grid2D = predictions.reshape(len(x), len(q2)) + self.interpolator = RectBivariateSpline(x, q2, grid2D) + + def fxq(self, x, q): + r""" + Compute the DIS structure function interpolating the grid. + + Parameters + ---------- + x : float + Bjorken x + Q : float + DIS hard scale + + Returns + ------- + F_{2,L}: float + Structure function F2 or FL + """ + return self.interpolator(x, q**2)[0, 0] + + +class F2LO(StructureFunction): + """Analytically compute the structure function F2 at leading order.""" + + def __init__(self, pdfs, theory): + self.pdfs = pdfs + self.thresh_c = theory["kcThr"] * theory["mc"] + self.thresh_b = theory["kbThr"] * theory["mb"] + self.thresh_t = theory["ktThr"] * theory["mt"] + if theory["MaxNfPdf"] <= 5: + self.thresh_t = np.inf + if theory["MaxNfPdf"] <= 4: + self.thresh_b = np.inf + if theory["MaxNfPdf"] <= 3: + self.thresh_c = np.inf + self.eq2 = [constants.ED2, constants.EU2] * 3 # d u s c b t + + def fxq(self, x, q): + r""" + Analytically compute F2LO. + + Parameters + ---------- + x : float + Bjorken x + Q : float + DIS hard scale + + Returns + ------- + F_2^{LO} : float + Structure function F2 at LO + """ + # at LO we use ZM-VFS + if q < self.thresh_c: + nf = 3 + elif q < self.thresh_b: + nf = 4 + elif q < self.thresh_t: + nf = 5 + else: + nf = 6 + res = 0 + pdfs_values = self.pdfs.xfxQ(x, q) + for i in range(1, nf + 1): + res += self.eq2[i - 1] * (pdfs_values[i] + pdfs_values[-i]) + return res diff --git a/validphys2/src/validphys/tests/photon/__init__.py b/validphys2/src/validphys/tests/photon/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py new file mode 100644 index 0000000000..30805be889 --- /dev/null +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -0,0 +1,157 @@ +from collections import namedtuple +from pathlib import Path + +import fiatlux +import numpy as np +from validphys.photon import structure_functions +from validphys.photon.compute import Photon, Alpha +from validphys.core import PDF as PDFset + +from ..conftest import PDF + + +class FakeTheory: + def __init__(self): + self.path = Path("/fake/path/") + + def get_description(self): + return { + "alphaqed": 0.01, + "Qref": 91.2, + "Qedref": 91.2, + "mc": 1.3, + "mb": 4.92, + "mt": 172.0, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "MaxNfAs": 5, + "MaxNfPdf": 5, + "MP": 0.938 + } + + +fiatlux_runcard = { + "luxset": PDFset(PDF), + "additional_errors": False, + "luxseed": 123456789 +} + +photon = namedtuple("photon", ["total", "elastic", "inelastic"]) + + +class FakeFiatlux: + def __init__(self, runcard): + self.runcard = runcard + self.alphaem = None + self.qref = None + self.trash1 = None + self.trash2 = None + self.f2 = None + self.fl = None + self.f2lo = None + self.res = photon(0., 0., 0.) + + def PlugAlphaQED(self, alphaem, qref): + self.alphaem = alphaem + self.qref = qref + + def InsertInelasticSplitQ(self, args): + self.trash1 = args[0] + self.trash2 = args[1] + + def PlugStructureFunctions(self, f2, fl, f2lo): + self.f2 = f2 + self.fl = fl + self.f2lo = f2lo + + def EvaluatePhoton(self, x, q): + return self.res + + +class FakeStructureFunction: + def __init__(self, path, pdfs): + self.path = path + self.pdfs = pdfs + self.q2_max = 1e8 + + def fxq(self): + return 0 + + +class FakeF2LO: + def __init__(self, pdfs, theory): + self.pdfs = pdfs + self.theory = theory + + def fxq(self): + return 0 + + +def test_parameters_init(monkeypatch): + monkeypatch.setattr( + structure_functions, "InterpStructureFunction", FakeStructureFunction + ) + monkeypatch.setattr(structure_functions, "F2LO", FakeF2LO) + monkeypatch.setattr(fiatlux, "FiatLux", FakeFiatlux) + monkeypatch.setattr(Photon, "compute_photon_array", lambda *args: np.zeros(196)) + + photon = Photon(FakeTheory(), fiatlux_runcard, [1, 2, 3]) + alpha = Alpha(FakeTheory().get_description()) + + np.testing.assert_equal(photon.replicas, [1, 2, 3]) + np.testing.assert_equal(photon.luxpdfset._name, fiatlux_runcard["luxset"].name) + np.testing.assert_equal(photon.additional_errors, fiatlux_runcard["additional_errors"]) + np.testing.assert_equal(photon.luxseed, fiatlux_runcard["luxseed"]) + np.testing.assert_almost_equal( + alpha.alpha_em_ref, FakeTheory().get_description()["alphaqed"] + ) + +def test_masses_init(): + alpha = Alpha(FakeTheory().get_description()) + np.testing.assert_equal(alpha.thresh_t, np.inf) + np.testing.assert_almost_equal(alpha.thresh_b, 4.92) + np.testing.assert_almost_equal(alpha.thresh_c, 1.3) + +def test_set_thresholds_alpha_em(monkeypatch): + monkeypatch.setattr( + structure_functions, "InterpStructureFunction", FakeStructureFunction + ) + monkeypatch.setattr(structure_functions, "F2LO", FakeF2LO) + + monkeypatch.setattr(fiatlux, "FiatLux", FakeFiatlux) + monkeypatch.setattr(Photon, "compute_photon_array", lambda *args: np.zeros(196)) + + alpha = Alpha(FakeTheory().get_description()) + + np.testing.assert_almost_equal(alpha.thresh[5], 91.2) + np.testing.assert_almost_equal(alpha.thresh[4], 4.92) + np.testing.assert_almost_equal(alpha.thresh[3], 1.3) + np.testing.assert_almost_equal(alpha.alpha_thresh[5], 0.01) + np.testing.assert_almost_equal( + alpha.alpha_thresh[4], alpha.alpha_em_fixed_flavor(4.92, 0.01, 91.2, 5) + ) + np.testing.assert_almost_equal( + alpha.alpha_thresh[3], + alpha.alpha_em_fixed_flavor(1.3, alpha.alpha_thresh[4], 4.92, 4), + ) + np.testing.assert_equal(len(alpha.alpha_thresh), 3) + np.testing.assert_equal(len(alpha.thresh), 3) + +def test_betas(): + alpha = Alpha(FakeTheory().get_description()) + vec_beta0 = [ + -0.5305164769729844, + -0.6719875374991137, + -0.7073553026306458, + -0.8488263631567751, + ] + vec_b1 = [ + 0.17507043740108488, + 0.1605510390839295, + 0.1538497783221655, + 0.1458920311675707, + ] + for nf in range(3, 6 + 1): + np.testing.assert_allclose(alpha.beta0[nf], vec_beta0[nf - 3], rtol=1e-7) + np.testing.assert_allclose(alpha.b1[nf], vec_b1[nf - 3], rtol=1e-7) diff --git a/validphys2/src/validphys/tests/photon/test_structurefunctions.py b/validphys2/src/validphys/tests/photon/test_structurefunctions.py new file mode 100644 index 0000000000..984eb2eecd --- /dev/null +++ b/validphys2/src/validphys/tests/photon/test_structurefunctions.py @@ -0,0 +1,127 @@ +import numpy as np +import pineappl +import validphys.photon.structure_functions as sf +from validphys.lhapdfset import LHAPDFSet + + +class ZeroPdfs: + def xfxQ(self, x, Q): + res = {} + for i in range(1, 6 + 1): + res[i] = res[-i] = 0.0 + return res + + +def test_zero_pdfs(): + + pdfs = ZeroPdfs() + + fake_theory = { + "mc": 1.3, + "mb": 5.0, + "mt": 172.0, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "MaxNfPdf": 5, + } + + f2lo = sf.F2LO(pdfs, fake_theory) + + np.testing.assert_equal(f2lo.thresh_t, np.inf) + + for x in np.geomspace(1e-4, 1.0, 10): + for Q in np.geomspace(10, 1000000, 10): + np.testing.assert_allclose(f2lo.fxq(x, Q), 0.0) + + +class FakeSet: + def get_entry(self, string): + return 0 + + +class ZeroFKTable: + def __init__(self, path): + self.path = path + self.xgrid = np.geomspace(1e-4, 1.0, 10) + self.qgrid = np.geomspace(1.65, 1000, 10) + + def bin_left(self, i): + if i == 1: + return self.xgrid + if i == 0: + return self.qgrid + else: + return 0 + + def convolute_with_one(self, pdgid, xfxQ2): + return np.zeros((10, 10)) + + +class OnePdf: + def __init__(self): + self.ao = 1 + + def xfxQ(self, x, Q): + return 1.0 + + def xfxQ2(self, x, Q): + return 1.0**2 + + def set(self): + return FakeSet() + + +class OneFKTable: + def __init__(self, path): + self.path = path + self.xgrid = np.geomspace(1e-4, 1.0, 10) + self.qgrid = np.geomspace(1.65, 1000, 10) + + def bin_left(self, i): + if i == 1: + return self.xgrid + if i == 0: + return self.qgrid + else: + return 0 + + def convolute_with_one(self, pdgid, xfxQ2): + return np.zeros((10, 10)) + + +class ZeroPdf: + def __init__(self): + self.ao = 1 + + def xfxQ(self, x, Q): + return 1.0 + + def xfxQ2(self, x, Q): + return 1.0**2 + + def set(self): + return FakeSet() + + +def test_F2(monkeypatch): + # grid put to 0 and pdf put to 1 + monkeypatch.setattr(pineappl.fk_table.FkTable, "read", ZeroFKTable) + structurefunc = sf.InterpStructureFunction("", OnePdf()) + for x in np.geomspace(1e-4, 1.0, 10): + for Q in np.geomspace(10, 1000000, 10): + np.testing.assert_allclose(structurefunc.fxq(x, Q), 0.0, rtol=1e-5) + + # grid put to 0 and real pdf + pdf = LHAPDFSet("NNPDF40_nnlo_as_01180", "replicas") + structurefunc = sf.InterpStructureFunction("", pdf.central_member) + for x in np.geomspace(1e-4, 1.0, 10): + for Q in np.geomspace(10, 1000000, 10): + np.testing.assert_allclose(structurefunc.fxq(x, Q), 0.0, rtol=1e-5) + + # grid put to 1 and pdf put to 0 + monkeypatch.setattr(pineappl.fk_table.FkTable, "read", OneFKTable) + structurefunc = sf.InterpStructureFunction("", ZeroPdf()) + for x in np.geomspace(1e-4, 1.0, 10): + for Q in np.geomspace(10, 1000000, 10): + np.testing.assert_allclose(structurefunc.fxq(x, Q), 0.0, rtol=1e-5)