diff --git a/n3fit/runcards/examples/Basic_runcard_qed.yml b/n3fit/runcards/examples/Basic_runcard_qed.yml index 309f5901d5..9fd98ce1b1 100644 --- a/n3fit/runcards/examples/Basic_runcard_qed.yml +++ b/n3fit/runcards/examples/Basic_runcard_qed.yml @@ -95,7 +95,7 @@ datacuts: ############################################################ theory: - theoryid: 522 # database id + theoryid: 523 # database id ############################################################ trvlseed: 1551864071 @@ -169,3 +169,4 @@ fiatlux: luxset: NNPDF40_nnlo_as_01180 additional_errors: true # should be set to true only for the last iteration luxseed: 1234567890 + eps_base: 1e-2 diff --git a/n3fit/src/n3fit/tests/regressions/quickcard_qed.yml b/n3fit/src/n3fit/tests/regressions/quickcard_qed.yml new file mode 100644 index 0000000000..84a880c44e --- /dev/null +++ b/n3fit/src/n3fit/tests/regressions/quickcard_qed.yml @@ -0,0 +1,93 @@ +# +# Configuration file for n3fit regression tests +# This runcard includes two DIS datasets, one Hadronic dataset +# and two positivity datasets +# + +############################################################ +description: n3fit regression test + +############################################################ +# frac: training fraction +# ewk: apply ewk k-factors +# sys: systematics treatment (see systypes) +dataset_inputs: +- { dataset: NMC, frac: 0.5 } +- { dataset: SLACP_dwsh, frac: 0.5} +- { dataset: CMSZDIFF12, frac: 0.5, cfac: ['QCD'], sys: 10 } +- { dataset: ATLASTTBARTOT8TEV, frac: 1.0, cfac: ['QCD'] } + +############################################################ +datacuts: + t0pdfset: NNPDF40_nnlo_as_01180 # PDF set to generate t0 covmat + q2min : 3.49 # Q2 minimum + w2min : 12.5 # W2 minimum + combocuts : NNPDF31 # NNPDF3.0 final kin. cuts + jetptcut_tev : 0 # jet pt cut for tevatron + jetptcut_lhc : 0 # jet pt cut for lhc + wptcut_lhc : 30.0 # Minimum pT for W pT diff distributions + jetycut_tev : 1e30 # jet rap. cut for tevatron + jetycut_lhc : 1e30 # jet rap. cut for lhc + dymasscut_min: 0 # dy inv.mass. min cut + dymasscut_max: 1e30 # dy inv.mass. max cut + jetcfactcut : 1e30 # jet cfact. cut + +############################################################ +theory: + theoryid: 398 # database id + +############################################################ +genrep: True # on = generate MC replicas, False = use real data +trvlseed: 3 +nnseed: 2 +mcseed: 1 + +load: "weights.h5" + +parameters: # This defines the parameter dictionary that is passed to the Model Trainer + nodes_per_layer: [15, 10, 8] + activation_per_layer: ['sigmoid', 'sigmoid', 'linear'] + initializer: 'glorot_normal' + optimizer: + optimizer_name: 'RMSprop' + learning_rate: 0.00001 + clipnorm: 1.0 + epochs: 1100 + positivity: + multiplier: 1.05 + initial: 1.5 + stopping_patience: 0.10 # percentage of the number of epochs + layer_type: 'dense' + dropout: 0.0 + threshold_chi2: 10.0 + +fitting: + fitbasis: NN31IC # EVOL (7), EVOLQED (8), etc. + basis: + - { fl: sng, smallx: [1.05,1.19], largex: [1.47,2.70] } + - { fl: g, smallx: [0.94,1.25], largex: [0.11,5.87] } + - { fl: v, smallx: [0.54,0.75], largex: [1.15,2.76] } + - { fl: v3, smallx: [0.21,0.57], largex: [1.35,3.08] } + - { fl: v8, smallx: [0.52,0.76], largex: [0.77,3.56] } + - { fl: t3, smallx: [-0.37,1.52], largex: [1.74,3.39] } + - { fl: t8, smallx: [0.56,1.29], largex: [1.45,3.03] } + - { fl: cp, smallx: [0.12,1.19], largex: [1.83,6.70] } + +############################################################ +positivity: + posdatasets: + - { dataset: POSF2U, maxlambda: 1e6 } # Positivity Lagrange Multiplier + - { dataset: POSDYS, maxlambda: 1e5 } + +integrability: + integdatasets: + - {dataset: INTEGXT8, maxlambda: 1e2} + +############################################################ +debug: true + +fiatlux: + luxset: NNPDF40_nnlo_as_01180 + additional_errors: true # should be set to true only for the last iteration + luxseed: 1234567890 + eps_base: 1e-2 diff --git a/n3fit/src/n3fit/tests/regressions/quickcard_qed_1.json b/n3fit/src/n3fit/tests/regressions/quickcard_qed_1.json new file mode 100644 index 0000000000..a960f7da9f --- /dev/null +++ b/n3fit/src/n3fit/tests/regressions/quickcard_qed_1.json @@ -0,0 +1,95 @@ +{ + "preprocessing": [ + { + "fl": "sng", + "smallx": 1.1307647228240967, + "largex": 2.6348154544830322, + "trainable": true + }, + { + "fl": "g", + "smallx": 1.1853630542755127, + "largex": 1.5627975463867188, + "trainable": true + }, + { + "fl": "v", + "smallx": 0.5399999022483826, + "largex": 2.004500150680542, + "trainable": true + }, + { + "fl": "v3", + "smallx": 0.3061824142932892, + "largex": 2.624323606491089, + "trainable": true + }, + { + "fl": "v8", + "smallx": 0.5774596929550171, + "largex": 2.120253801345825, + "trainable": true + }, + { + "fl": "t3", + "smallx": 1.3441987037658691, + "largex": 1.7566683292388916, + "trainable": true + }, + { + "fl": "t8", + "smallx": 1.04995858669281, + "largex": 1.945939064025879, + "trainable": true + }, + { + "fl": "cp", + "smallx": 0.7400740385055542, + "largex": 3.461853504180908, + "trainable": true + } + ], + "stop_epoch": 1100, + "best_epoch": 1099, + "erf_tr": 31.486101150512695, + "erf_vl": 28.4218692779541, + "chi2": 19.611825942993164, + "pos_state": "POS_VETO", + "arc_lengths": [ + 1.035533162554303, + 1.1953713068471692, + 1.0881025884095246, + 1.3414978876721764, + 1.0839843290638607 + ], + "integrability": [ + 0.002653829054906187, + 0.002653829054905521, + 0.0002567724250169823, + 3.2896786332130423, + 0.003927801561079747 + ], + "timing": { + "walltime": { + "Total": 58.37841296195984, + "start": 0.0, + "replica_set": 0.3208029270172119, + "replica_fitted": 58.37802076339722, + "replica_set_to_replica_fitted": 58.057217836380005 + }, + "cputime": { + "Total": 61.517351913999995, + "start": 0.0, + "replica_set": 2.813359167999998, + "replica_fitted": 61.516948647, + "replica_set_to_replica_fitted": 58.703589479 + } + }, + "version": { + "keras": "2.11.0", + "tensorflow": "2.11.0, mkl=False", + "numpy": "1.22.4", + "nnpdf": "4.0.6.846+g0a926fdf1", + "validphys": "4.0.6.846+g0a926fdf1" + } +} diff --git a/n3fit/src/n3fit/tests/regressions/quickcard_qed_2.json b/n3fit/src/n3fit/tests/regressions/quickcard_qed_2.json new file mode 100644 index 0000000000..18e6a605ae --- /dev/null +++ b/n3fit/src/n3fit/tests/regressions/quickcard_qed_2.json @@ -0,0 +1,95 @@ +{ + "preprocessing": [ + { + "fl": "sng", + "smallx": 1.1090763807296753, + "largex": 2.683891534805298, + "trainable": true + }, + { + "fl": "g", + "smallx": 0.9399998784065247, + "largex": 1.6250606775283813, + "trainable": true + }, + { + "fl": "v", + "smallx": 0.7499998211860657, + "largex": 1.7267229557037354, + "trainable": true + }, + { + "fl": "v3", + "smallx": 0.21201331913471222, + "largex": 1.3532426357269287, + "trainable": true + }, + { + "fl": "v8", + "smallx": 0.7599998712539673, + "largex": 2.390087127685547, + "trainable": true + }, + { + "fl": "t3", + "smallx": 1.4388009309768677, + "largex": 2.2958621978759766, + "trainable": true + }, + { + "fl": "t8", + "smallx": 1.0386217832565308, + "largex": 1.7531665563583374, + "trainable": true + }, + { + "fl": "cp", + "smallx": 0.24638080596923828, + "largex": 2.7976953983306885, + "trainable": true + } + ], + "stop_epoch": 1100, + "best_epoch": 1099, + "erf_tr": 3.5809452533721924, + "erf_vl": 3.6826603412628174, + "chi2": 2.141340732574463, + "pos_state": "POS_VETO", + "arc_lengths": [ + 1.3194883264768875, + 1.1995146017416334, + 1.054266019685804, + 5.074692492247958, + 1.0689068380364566 + ], + "integrability": [ + 0.02917606895789332, + 0.029176068957894374, + 0.0003924771135637162, + 12.764093160629272, + 0.02982003148645135 + ], + "timing": { + "walltime": { + "Total": 57.77593541145325, + "start": 0.0, + "replica_set": 0.31868839263916016, + "replica_fitted": 57.77560114860535, + "replica_set_to_replica_fitted": 57.45691275596619 + }, + "cputime": { + "Total": 60.771742980000006, + "start": 0.0, + "replica_set": 2.6913169609999983, + "replica_fitted": 60.771408723, + "replica_set_to_replica_fitted": 58.080091762 + } + }, + "version": { + "keras": "2.11.0", + "tensorflow": "2.11.0, mkl=False", + "numpy": "1.22.4", + "nnpdf": "4.0.6.846+g0a926fdf1", + "validphys": "4.0.6.846+g0a926fdf1" + } +} diff --git a/n3fit/src/n3fit/tests/test_evolven3fit.py b/n3fit/src/n3fit/tests/test_evolven3fit.py index 8cab0b4927..9dd372c60e 100644 --- a/n3fit/src/n3fit/tests/test_evolven3fit.py +++ b/n3fit/src/n3fit/tests/test_evolven3fit.py @@ -161,12 +161,13 @@ def test_eko_utils(tmp_path): assert_allclose(list(eko_op.operator_card.raw["mugrid"]), op_card_dict["mugrid"]) -@pytest.mark.parametrize("fitname", ["Basic_runcard_3replicas_lowprec_399"]) +@pytest.mark.parametrize("fitname", ["Basic_runcard_3replicas_lowprec_399", "Basic_runcard_qed_3replicas_lowprec_398"]) def test_perform_evolution(tmp_path, fitname): """Test that evolven3fit_new is able to utilize the current eko in the respective theory. In addition checks that the generated .info files are correct """ fit = API.fit(fit=fitname) + _ = API.theoryid(theoryid=fit.as_input()['theory']['theoryid']) # Move the fit to a temporary folder tmp_fit = tmp_path / fitname shutil.copytree(fit.path, tmp_fit) @@ -185,3 +186,4 @@ def test_perform_evolution(tmp_path, fitname): info = check_lhapdf_info(tmp_info) for datpath in tmp_nnfit.glob("replica_*/*.dat"): check_lhapdf_dat(datpath, info) + diff --git a/n3fit/src/n3fit/tests/test_fit.py b/n3fit/src/n3fit/tests/test_fit.py index 80a329faa0..97d7d55175 100644 --- a/n3fit/src/n3fit/tests/test_fit.py +++ b/n3fit/src/n3fit/tests/test_fit.py @@ -33,6 +33,7 @@ log = logging.getLogger(__name__) REGRESSION_FOLDER = pathlib.Path(__file__).with_name("regressions") QUICKNAME = "quickcard" +QUICKNAME_QED = "quickcard_qed" EXE = "n3fit" REPLICA = "1" EXPECTED_MAX_FITTIME = 130 # seen mac ~ 180 and linux ~ 90 @@ -64,22 +65,22 @@ def test_initialize_seeds(): assert len({replica_mcseed(rep, 1, True) for rep in same_replicas}) == 1 -def auxiliary_performfit(tmp_path, replica=1, timing=True, rel_error=2e-3): +def auxiliary_performfit(tmp_path, runcard=QUICKNAME, replica=1, timing=True, rel_error=2e-3): """Fits quickcard and checks the json file to ensure the results have not changed. """ - quickcard = f"{QUICKNAME}.yml" + quickcard = f"{runcard}.yml" # Prepare the runcard quickpath = REGRESSION_FOLDER / quickcard weightpath = REGRESSION_FOLDER / f"weights_{replica}.h5" # read up the previous json file for the given replica - old_json = load_data(REGRESSION_FOLDER / f"{QUICKNAME}_{replica}.json") + old_json = load_data(REGRESSION_FOLDER / f"{runcard}_{replica}.json") # cp runcard and weights to tmp folder shutil.copy(quickpath, tmp_path) shutil.copy(weightpath, tmp_path / "weights.h5") # run the fit sp.run(f"{EXE} {quickcard} {replica}".split(), cwd=tmp_path, check=True) # read up json files - full_json = tmp_path / f"{QUICKNAME}/nnfit/replica_{replica}/{QUICKNAME}.json" + full_json = tmp_path / f"{runcard}/nnfit/replica_{replica}/{runcard}.json" new_json = load_data(full_json) # Now compare to regression results, taking into account precision won't be 100% equal_checks = ["stop_epoch", "pos_state"] @@ -101,14 +102,16 @@ def auxiliary_performfit(tmp_path, replica=1, timing=True, rel_error=2e-3): @pytest.mark.darwin -def test_performfit(tmp_path): - auxiliary_performfit(tmp_path, replica=2, timing=False, rel_error=1e-1) +@pytest.mark.parametrize("runcard", [QUICKNAME, QUICKNAME_QED]) +def test_performfit(tmp_path, runcard): + auxiliary_performfit(tmp_path, runcard=runcard, replica=2, timing=False, rel_error=1e-1) @pytest.mark.linux @pytest.mark.parametrize("replica", [1, 2]) -def test_performfit_and_timing(tmp_path, replica): - auxiliary_performfit(tmp_path, replica=replica, timing=True) +@pytest.mark.parametrize("runcard", [QUICKNAME, QUICKNAME_QED]) +def test_performfit_and_timing(tmp_path, runcard, replica): + auxiliary_performfit(tmp_path, runcard=runcard, replica=replica, timing=True) @pytest.mark.skip(reason="Still not implemented in parallel mode") diff --git a/nnpdfcpp/data/theory.db b/nnpdfcpp/data/theory.db index 4bc4113daa..fa202758df 100644 Binary files a/nnpdfcpp/data/theory.db and b/nnpdfcpp/data/theory.db differ diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 657e807ba2..52c3b04c78 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -17,9 +17,9 @@ log = logging.getLogger(__name__) +# not the complete fiatlux runcard since some parameters are set in the code 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 @@ -59,6 +59,15 @@ def __init__(self, theoryid, lux_params, replicas): # This is going to be changed in favor of a bool em_running # in the runcard fiatlux_runcard["mproton"] = theory["MP"] + + # precision on final integration of double integral + if "eps_base" in lux_params: + fiatlux_runcard["eps_base"] = lux_params["eps_base"] + log.warning(f"Using fiatlux parameter eps_base from runcard") + else: + fiatlux_runcard["eps_base"] = 1e-5 + log.info(f"Using default value for fiatlux parameter eps_base") + self.replicas = replicas # structure functions @@ -112,7 +121,7 @@ def __init__(self, theoryid, lux_params, replicas): 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.interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic")) self.integral.append(trapezoid(photon_array, XGRID)) def compute_photon_array(self, replica): diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index f61b13a6bd..9e3226062e 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -1,158 +1,82 @@ -from collections import namedtuple -from pathlib import Path +import tempfile import fiatlux import numpy as np -from validphys.photon import structure_functions -from validphys.photon.compute import Photon, Alpha +import yaml + +from eko.io import EKO +from n3fit.io.writer import XGRID +from validphys.api import API from validphys.core import PDF as PDFset +from validphys.loader import FallbackLoader +from validphys.photon import structure_functions as sf +from validphys.photon.compute import FIATLUX_DEFAULT, Alpha, Photon from ..conftest import PDF -from eko.io import EKO +TEST_THEORY = API.theoryid(theoryid=398) -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 = { +FIATLUX_RUNCARD = { "luxset": PDFset(PDF), - "additional_errors": False, - "luxseed": 123456789 + # check if "LUXqed17_plus_PDF4LHC15_nnlo_100" is installed + "additional_errors": FallbackLoader().check_pdf("LUXqed17_plus_PDF4LHC15_nnlo_100"), + "luxseed": 123456789, + "eps_base": 1e-2, # using low precision to speed up tests } -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 FakeEKO: - def __init__(self, path): - self.path = path - self.mu20 = 100**2 - - def __enter__(self): - return self - - def __exit__(self, exc_type: type, _exc_value, _traceback): - pass - - -class FakeStructureFunction: - def __init__(self, path, pdfs): - self.path = path - self.pdfs = pdfs - self.q2_max = 1e8 - def fxq(self): - return 0 +def test_parameters_init(): + "test initailization of the parameters from Photon class" + fiatlux_runcard = FIATLUX_RUNCARD.copy() + # we are not testing the photon here so we make it faster + fiatlux_runcard['eps_base'] = 1e-1 -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)) - monkeypatch.setattr(EKO, "read", FakeEKO) - - photon = Photon(FakeTheory(), fiatlux_runcard, [1, 2, 3]) - alpha = Alpha(FakeTheory().get_description()) + photon = Photon(TEST_THEORY, fiatlux_runcard, [1, 2, 3]) 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"] - ) + np.testing.assert_equal(photon.luxpdfset._name, FIATLUX_RUNCARD["luxset"].name) + np.testing.assert_equal(photon.additional_errors.name, "LUXqed17_plus_PDF4LHC15_nnlo_100") + np.testing.assert_equal(photon.luxseed, FIATLUX_RUNCARD["luxseed"]) + np.testing.assert_equal(photon.path_to_eko_photon, TEST_THEORY.path / "eko_photon.tar") + np.testing.assert_equal(photon.q_in, 100.0) + def test_masses_init(): - alpha = Alpha(FakeTheory().get_description()) + "test thresholds values in Alpha class" + theory = TEST_THEORY.get_description() + alpha = Alpha(theory) 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) + np.testing.assert_almost_equal(alpha.thresh_b, theory["mb"]) + np.testing.assert_almost_equal(alpha.thresh_c, theory["mc"]) -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)) +def test_set_thresholds_alpha_em(): + "test value of alpha_em at threshold values" + theory = TEST_THEORY.get_description() - alpha = Alpha(FakeTheory().get_description()) + alpha = Alpha(theory) - 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_em_ref, theory["alphaqed"]) + np.testing.assert_almost_equal(alpha.thresh[5], theory["Qedref"]) + np.testing.assert_almost_equal(alpha.thresh[4], theory["mb"]) + np.testing.assert_almost_equal(alpha.thresh[3], theory["mc"]) + np.testing.assert_almost_equal(alpha.alpha_thresh[5], theory["alphaqed"]) np.testing.assert_almost_equal( - alpha.alpha_thresh[4], alpha.alpha_em_fixed_flavor(4.92, 0.01, 91.2, 5) + alpha.alpha_thresh[4], + alpha.alpha_em_fixed_flavor(theory["mb"], theory["alphaqed"], theory["Qedref"], 5), ) np.testing.assert_almost_equal( alpha.alpha_thresh[3], - alpha.alpha_em_fixed_flavor(1.3, alpha.alpha_thresh[4], 4.92, 4), + alpha.alpha_em_fixed_flavor(theory["mc"], alpha.alpha_thresh[4], theory["mb"], 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()) + "test betas for different nf" + alpha = Alpha(TEST_THEORY.get_description()) vec_beta0 = [ -0.5305164769729844, -0.6719875374991137, @@ -168,3 +92,74 @@ def test_betas(): 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) + + +def test_photon(): + """test that photon coming out of Photon interpolator matches the photon array + for XGRID points + """ + fiatlux_runcard = FIATLUX_RUNCARD.copy() + fiatlux_runcard["additional_errors"] = False + theory = TEST_THEORY.get_description() + + for replica in [1, 2, 3]: + photon = Photon(TEST_THEORY, fiatlux_runcard, [replica]) + + # set up fiatlux + path_to_F2 = TEST_THEORY.path / "fastkernel/fiatlux_dis_F2.pineappl.lz4" + path_to_FL = TEST_THEORY.path / "fastkernel/fiatlux_dis_FL.pineappl.lz4" + pdfs = FIATLUX_RUNCARD["luxset"].load() + f2 = sf.InterpStructureFunction(path_to_F2, pdfs.members[replica]) + fl = sf.InterpStructureFunction(path_to_FL, pdfs.members[replica]) + f2lo = sf.F2LO(pdfs.members[replica], theory) + + # runcard + fiatlux_default = FIATLUX_DEFAULT.copy() + fiatlux_default['mproton'] = theory['MP'] + fiatlux_default["qed_running"] = bool(np.isclose(theory["Qedref"], theory["Qref"])) + fiatlux_default["q2_max"] = float(f2.q2_max) + fiatlux_default["eps_base"] = FIATLUX_RUNCARD["eps_base"] + + # load fiatlux + with tempfile.NamedTemporaryFile(mode="w") as tmp: + with tmp.file as tmp_file: + tmp_file.write(yaml.dump(FIATLUX_DEFAULT)) + lux = fiatlux.FiatLux(tmp.name) + + alpha = Alpha(theory) + + lux.PlugAlphaQED(alpha.alpha_em, alpha.qref) + lux.InsertInelasticSplitQ( + [ + theory["kbThr"] * theory["mb"], + theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100, + ] + ) + lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq) + path_to_eko_photon = TEST_THEORY.path / "eko_photon.tar" + with EKO.read(path_to_eko_photon) as eko: + photon_fiatlux_qin = np.array([lux.EvaluatePhoton(x, eko.mu20).total for x in XGRID]) + photon_fiatlux_qin /= XGRID + # construct PDFs + pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID))) + for j, pid in enumerate(eko.bases.inputpids): + if pid == 22: + pdfs_init[j] = photon_fiatlux_qin + ph_id = j + else: + if pid not in pdfs.flavors: + continue + pdfs_init[j] = np.array( + [pdfs.xfxQ(x, np.sqrt(eko.mu20), replica, pid) / x for x in XGRID] + ) + + # Apply EKO to PDFs + for _, elem in eko.items(): + pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init) + + photon_Q0 = pdfs_final[ph_id] + photon_fiatlux = XGRID * photon_Q0 + + photon_validphys = photon(XGRID[np.newaxis, :, np.newaxis])[0][0, :, 0] + + np.testing.assert_allclose(photon_fiatlux, photon_validphys, rtol=1e-7) diff --git a/validphys2/src/validphys/tests/photon/test_structurefunctions.py b/validphys2/src/validphys/tests/photon/test_structurefunctions.py index 984eb2eecd..571ae765f2 100644 --- a/validphys2/src/validphys/tests/photon/test_structurefunctions.py +++ b/validphys2/src/validphys/tests/photon/test_structurefunctions.py @@ -1,7 +1,13 @@ import numpy as np import pineappl + +from validphys.api import API +from validphys.core import PDF as PDFset import validphys.photon.structure_functions as sf -from validphys.lhapdfset import LHAPDFSet + +from ..conftest import PDF + +TEST_THEORY = API.theoryid(theoryid=398) class ZeroPdfs: @@ -9,35 +15,12 @@ def xfxQ(self, x, Q): res = {} for i in range(1, 6 + 1): res[i] = res[-i] = 0.0 + res[21] = 0.0 + res[22] = 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 + def xfxQ2(self, i, x, Q2): + return self.xfxQ(x, np.sqrt(Q2))[i] class ZeroFKTable: @@ -58,70 +41,69 @@ 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 test_zero_pdfs(): + "test that a zero PDF gives a zero structure function" + pdfs = ZeroPdfs() + theory = TEST_THEORY.get_description() + path_to_F2 = TEST_THEORY.path / "fastkernel/fiatlux_dis_F2.pineappl.lz4" + path_to_FL = TEST_THEORY.path / "fastkernel/fiatlux_dis_FL.pineappl.lz4" - def xfxQ(self, x, Q): - return 1.0 + f2 = sf.InterpStructureFunction(path_to_F2, pdfs) + fl = sf.InterpStructureFunction(path_to_FL, pdfs) + f2lo = sf.F2LO(pdfs, theory) - def xfxQ2(self, x, Q): - return 1.0**2 + np.testing.assert_equal(f2lo.thresh_t, np.inf) - def set(self): - return FakeSet() + 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) + np.testing.assert_allclose(f2.fxq(x, Q), 0.0) + np.testing.assert_allclose(fl.fxq(x, Q), 0.0) -def test_F2(monkeypatch): - # grid put to 0 and pdf put to 1 +def test_zero_grid(monkeypatch): + "test that a zero grid gives a zero structure function" + # patching pineappl.fk_table.FkTable to use ZeroFKTable monkeypatch.setattr(pineappl.fk_table.FkTable, "read", ZeroFKTable) - structurefunc = sf.InterpStructureFunction("", OnePdf()) + pdfs = PDFset(PDF).load() + structurefunc = sf.InterpStructureFunction("", pdfs.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 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) +def test_params(): + "test initialization of parameters" + pdfs = PDFset(PDF).load() + replica = 1 + theory = TEST_THEORY.get_description() + for channel in ["F2", "FL"]: + tmp = "fastkernel/fiatlux_dis_" + channel + ".pineappl.lz4" + path_to_fktable = TEST_THEORY.path / tmp + struct_func = sf.InterpStructureFunction(path_to_fktable, pdfs.members[replica]) + np.testing.assert_allclose(struct_func.q2_max, 1e8) + f2lo = sf.F2LO(pdfs.members[replica], theory) + np.testing.assert_allclose(f2lo.thresh_c, theory["mc"]) + np.testing.assert_allclose(f2lo.thresh_b, theory["mb"]) + np.testing.assert_allclose(f2lo.thresh_t, np.inf) + + +def test_interpolation_grid(): + """test that the values coming out of InterpStructureFunction match the grid ones""" + pdfs = PDFset(PDF).load() + for replica in [1, 2, 3]: + for channel in ["F2", "FL"]: + tmp = "fastkernel/fiatlux_dis_" + channel + ".pineappl.lz4" + path_to_fktable = TEST_THEORY.path / tmp + fktable = pineappl.fk_table.FkTable.read(path_to_fktable) + x = np.unique(fktable.bin_left(1)) + q2 = np.unique(fktable.bin_left(0)) + predictions = fktable.convolute_with_one(2212, pdfs.members[replica].xfxQ2) + grid2D = predictions.reshape(len(x), len(q2)) + + struct_func = sf.InterpStructureFunction(path_to_fktable, pdfs.members[replica]) + for i, x_ in enumerate(x): + for j, q2_ in enumerate(q2): + np.testing.assert_allclose( + struct_func.fxq(x_, np.sqrt(q2_)), grid2D[i, j], rtol=1e-5 + )