diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index fb9edea8bf..d2c41ea803 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -51,7 +51,7 @@ requirements: - sphinxcontrib-bibtex - curio >=1.0 - pineappl >=0.5.8 - - eko >=0.13.4,<0.14 + - eko >=0.13.5,<0.14 - banana-hep >=0.6.8 - fiatlux 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/evolven3fit_new/eko_utils.py b/n3fit/src/evolven3fit_new/eko_utils.py index 7b6be83110..d636b82354 100644 --- a/n3fit/src/evolven3fit_new/eko_utils.py +++ b/n3fit/src/evolven3fit_new/eko_utils.py @@ -22,7 +22,7 @@ } EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = { - "ev_op_iterations": 10, + "ev_op_iterations": 30, "ev_op_max_order": (1, 0), "evolution_method": "iterate-exact", "inversion_method": "exact", @@ -30,7 +30,6 @@ } NFREF_DEFAULT = 5 -NF0_DEFAULT = 4 def construct_eko_cards( @@ -58,46 +57,77 @@ def construct_eko_cards( theory = Loader().check_theoryID(theoryID).get_description() theory.pop("FNS") theory.update(theory_card_dict) - if "nfref" not in theory: - theory["nfref"] = NFREF_DEFAULT - if "nf0" not in theory: - theory["nf0"] = NF0_DEFAULT - + # Prepare the thresholds according to MaxNfPdf thresholds = {"c": theory["kcThr"], "b": theory["kbThr"], "t": theory["ktThr"]} if theory["MaxNfPdf"] < 5: thresholds["b"] = np.inf if theory["MaxNfPdf"] < 6: thresholds["t"] = np.inf - + + if "nfref" not in theory: + theory["nfref"] = NFREF_DEFAULT + + # Set nf_0 according to the fitting scale unless set explicitly + mu0 = theory["Q0"] + if "nf0" not in theory: + if mu0 < theory["mc"] * thresholds["c"]: + theory["nf0"] = 3 + elif mu0 < theory["mb"] * thresholds["b"]: + theory["nf0"] = 4 + elif mu0 < theory["mt"] * thresholds["t"]: + theory["nf0"] = 5 + else: + theory["nf0"] = 6 + # Setting the thresholds in the theory card to inf if necessary - theory.update({"kbThr":thresholds["b"], "ktThr":thresholds["t"] }) - + theory.update({"kbThr": thresholds["b"], "ktThr": thresholds["t"]}) + # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory legacy_class = runcards.Legacy(theory, {}) theory_card = legacy_class.new_theory - # construct operator card + # Generate the q2grid, if q_fin and q_points are None, use `nf0` to select a default q2_grid = utils.generate_q2grid( - theory["Q0"], + mu0, q_fin, q_points, - {theory["mb"]: thresholds["b"], theory["mt"]: thresholds["t"]}, + { + theory["mc"]: thresholds["c"], + theory["mb"]: thresholds["b"], + theory["mt"]: thresholds["t"], + }, + theory["nf0"], ) + + # construct operator card op_card = default_op_card masses = np.array([theory["mc"], theory["mb"], theory["mt"]]) ** 2 thresholds_ratios = np.array([thresholds["c"], thresholds["b"], thresholds["t"]]) ** 2 atlas = Atlas( matching_scales=MatchingScales(masses * thresholds_ratios), - origin=(theory["Q0"] ** 2, theory["nf0"]), - ) - op_card.update( - { - "mu0": theory["Q0"], - "mugrid": [(float(np.sqrt(q2)), int(nf_default(q2, atlas))) for q2 in q2_grid], - } + origin=(mu0**2, theory["nf0"]), ) + + # Create the eko operator q2grid + # This is a grid which contains information on (q, nf) + # in VFNS values at the matching scales need to be doubled so that they are considered in both sides + + ep = 1e-4 + mugrid = [] + for q2 in q2_grid: + q = float(np.sqrt(q2)) + if nf_default(q2 + ep, atlas) != nf_default(q2 - ep, atlas): + nf_l = int(nf_default(q2 - ep, atlas)) + nf_u = int(nf_default(q2 + ep, atlas)) + mugrid.append((q, nf_l)) + mugrid.append((q, nf_u)) + else: + mugrid.append((q, int(nf_default(q2, atlas)))) + + op_card.update({"mu0": theory["Q0"], "mugrid": mugrid}) + op_card["xgrid"] = x_grid # Specific defaults for evolven3fit evolution if theory["ModEv"] == "TRN": @@ -116,22 +146,3 @@ def construct_eko_cards( op_card = runcards.OperatorCard.from_dict(op_card) return theory_card, op_card - - -def split_evolgrid(evolgrid): - """Split the evolgrid in blocks according to the number of flavors and repeating the last entry of one block in the first entry of the next.""" - evolgrid_index_list = [] - evolgrid.sort() - starting_nf = evolgrid[0][1] - for evo_point in evolgrid: - current_nf = evo_point[1] - if current_nf != starting_nf: - evolgrid_index_list.append(evolgrid.index(evo_point)) - starting_nf = current_nf - start_index = 0 - evolgrid_list = [] - for index in evolgrid_index_list: - evolgrid_list.append(evolgrid[start_index : index + 1]) - start_index = index - evolgrid_list.append(evolgrid[start_index:]) - return evolgrid_list diff --git a/n3fit/src/evolven3fit_new/evolve.py b/n3fit/src/evolven3fit_new/evolve.py index cddb1eb308..8caab4e3e2 100644 --- a/n3fit/src/evolven3fit_new/evolve.py +++ b/n3fit/src/evolven3fit_new/evolve.py @@ -1,3 +1,4 @@ +from collections import defaultdict import logging import pathlib import sys @@ -145,7 +146,7 @@ def load_fit(usr_path): """ nnfitpath = usr_path / "nnfit" pdf_dict = {} - for yaml_file in nnfitpath.glob("replica_*/*.exportgrid"): + for yaml_file in nnfitpath.glob(f"replica_*/{usr_path.name}.exportgrid"): data = yaml.safe_load(yaml_file.read_text(encoding="UTF-8")) pdf_dict[yaml_file.parent.stem] = data return pdf_dict @@ -178,19 +179,27 @@ def evolve_exportgrid(exportgrid, eko, x_grid, qed): # generate block to dump targetgrid = eko.bases.targetgrid.tolist() - def ev_pdf(pid, x, Q2): - return x * evolved_pdf[Q2]["pdfs"][pid][targetgrid.index(x)] + # Finally separate by nf block (and order per nf/q) + by_nf = defaultdict(list) + for q, nf in sorted(eko.evolgrid, key=lambda ep: ep[1]): + by_nf[nf].append(q) + q2block_per_nf = {nf: sorted(qs) for nf, qs in by_nf.items()} - evolgrid_list = eko_utils.split_evolgrid(eko.evolgrid) blocks = [] - for evgrid in evolgrid_list: + for nf, q2grid in q2block_per_nf.items(): + + def pdf_xq2(pid, x, Q2): + x_idx = targetgrid.index(x) + return x * evolved_pdf[(Q2, nf)]["pdfs"][pid][x_idx] + block = genpdf.generate_block( - ev_pdf, + pdf_xq2, xgrid=targetgrid, - evolgrid=evgrid, + sorted_q2grid=q2grid, pids=basis_rotation.flavor_basis_pids, ) blocks.append(block) + return blocks diff --git a/n3fit/src/evolven3fit_new/utils.py b/n3fit/src/evolven3fit_new/utils.py index f1303fcbd7..4df9274459 100644 --- a/n3fit/src/evolven3fit_new/utils.py +++ b/n3fit/src/evolven3fit_new/utils.py @@ -7,7 +7,7 @@ from reportengine.compat import yaml from validphys.pdfbases import PIDS_DICT -DEFAULT_Q2GRID = ( +Q2GRID_Nf04 = ( np.array( [ 1.6500000e00, @@ -22,7 +22,6 @@ 3.8800751e00, 4.3584516e00, 4.9200000e00, - 4.9200000e00, 5.5493622e00, 6.2897452e00, 7.1650687e00, @@ -65,6 +64,62 @@ ** 2 ) +Q2GRID_Nf03 = ( + np.array( + [ + 1.0000000e00, + 1.0768843e00, + 1.1642787e00, + 1.2640247e00, + 1.3783565e00, + 1.5100000e00, + 1.6573843e00, + 1.8279487e00, + 2.0263188e00, + 2.2582323e00, + 2.5308507e00, + 2.8531703e00, + 3.2365690e00, + 3.6955380e00, + 4.2486693e00, + 4.9200000e00, + 5.6571821e00, + 6.5475141e00, + 7.6300446e00, + 8.9555329e00, + 1.0590474e01, + 1.2622686e01, + 1.5169120e01, + 1.8386905e01, + 2.2489085e01, + 2.7767274e01, + 3.4624624e01, + 4.3624282e01, + 5.5561424e01, + 7.1571582e01, + 9.3295496e01, + 1.2313315e02, + 1.6464038e02, + 2.2315640e02, + 3.0681103e02, + 4.2816505e02, + 6.0692308e02, + 8.7449251e02, + 1.2817733e03, + 1.9127020e03, + 2.9082314e03, + 4.5095982e03, + 7.1379509e03, + 1.1543948e04, + 1.9094934e04, + 3.2338760e04, + 5.6137084e04, + 1.0000000e05, + ] + ) + ** 2 +) + class LhapdfLike: """ @@ -123,7 +178,7 @@ def get_theoryID_from_runcard(usr_path): return my_runcard["theory"]["theoryid"] -def generate_q2grid(Q0, Qfin, Q_points, match_dict): +def generate_q2grid(Q0, Qfin, Q_points, match_dict, nf0=None): """Generate the q2grid used in the final evolved pdfs or use the default grid if Qfin or Q_points is not provided. @@ -131,7 +186,14 @@ def generate_q2grid(Q0, Qfin, Q_points, match_dict): in order to obtain the relative matching scale. """ if Qfin is None and Q_points is None: - return DEFAULT_Q2GRID + if nf0 == 4: + return Q2GRID_Nf04 + elif nf0 == 3: + return Q2GRID_Nf03 + elif nf0 is None: + raise ValueError("In order to use a default grid, a value of nf0 must be provided") + else: + raise NotImplementedError(f"No default grid in Q available for {nf0=}") elif Qfin is None or Q_points is None: raise ValueError("q_fin and q_points must be specified either both or none of them") else: @@ -148,7 +210,9 @@ def generate_q2grid(Q0, Qfin, Q_points, match_dict): frac_of_point = np.log(match_scale / Q_ini) / np.log(Qfin / Q0) num_points = int(Q_points * frac_of_point) num_points_list.append(num_points) - grids.append(np.geomspace(Q_ini**2, match_scale**2, num=num_points)) + grids.append( + np.geomspace(Q_ini**2, match_scale**2, num=num_points, endpoint=False) + ) Q_ini = match_scale num_points = Q_points - sum(num_points_list) grids.append(np.geomspace(Q_ini**2, Qfin**2, num=num_points)) 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 1796095df1..9dd372c60e 100644 --- a/n3fit/src/n3fit/tests/test_evolven3fit.py +++ b/n3fit/src/n3fit/tests/test_evolven3fit.py @@ -23,11 +23,6 @@ def assert_sorted(arr, title): raise ValueError(f"The values of {title} are not sorted!") -def check_consecutive_members(grid, value): - """Check if the first occurrence of value in grid is followed by value again""" - return np.allclose(grid[list(grid).index(value) + 1], value) - - def check_lhapdf_info(info_path): """Check the LHAPDF info file is correct""" info = yaml.load(info_path.open("r", encoding="utf-8")) @@ -77,24 +72,40 @@ def check_lhapdf_dat(dat_path, info): # Use allclose here to avoid failing because of a different in the 7th place np.testing.assert_allclose(q[-1], info["QMax"]) + +def test_generate_q2grid(): + """Tests the creation of the default grids for different values of nf + and whether the matched grid is generating points in the desired locations + """ + # nf 3, q0 = 1.0 + grid = utils.generate_q2grid(None, None, None, {}, 3) + assert grid[0] == 1.0**2 + # nf 4, q0 = 1.65 + grid = utils.generate_q2grid(None, None, None, {}, 4) + assert grid[0] == 1.65**2 + + for nf in [1, 2, 5, 6]: + with pytest.raises(NotImplementedError): + grid = utils.generate_q2grid(None, None, None, {}, nf) + + with pytest.raises(ValueError): + grid = utils.generate_q2grid(None, None, None, {}) -def test_utils(): - # Testing the default grid - grid = utils.generate_q2grid(1.65, None, None, {}) - assert_allclose(1.65**2, grid[0]) - assert len(grid) == 50 - # We expect the bottom mass to be repeated twice because it is intended once in 4 flavors and once in 5 flavors. - assert check_consecutive_members(grid, 4.92**2) - # Testing if the points of the matching are correctly repeated twice matched_grid = utils.generate_q2grid(1.65, 1.0e5, 100, {4.92: 2.0, 100: 1.0}) - assert len(matched_grid) == 100 + t1 = 4.92 * 2.0 + t2 = 100.0 * 1.0 + + assert_allclose((1.65) ** 2, matched_grid[0]) assert_allclose((1.0e5) ** 2, matched_grid[-1]) - assert check_consecutive_members(matched_grid, (4.92 * 2.0) ** 2) - assert check_consecutive_members(matched_grid, (100.0 * 1.0) ** 2) + assert t1**2 in matched_grid + assert t2**2 in matched_grid + + +def test_utils(): # Testing the fake LHAPDF class q20 = 1.65**2 x_grid = np.geomspace(1.0e-7, 1.0, 30) - fake_grids = [[x * (1.0 - x) for x in x_grid] for pid in PIDS_DICT.keys()] + fake_grids = [[x * (1.0 - x) for x in x_grid] for _ in PIDS_DICT.keys()] pdf_grid = dict([(pid, v) for pid, v in zip(range(len(PIDS_DICT)), fake_grids)]) my_PDF = utils.LhapdfLike(pdf_grid, q20, x_grid) assert my_PDF.hasFlavor(6) @@ -133,10 +144,15 @@ def test_eko_utils(tmp_path): t_card_dict["order"][0] == pto + 1 ) # This is due to a different convention in eko orders due to QED assert_allclose(op_card_dict["xgrid"], x_grid) - assert_allclose(op_card_dict["mugrid"][0], (1.65, 4)) + # In theory 162 the charm threshold is at 1.51 + # and we should find two entries, one for nf=3 and another one for nf=4 + assert_allclose(op_card_dict["mugrid"][0], (1.51, 3)) + assert_allclose(op_card_dict["mugrid"][1], (1.51, 4)) + # Then (with the number of points we chosen it will happen in position 2,3 + # we will find the bottom threshold at two different nf + assert_allclose(op_card_dict["mugrid"][2], (4.92, 4)) + assert_allclose(op_card_dict["mugrid"][3], (4.92, 5)) assert_allclose(op_card_dict["mugrid"][-1], (q_fin, 5)) - # In this case there are not enough points to have twice the bottom matching scale - assert_allclose(op_card_dict["mugrid"][1], (4.92, 5)) # Testing computation of eko save_path = tmp_path / "ekotest.tar" runner.solve(t_card, op_card, save_path) @@ -145,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) @@ -169,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/commondata/PLOTTING_ATLAS_WZ_TOT_13TEV.yaml b/nnpdfcpp/data/commondata/PLOTTING_ATLAS_WZ_TOT_13TEV.yaml index c05aa86e87..73bbec06cb 100644 --- a/nnpdfcpp/data/commondata/PLOTTING_ATLAS_WZ_TOT_13TEV.yaml +++ b/nnpdfcpp/data/commondata/PLOTTING_ATLAS_WZ_TOT_13TEV.yaml @@ -3,8 +3,7 @@ dataset_label: "ATLAS $W,Z$ inclusive 13 TeV" x: " " y_label: $\sigma^{fid}$ (fb) -figure_by: - - boson + extra_labels: " ": - $W^-$ diff --git a/nnpdfcpp/data/theory.db b/nnpdfcpp/data/theory.db index 1f953e2459..fa202758df 100644 Binary files a/nnpdfcpp/data/theory.db and b/nnpdfcpp/data/theory.db differ diff --git a/validphys2/examples/pdf_lumi_plots.yaml b/validphys2/examples/pdf_lumi_plots.yaml index 67d046dba7..c570249a83 100644 --- a/validphys2/examples/pdf_lumi_plots.yaml +++ b/validphys2/examples/pdf_lumi_plots.yaml @@ -13,7 +13,8 @@ pdf: {id: "NNPDF40_nlo_as_01180", label: "4.0 NLO"} sqrts: 13000 # GeV lumi_channel: "gg" # one of [gg, gq, qqbar, qq, ddbar, uubar, ssbar, - # ccbar, bbbar, dubar, udbar, scbar, csbar, pp, gp] + # ccbar, bbbar, dubar, udbar, scbar, csbar, pp, gp, + # wlum1, zlum1] PDFscalespecs: - xscale: log diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py index 8ae5fe1012..7a8088ed31 100644 --- a/validphys2/src/validphys/dataplots.py +++ b/validphys2/src/validphys/dataplots.py @@ -387,6 +387,7 @@ def _plot_fancy_impl( lb = labellist[normalize_to] ax.set_ylabel(f"Ratio to {lb if lb else norm_result.label}") + ax.legend().set_zorder(100000) ax.set_xlabel(info.xlabel) fig.tight_layout() @@ -879,7 +880,6 @@ def plot_smpdf(pdf, dataset, obs_pdf_correlations, mark_threshold: float = 0.9): info = get_info(dataset) table = kitable(dataset, info) - figby = sane_groupby_iter(table, info.figure_by) basis = obs_pdf_correlations.basis @@ -891,54 +891,83 @@ def plot_smpdf(pdf, dataset, obs_pdf_correlations, mark_threshold: float = 0.9): plotting_var = info.get_xcol(table) - # TODO: vmin vmax should be global or by figure? - vmin, vmax = min(plotting_var), max(plotting_var) - if info.x_scale == 'log': - norm = mcolors.LogNorm(vmin, vmax) + categorical = not np.issubdtype(plotting_var.dtype, np.number) + if categorical: + # Plot lines using a categorical color map (for a reasonable number of + # categories), and set up the categorical labels (used below). + categorical_keys, values = np.unique(plotting_var, return_inverse=True) + plotting_var = values + num_categories = len(categorical_keys) + if num_categories <= len(cm.Set2.colors): + cmap = mcolors.ListedColormap(cm.Set2.colors[:num_categories]) + else: + cmap = cm.viridis.resample(num_categories) + bins = np.linspace(0, num_categories, num_categories + 1) + norm = mcolors.BoundaryNorm(bins, num_categories) + else: - norm = mcolors.Normalize(vmin, vmax) - # http://stackoverflow.com/a/11558629/1007990 - sm = cm.ScalarMappable(cmap=cm.viridis, norm=norm) + cmap = cm.viridis + #TODO: vmin vmax should be global or by figure? + vmin, vmax = min(plotting_var), max(plotting_var) + if info.x_scale == 'log': + norm = mcolors.LogNorm(vmin, vmax) + else: + norm = mcolors.Normalize(vmin, vmax) + + + table["__plotting_var"] = plotting_var + sm = cm.ScalarMappable(cmap=cmap, norm=norm) + + figby = sane_groupby_iter(table, info.figure_by) for same_vals, fb in figby: - grid = fullgrid[np.asarray(fb.index), ...] + grid = fullgrid[ np.asarray(fb.index),...] - # Use the maximum absolute correlation for plotting purposes + + #Use the maximum absolute correlation for plotting purposes absgrid = np.max(np.abs(grid), axis=0) - mark_mask = absgrid > np.max(absgrid) * mark_threshold + mark_mask = absgrid > np.max(absgrid)*mark_threshold label = info.group_label(same_vals, info.figure_by) - # TODO: PY36ScalarMappable - # TODO Improve title? - title = "%s %s\n[%s]" % (info.dataset_label, '(%s)' % label if label else '', pdf.label) - - # Start plotting - w, h = mpl.rcParams["figure.figsize"] - h *= 2.5 - fig, axes = plotutils.subplots(nrows=nf, sharex=True, figsize=(w, h), sharey=True) + #TODO: PY36ScalarMappable + #TODO Improve title? + title = f"{info.dataset_label} {label if label else ''}\n[{pdf.label}]" + + #Start plotting + w,h = mpl.rcParams["figure.figsize"] + h*=2.5 + fig, axes = plotutils.subplots(nrows=nf, sharex=True, figsize=(w,h), sharey=True) fig.suptitle(title) - colors = sm.to_rgba(info.get_xcol(fb)) + colors = sm.to_rgba(fb["__plotting_var"]) for flindex, (ax, fl) in enumerate(zip(axes, fls)): - for i, color in enumerate(colors): - ax.plot(x, grid[i, flindex, :].T, color=color) + for i,color in enumerate(colors): + ax.plot(x, grid[i,flindex,:].T, color=color) + - flmask = mark_mask[flindex, :] + flmask = mark_mask[flindex,:] ranges = split_ranges(x, flmask, filter_falses=True) for r in ranges: ax.axvspan(r[0], r[-1], color='#eeeeff') - ax.set_ylabel("$%s$" % basis.elementlabel(fl)) + ax.set_ylabel("$%s$"%basis.elementlabel(fl)) ax.set_xscale(scale_from_grid(obs_pdf_correlations)) - ax.set_ylim(-1, 1) + ax.set_ylim(-1,1) ax.set_xlim(x[0], x[-1]) ax.set_xlabel('$x$') - # fig.subplots_adjust(hspace=0) - fig.colorbar(sm, ax=axes.ravel().tolist(), label=info.xlabel, aspect=100) - # TODO: Fix title for this - # fig.tight_layout() - yield fig + cbar = fig.colorbar( + sm, + ax=axes.ravel().tolist(), + label=info.xlabel, + aspect=100, + ) + if categorical: + cbar.set_ticks(np.linspace(0.5, num_categories - 0.5, num_categories)) + cbar.ax.set_yticklabels(categorical_keys) + #TODO: Fix title for this + #fig.tight_layout() + yield fig @figure def plot_obscorrs(corrpair_datasets, obs_obs_correlations, pdf): @@ -1225,14 +1254,14 @@ def plot_xq2( highlight_datasets = set() def next_options(): - # Get the colors + #Get the colors prop_settings = mpl.rcParams['axes.prop_cycle'] - # Apparently calling the object gives us an infinite cycler + #Apparently calling the object gives us an infinite cycler settings_cycler = prop_settings() - # So far, I don't understand how this is done with mpl "cycler" - # objects, or wether I like it. So far this is godd enough - for markeropts, settings in zip(plotutils.marker_iter_plot(), settings_cycler): - # Override last with first + #So far, I don't understand how this is done with mpl "cycler" + #objects, or wether I like it. So far this is godd enough + for markeropts, settings in zip(plotutils.marker_iter_plot(), settings_cycler): + #Override last with first options = { 'linestyle': 'none', **markeropts, @@ -1243,7 +1272,7 @@ def next_options(): next_opts = next_options() key_options = {} - for experiment, commondata, fitted, masked, group in dataset_inputs_by_groups_xq2map: + for (experiment, commondata, fitted, masked, group) in dataset_inputs_by_groups_xq2map: info = get_info(commondata) if marker_by == 'process type': key = info.process_description diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py index 91c80ce2fb..ad3e300070 100644 --- a/validphys2/src/validphys/filters.py +++ b/validphys2/src/validphys/filters.py @@ -110,7 +110,7 @@ def filter_closure_data(filter_path, data, fakepdf, fakenoise, filterseed, sep_m def filter_closure_data_by_experiment( - filter_path, experiments_data, fakepdf, fakenoise, filterseed, experiments_index, sep_mult + filter_path, experiments_data, fakepdf, fakenoise, filterseed, data_index, sep_mult ): """ Like :py:func:`filter_closure_data` except filters data by experiment. @@ -124,7 +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 = data_index[data_index.isin([exp.name], level=0)] res.append( _filter_closure_data( filter_path, exp, fakepdf, fakenoise, filterseed, experiment_index, sep_mult @@ -179,7 +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, data_index, sep_mult ): """ This function is accessed within a closure test only, that is, the fakedata @@ -213,7 +213,7 @@ def _filter_closure_data( random noise added to Level 0 data - experiments_index : pandas.MultiIndex + data_index : pandas.MultiIndex Returns @@ -241,7 +241,7 @@ def _filter_closure_data( if fakenoise: # ======= Level 1 closure test =======# - closure_data = make_level1_data(data, closure_data, filterseed, experiments_index, sep_mult) + closure_data = make_level1_data(data, closure_data, filterseed, data_index, sep_mult) # ====== write commondata and systype files ======# if fakenoise: diff --git a/validphys2/src/validphys/gridvalues.py b/validphys2/src/validphys/gridvalues.py index 83628a9997..96d2e3aacb 100644 --- a/validphys2/src/validphys/gridvalues.py +++ b/validphys2/src/validphys/gridvalues.py @@ -32,6 +32,8 @@ 'csbar': r'c\bar{s}', 'pp': r'\gamma\gamma', 'gp': r'g\gamma', + 'zlum1': r'u\bar{u} + d\bar{d}', + 'wlum1': r'u\bar{d} + d\bar{u}', } QUARK_COMBINATIONS = { @@ -115,6 +117,15 @@ def central_grid_values(pdf: PDF, flmat, xmat, qmat): # TODO: Investigate writting these in cython/cffi/numba/... +def _parton_pair_lumi_inner(pdf_set, n, mx, x1, x2, i, j): + """Helper to evaluate lumis for pairs of partons.""" + # fmt: off + return ( + pdf_set.xfxQ(x1, mx, n, i)*pdf_set.xfxQ(x2, mx, n, j) + + pdf_set.xfxQ(x1, mx, n, j)*pdf_set.xfxQ(x2, mx, n, i) + ) + + def evaluate_luminosity( pdf_set: LHAPDFSet, n: int, s: float, mx: float, x1: float, x2: float, channel ): @@ -155,11 +166,24 @@ def evaluate_luminosity( # as in the second of Eq.(4) in arXiv:1607.01831 res = sum(a*b for a,b in itertools.product(r1,r2)) + elif channel == 'zlum1': + u, ubar = 2, -2 + d, dbar = -1, 1 + res = ( + _parton_pair_lumi_inner(pdf_set=pdf_set, n=n, mx=mx, x1=x1, x2=x2, i=u, j=ubar) + + _parton_pair_lumi_inner(pdf_set=pdf_set, n=n, mx=mx, x1=x1, x2=x2, i=d, j=dbar) + ) + elif channel == 'wlum1': + u, dbar = 2, -1 + d, ubar = 1, -2 + res = ( + _parton_pair_lumi_inner(pdf_set=pdf_set, n=n, mx=mx, x1=x1, x2=x2, i=u, j=dbar) + + _parton_pair_lumi_inner(pdf_set=pdf_set, n=n, mx=mx, x1=x1, x2=x2, i=d, j=ubar) + ) + elif channel in QUARK_COMBINATIONS.keys(): i, j = QUARK_COMBINATIONS[channel] - res = (pdf_set.xfxQ(x1, mx, n, i) * pdf_set.xfxQ(x2, mx, n, j) - + pdf_set.xfxQ(x1, mx, n, j) * pdf_set.xfxQ(x2, mx, n, i)) - + res = _parton_pair_lumi_inner(pdf_set=pdf_set, n=n, mx=mx, x1=x1, x2=x2, i=i, j=j) else: raise ValueError("Bad channel") # fmt: on diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 4dde4d952c..c9a5805f07 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -17,10 +17,9 @@ log = logging.getLogger(__name__) -Q_IN = 100 +# 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 @@ -60,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 @@ -71,6 +79,8 @@ def __init__(self, theoryid, lux_params, replicas): 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" + with EKO.read(self.path_to_eko_photon) as eko: + self.q_in = np.sqrt(eko.mu20) # set fiatlux self.lux = {} @@ -110,7 +120,9 @@ 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): @@ -129,7 +141,9 @@ def compute_photon_array(self, replica): """ # 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 = np.array( + [self.lux[replica].EvaluatePhoton(x, self.q_in**2).total for x in XGRID] + ) photon_qin += self.generate_errors(replica) # fiatlux computes x * gamma(x) photon_qin /= XGRID @@ -141,8 +155,8 @@ def compute_photon_array(self, replica): # 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): + pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID))) + for j, pid in enumerate(eko.bases.inputpids): if pid == 22: pdfs_init[j] = photon_qin ph_id = j @@ -150,12 +164,11 @@ def compute_photon_array(self, replica): 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] + [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID] ) # Apply EKO to PDFs - q2 = eko.mu2grid[0] - with eko.operator(q2) as elem: + for _, elem in eko.items(): pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init) photon_Q0 = pdfs_final[ph_id] @@ -188,7 +201,7 @@ def error_matrix(self): if not self.additional_errors: return None extra_set = self.additional_errors.load() - qs = [Q_IN] * len(XGRID) + qs = [self.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): diff --git a/validphys2/src/validphys/pseudodata.py b/validphys2/src/validphys/pseudodata.py index f5d6fb0389..1fd0a5233d 100644 --- a/validphys2/src/validphys/pseudodata.py +++ b/validphys2/src/validphys/pseudodata.py @@ -293,7 +293,7 @@ def level0_commondata_wc(data, fakepdf): return level0_commondata_instances_wc -def make_level1_data(data, level0_commondata_wc, filterseed, experiments_index, sep_mult): +def make_level1_data(data, level0_commondata_wc, filterseed, data_index, sep_mult): """ Given a list of Level 0 commondata instances, return the same list with central values replaced by Level 1 data. @@ -331,7 +331,7 @@ def make_level1_data(data, level0_commondata_wc, filterseed, experiments_index, filterseed : int random seed used for the generation of Level 1 data - experiments_index : pandas.MultiIndex + data_index : pandas.MultiIndex Returns ------- @@ -367,16 +367,21 @@ def make_level1_data(data, level0_commondata_wc, filterseed, experiments_index, level0_commondata_wc, filterseed, covmat, sep_mult=sep_mult, genrep=True ) - indexed_level1_data = indexed_make_replica(experiments_index, level1_data) + indexed_level1_data = indexed_make_replica(data_index, level1_data) + + dataset_order = {cd.setname: i for i, cd in enumerate(level0_commondata_wc)} # ===== create commondata instances with central values given by pseudo_data =====# level1_commondata_dict = {c.setname: c for c in level0_commondata_wc} level1_commondata_instances_wc = [] + # note that groupby sorts alphabetically for xx, grp in indexed_level1_data.groupby('dataset'): level1_commondata_instances_wc.append( level1_commondata_dict[xx].with_central_value(grp.values) ) + # sort back so as to mantain same order as in level0_commondata_wc + level1_commondata_instances_wc.sort(key=lambda x: dataset_order[x.setname]) return level1_commondata_instances_wc diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index 393891d223..e5a23c9dd4 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -26,6 +26,7 @@ ) from validphys.convolution import PredictionsRequireCutsError, predictions from validphys.core import PDF, DataGroupSpec, DataSetSpec, Stats +from validphys.plotoptions import get_info log = logging.getLogger(__name__) @@ -151,6 +152,36 @@ def from_convolution(cls, pdf, posset): return cls(stats) +def data_index(data): + """ + Given a core.DataGroupSpec instance, return pd.MultiIndex + with the following levels: + + 1. experiment + 2. datasets + 3. datapoints indices (cuts already applied to) + + + Parameters + ---------- + data: core.DataGroupSpec + + Returns + ------- + pd.MultiIndex + + """ + tuples = [] + + for ds in data.datasets: + experiment = get_info(ds).experiment + + for i in ds.cuts.load(): + tp = (experiment, ds.name, i) + tuples.append(tp) + return pd.MultiIndex.from_tuples(tuples, names=('experiment', 'dataset', 'id')) + + # TODO: finish deprecating all dependencies on this index largely in theorycovmat module groups_data = collect("data", ("group_dataset_inputs_by_metadata",)) diff --git a/validphys2/src/validphys/tests/baseline/test_plot_smpdf.png b/validphys2/src/validphys/tests/baseline/test_plot_smpdf.png index 7274fa0469..883dc4d077 100644 Binary files a/validphys2/src/validphys/tests/baseline/test_plot_smpdf.png and b/validphys2/src/validphys/tests/baseline/test_plot_smpdf.png differ diff --git a/validphys2/src/validphys/tests/baseline/test_plot_smpdf_categorical.png b/validphys2/src/validphys/tests/baseline/test_plot_smpdf_categorical.png new file mode 100644 index 0000000000..42349bad95 Binary files /dev/null and b/validphys2/src/validphys/tests/baseline/test_plot_smpdf_categorical.png differ diff --git a/validphys2/src/validphys/tests/conftest.py b/validphys2/src/validphys/tests/conftest.py index 0f36cd2455..75d35c718c 100644 --- a/validphys2/src/validphys/tests/conftest.py +++ b/validphys2/src/validphys/tests/conftest.py @@ -27,6 +27,8 @@ def tmp(tmpdir): SINGLE_DATASET = {'dataset': 'NMC'} +SINGLE_CATEGORICAL = {"dataset": "ATLAS_WZ_TOT_13TEV", 'cfac': ["QCD"]} + DATA = [ {'dataset': 'NMC'}, {'dataset': 'ATLASTTBARTOT', 'cfac':['QCD']}, @@ -54,6 +56,7 @@ def tmp(tmpdir): HESSIAN_PDF = "NNPDF40_nnlo_as_01180_hessian" THEORYID = 162 THEORYID_NEW = 399 +THEORY_QED = 398 FIT = "NNPDF40_nnlo_lowprecision" FIT_3REPLICAS = "Basic_runcard_3replicas_lowprec_221130" FIT_3REPLICAS_DCUTS = "Basic_runcard_3replicas_diffcuts_230221" @@ -95,6 +98,16 @@ def single_data_internal_cuts_config(data_internal_cuts_config): config_dict.update(dataset_input=DATA[0]) return config_dict +@pytest.fixture(scope='module') +def single_data_categorical_internal_cuts_config(data_internal_cuts_config): + """Test dataset with categorical plotting variables""" + return { + **data_internal_cuts_config, + 'dataset_input': SINGLE_CATEGORICAL, + # NOTE: The old theory is currently bugged for this dataset + 'theoryid': THEORYID_NEW, + } + @pytest.fixture(scope='module') def single_data_single_point_internal_cuts_config(single_data_internal_cuts_config): config_dict = dict(single_data_internal_cuts_config) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index 30805be889..82787e7492 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -1,145 +1,86 @@ -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 -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] +import yaml - 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 +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 - def fxq(self): - return 0 +from ..conftest import PDF, THEORY_QED -class FakeF2LO: - def __init__(self, pdfs, theory): - self.pdfs = pdfs - self.theory = theory +def generate_fiatlux_runcard(): + return { + "luxset": PDFset(PDF), + # 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 + } - def fxq(self): - return 0 +def test_parameters_init(): + "test initailization of the parameters from Photon class" + fiatlux_runcard = generate_fiatlux_runcard() + test_theory = API.theoryid(theoryid=THEORY_QED) -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)) + # we are not testing the photon here so we make it faster + fiatlux_runcard['eps_base'] = 1e-1 - 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.additional_errors.name, "LUXqed17_plus_PDF4LHC15_nnlo_100") 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.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" + test_theory = API.theoryid(theoryid=THEORY_QED) + 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" + test_theory = API.theoryid(theoryid=THEORY_QED) + 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" + test_theory = API.theoryid(theoryid=THEORY_QED) + alpha = Alpha(test_theory.get_description()) vec_beta0 = [ -0.5305164769729844, -0.6719875374991137, @@ -155,3 +96,75 @@ 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 = generate_fiatlux_runcard() + fiatlux_runcard["additional_errors"] = False + test_theory = API.theoryid(theoryid=THEORY_QED) + 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..5e83e6458f 100644 --- a/validphys2/src/validphys/tests/photon/test_structurefunctions.py +++ b/validphys2/src/validphys/tests/photon/test_structurefunctions.py @@ -1,7 +1,11 @@ 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, THEORY_QED class ZeroPdfs: @@ -9,35 +13,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 +39,72 @@ 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() + test_theory = API.theoryid(theoryid=THEORY_QED) + 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 + test_theory = API.theoryid(theoryid=THEORY_QED) + 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() + test_theory = API.theoryid(theoryid=THEORY_QED) + 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 + ) diff --git a/validphys2/src/validphys/tests/test_plots.py b/validphys2/src/validphys/tests/test_plots.py index 004514c71b..da85cc5cd5 100644 --- a/validphys2/src/validphys/tests/test_plots.py +++ b/validphys2/src/validphys/tests/test_plots.py @@ -40,6 +40,11 @@ def test_dataspecschi2(): def test_plot_smpdf(single_data_internal_cuts_config): return next(API.plot_smpdf(**single_data_internal_cuts_config)) +@pytest.mark.linux +@pytest.mark.mpl_image_compare +def test_plot_smpdf_categorical(single_data_categorical_internal_cuts_config): + return next(API.plot_smpdf(**single_data_categorical_internal_cuts_config)) + @pytest.mark.linux @pytest.mark.mpl_image_compare def test_plot_obscorrs(single_data_internal_cuts_config):