diff --git a/n3fit/runcards/examples/Basic_runcard_pc_covmat.yml b/n3fit/runcards/examples/Basic_runcard_pc_covmat.yml new file mode 100644 index 0000000000..991e5a7cba --- /dev/null +++ b/n3fit/runcards/examples/Basic_runcard_pc_covmat.yml @@ -0,0 +1,151 @@ +# +# Configuration file for n3fit +# +###################################################################################### +description: NNPDF4.0 ht with TCM - DIS (NC & CC) only + +###################################################################################### +dataset_inputs: +- {dataset: NMC_NC_NOTFIXED_EM-F2, frac: 0.75, variant: legacy_dw} +- {dataset: NMC_NC_NOTFIXED_P_EM-SIGMARED, frac: 0.75, variant: legacy} +- {dataset: SLAC_NC_NOTFIXED_D_EM-F2, frac: 0.75, variant: legacy_dw} +- {dataset: BCDMS_NC_NOTFIXED_P_EM-F2, frac: 0.75, variant: legacy_dw} +- {dataset: CHORUS_CC_NOTFIXED_PB_NU-SIGMARED, frac: 0.75, variant: legacy_dw} +- {dataset: NUTEV_CC_NOTFIXED_FE_NB-SIGMARED, cfac: [MAS], frac: 0.75, variant: legacy_dw} +- {dataset: HERA_CC_318GEV_EP-SIGMARED, frac: 0.75, variant: legacy} +- {dataset: HERA_NC_318GEV_EAVG_CHARM-SIGMARED, frac: 0.75, variant: legacy} +- {dataset: HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED, frac: 0.75, variant: legacy} +- {dataset: DYE866_Z0_800GEV_DW_RATIO_PDXSECRATIO, frac: 0.75, variant: legacy} +- {dataset: CDF_Z0_1P96TEV_ZRAP, frac: 0.75, variant: legacy} +- {dataset: ATLAS_Z0J_8TEV_PT-Y, frac: 0.75, variant: legacy_10} +- {dataset: ATLAS_1JET_8TEV_R06_PTY, frac: 0.75, variant: legacy_decorrelated} +- {dataset: ATLAS_2JET_7TEV_R06_M12Y, frac: 0.75, variant: legacy} +- {dataset: CMS_2JET_7TEV_M12Y, frac: 0.75} +- {dataset: CMS_1JET_8TEV_PTY, frac: 0.75, variant: legacy} +- {dataset: LHCB_Z0_13TEV_DIELECTRON-Y, frac: 0.75} + +################################################################################ +datacuts: + t0pdfset: 240701-02-rs-nnpdf40-baseline + q2min: 2.5 + w2min: 3.24 + +################################################################################ +# NNLO QCD TRN evolution +theory: + theoryid: 708 + +theorycovmatconfig: + point_prescriptions: ["9 point", "power corrections"] + pc_parameters: + H2p: {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]} + H2d: {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]} + HLp: {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]} + HLd: {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]} + H3p: {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]} + H3d: {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]} + Hj: {yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]} + H2j_ATLAS: {yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]} + H2j_CMS: {yshift: [2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25]} + pc_included_procs: ["JETS", "DIJET", "DIS NC", "DIS CC"] + pc_excluded_exps: [HERA_NC_318GEV_EAVG_CHARM-SIGMARED, + HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED,] + pdf: 210619-n3fit-001 + use_thcovmat_in_fitting: true + use_thcovmat_in_sampling: true +resample_negative_pseudodata: false + +# For fits <= 4.0 multiplicative and additive uncertainties were sampled separately +# and thus the flag `separate_multiplicative` needs to be set to True +# sampling: +# separate_multiplicative: True + +################################################################################ +trvlseed: 591866982 +nnseed: 945709987 +mcseed: 519562661 +genrep: true + +################################################################################ +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: 3000 + 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 + savepseudodata: True + basis: + - {fl: sng, trainable: false, smallx: [1.089, 1.119], largex: [1.475, 3.119]} + - {fl: g, trainable: false, smallx: [0.7504, 1.098], largex: [2.814, 5.669]} + - {fl: v, trainable: false, smallx: [0.479, 0.7384], largex: [1.549, 3.532]} + - {fl: v3, trainable: false, smallx: [0.1073, 0.4397], largex: [1.733, 3.458]} + - {fl: v8, trainable: false, smallx: [0.5507, 0.7837], largex: [1.516, 3.356]} + - {fl: t3, trainable: false, smallx: [-0.4506, 0.9305], largex: [1.745, 3.424]} + - {fl: t8, trainable: false, smallx: [0.5877, 0.8687], largex: [1.522, 3.515]} + - {fl: t15, trainable: false, smallx: [1.089, 1.141], largex: [1.492, 3.222]} + +################################################################################ +positivity: + posdatasets: + # Positivity Lagrange Multiplier + - {dataset: NNPDF_POS_2P24GEV_F2U, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_F2D, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_F2S, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_FLL, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_DYU, maxlambda: 1e10} + - {dataset: NNPDF_POS_2P24GEV_DYD, maxlambda: 1e10} + - {dataset: NNPDF_POS_2P24GEV_DYS, maxlambda: 1e10} + - {dataset: NNPDF_POS_2P24GEV_F2C, maxlambda: 1e6} + # Positivity of MSbar PDFs + - {dataset: NNPDF_POS_2P24GEV_XUQ, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_XUB, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_XDQ, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_XDB, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_XSQ, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_XSB, maxlambda: 1e6} + - {dataset: NNPDF_POS_2P24GEV_XGL, maxlambda: 1e6} + +added_filter_rules: + - dataset: NNPDF_POS_2P24GEV_FLL + rule: "x > 5.0e-7" + - dataset: NNPDF_POS_2P24GEV_F2C + rule: "x < 0.74" + - dataset: NNPDF_POS_2P24GEV_XGL + rule: "x > 0.1" + - dataset: NNPDF_POS_2P24GEV_XUQ + rule: "x > 0.1" + - dataset: NNPDF_POS_2P24GEV_XUB + rule: "x > 0.1" + - dataset: NNPDF_POS_2P24GEV_XDQ + rule: "x > 0.1" + - dataset: NNPDF_POS_2P24GEV_XDB + rule: "x > 0.1" + - dataset: NNPDF_POS_2P24GEV_XSQ + rule: "x > 0.1" + - dataset: NNPDF_POS_2P24GEV_XSB + rule: "x > 0.1" + +integrability: + integdatasets: + - {dataset: NNPDF_INTEG_3GEV_XT8, maxlambda: 1e2} + - {dataset: NNPDF_INTEG_3GEV_XT3, maxlambda: 1e2} + +################################################################################ +debug: false +maxcores: 8 diff --git a/n3fit/src/n3fit/performfit.py b/n3fit/src/n3fit/performfit.py index e0fd1af9b5..8cdc8fd10f 100644 --- a/n3fit/src/n3fit/performfit.py +++ b/n3fit/src/n3fit/performfit.py @@ -41,7 +41,7 @@ def performfit( debug=False, maxcores=None, double_precision=False, - parallel_models=True, + parallel_models=False, ): """ This action will (upon having read a validcard) process a full PDF fit diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py index d63429cc70..3da207cf0e 100644 --- a/n3fit/src/n3fit/scripts/vp_setupfit.py +++ b/n3fit/src/n3fit/scripts/vp_setupfit.py @@ -1,8 +1,8 @@ #!/usr/bin/env python """ - setup-fit - prepare and apply data cuts before fit - setup-fit constructs the fit [results] folder where data used by nnfit - will be stored. +setup-fit - prepare and apply data cuts before fit +setup-fit constructs the fit [results] folder where data used by nnfit +will be stored. """ # Implementation notes diff --git a/validphys2/src/validphys/checks.py b/validphys2/src/validphys/checks.py index 4f958b0dc4..d0d608c7b4 100644 --- a/validphys2/src/validphys/checks.py +++ b/validphys2/src/validphys/checks.py @@ -361,3 +361,15 @@ def check_darwin_single_process(NPROC): """ if platform.system() == "Darwin" and NPROC != 1: raise CheckError("NPROC must be set to 1 on OSX, because multithreading is not supported.") + + +@make_argcheck +def check_pc_parameters(pc_parameters): + """Check that the parameters for the PC method are set correctly""" + for par in pc_parameters.values(): + # Check that the length of shifts is the same as the length of nodes + if (len(par['yshift']) != len(par['nodes'])): + raise ValueError( + f"The length of nodes does not match that of the list in {par['ht']}." + f"Check the runcard. Got {len(par['yshift'])} != {len(par['nodes'])}" + ) diff --git a/validphys2/src/validphys/config.py b/validphys2/src/validphys/config.py index 04f306a858..146846b337 100644 --- a/validphys2/src/validphys/config.py +++ b/validphys2/src/validphys/config.py @@ -1276,6 +1276,21 @@ def produce_nnfit_theory_covmat( f = user_covmat_fitting return f + + def produce_mult_factor_user_covmat(self, mult_factor: float = None, user_covmat_path: str = None): + """ + Multiplicative factors for the user covmat provided by mult_factors_user_covmat in the runcard. + If no factors are provided, returns None. + For use in theorycovariance.construction.user_covmat. + """ + # Check that if mult_factors is provided, user_covmat_paths is also provided + if mult_factor is not None and user_covmat_path is None: + raise ConfigError("If mult_factors is provided, user_covmat_paths must also be provided.") + + if mult_factor is None: + return 1.0 if user_covmat_path is not None else None + else: + return mult_factor def produce_fitthcovmat( self, use_thcovmat_if_present: bool = False, fit: (str, type(None)) = None @@ -1799,6 +1814,8 @@ def produce_theoryids(self, t0id, point_prescription): prescription. The options for the latter are defined in pointprescriptions.yaml. This hard codes the theories needed for each prescription to avoid user error.""" th = t0id.id + if point_prescription == 'power corrections': + return NSList([t0id], nskey="theoryid") lsv = yaml_safe.load(read_text(validphys.scalevariations, "scalevariationtheoryids.yaml")) @@ -1881,7 +1898,18 @@ def produce_filter_data( if not fakedata: return validphys.filters.filter_real_data else: - if inconsistent_fakedata: + # TODO we don't want to sample from the theory covmat for L1 data, + # but we do want to use the theory covmat for L2 data + if theorycovmatconfig is not None and theorycovmatconfig.get( + "use_thcovmat_in_fakedata_sampling" + ): + # NOTE: By the time we run theory covmat closure tests, + # hopefully the generation of pseudodata will be done in python. + raise ConfigError( + "Generating L1 closure test data which samples from the theory " + "covariance matrix has not been implemented yet." + ) + elif inconsistent_fakedata: log.info("Using filter for inconsistent closure data") return validphys.filters.filter_inconsistent_closure_data_by_experiment @@ -1909,6 +1937,19 @@ def produce_total_phi_data(self, fitthcovmat): return validphys.results.total_phi_data_from_experiments return validphys.results.dataset_inputs_phi_data + @configparser.explicit_node + def produce_covs_pt_prescrip(self, point_prescription): + if point_prescription != 'power corrections': + from validphys.theorycovariance.construction import covs_pt_prescrip_mhou + + f = covs_pt_prescrip_mhou + else: + from validphys.theorycovariance.construction import covs_pt_prescrip_pc + + f = covs_pt_prescrip_pc + + return f + class Config(report.Config, CoreConfig): """The effective configuration parser class.""" diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py index 961c272406..2829fa2f58 100644 --- a/validphys2/src/validphys/dataplots.py +++ b/validphys2/src/validphys/dataplots.py @@ -2,8 +2,6 @@ Plots of relations between data PDFs and fits. """ -from __future__ import generator_stop - from collections import defaultdict from collections.abc import Sequence import itertools @@ -28,7 +26,7 @@ from validphys.core import CutsPolicy, MCStats, cut_mask from validphys.plotoptions.core import get_info, kitable, transform_result from validphys.results import chi2_stat_labels, chi2_stats -from validphys.sumrules import POL_LIMS, partial_polarized_sum_rules +from validphys.sumrules import POL_LIMS from validphys.utils import sane_groupby_iter, scale_from_grid, split_ranges log = logging.getLogger(__name__) @@ -301,9 +299,7 @@ def _plot_fancy_impl( min_vals = [] max_vals = [] fig, ax = plotutils.subplots() - ax.set_title( - "{} {}".format(info.dataset_label, info.group_label(samefig_vals, info.figure_by)) - ) + ax.set_title(f"{info.dataset_label} {info.group_label(samefig_vals, info.figure_by)}") lineby = sane_groupby_iter(fig_data, info.line_by) @@ -1287,7 +1283,7 @@ def _check_display_cuts_requires_use_cuts(display_cuts, use_cuts): @make_argcheck def _check_marker_by(marker_by): - markers = ('process type', 'experiment', 'dataset', 'group') + markers = ('process type', 'experiment', 'dataset', 'group', 'kinematics') if marker_by not in markers: raise CheckError("Unknown marker_by value", marker_by, markers) @@ -1346,7 +1342,8 @@ def plot_xq2( will be displaed and marked. The points are grouped according to the `marker_by` option. The possible - values are: "process type", "experiment", "group" or "dataset". + values are: "process type", "experiment", "group" or "dataset" for discrete + colors, or "kinematics" for coloring by 1/(Q2(1-x)) Some datasets can be made to appear highlighted in the figure: Define a key called ``highlight_datasets`` containing the names of the datasets to be @@ -1477,6 +1474,7 @@ def plot_xq2( xh = defaultdict(list) q2h = defaultdict(list) + cvdict = defaultdict(list) if not highlight_datasets: highlight_datasets = set() @@ -1507,6 +1505,8 @@ def next_options(): elif marker_by == "group": # if group is None then make sure that shows on legend. key = str(group) + elif marker_by == "kinematics": + key = None else: raise ValueError('Unknown marker_by value') @@ -1522,6 +1522,7 @@ def next_options(): xdict = x q2dict = q2 + cvdict[key].append(commondata.load().get_cv()) xdict[key].append(fitted[0]) q2dict[key].append(fitted[1]) if display_cuts: @@ -1536,6 +1537,13 @@ def next_options(): else: # This is to get the label key coords = [], [] + if marker_by == "kinematics": + ht_magnitude = np.concatenate(cvdict[key]) / (coords[1] * (1 - coords[0])) + out = ax.scatter( + *coords, marker='.', c=ht_magnitude, cmap="viridis", norm=mcolors.LogNorm() + ) + clb = fig.colorbar(out) + clb.ax.set_title(r'$F_\mathrm{exp}\frac{1}{Q^2(1-x)}$') ax.plot(*coords, label=key, markeredgewidth=1, markeredgecolor=None, **key_options[key]) # Iterate again so highlights are printed on top. diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py index add8a3bf4b..f02d843cb9 100644 --- a/validphys2/src/validphys/filters.py +++ b/validphys2/src/validphys/filters.py @@ -125,7 +125,6 @@ def to_dict(self): class FilterRule: """ Dataclass which carries the filter rule information. - """ dataset: str = None @@ -172,6 +171,9 @@ def default_filter_rules_input(): are unique, i.d. if there are no multiple rules for the same dataset of process with the same rule (`reason` and `local_variables` are not hashed). """ + # TODO: This should be done using a more sophisticated comparison + # that checks if two rules are actually the same, regardless of the + # order in which the cuts are defined. list_rules = yaml_safe.load(read_text(validphys.cuts, "filters.yaml")) unique_rules = set(FilterRule(**rule) for rule in list_rules) if len(unique_rules) != len(list_rules): diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index bfc3279782..596488f90e 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -239,7 +239,7 @@ def data_index(data): experiments_data = collect("data", ("group_dataset_inputs_by_experiment",)) -procs_data = collect("data", ("group_dataset_inputs_by_process",)) +groups_data_by_process = collect("data", ("group_dataset_inputs_by_process",)) def groups_index(groups_data, diagonal_basis=False): @@ -279,12 +279,12 @@ def groups_index(groups_data, diagonal_basis=False): return df.index -def experiments_index(experiments_data, diagonal_basis=False): - return groups_index(experiments_data, diagonal_basis) +def experiments_index(experiments_data): + return groups_index(experiments_data) -def procs_index(procs_data): - return groups_index(procs_data) +def procs_index(groups_data_by_process): + return groups_index(groups_data_by_process) def groups_data_values(group_result_table): @@ -817,10 +817,12 @@ def groups_chi2_table(groups_data, pdf, groups_chi2, groups_each_dataset_chi2): @table -def procs_chi2_table(procs_data, pdf, groups_chi2_by_process, groups_each_dataset_chi2_by_process): +def procs_chi2_table( + groups_data_by_process, pdf, groups_chi2_by_process, groups_each_dataset_chi2_by_process +): """Same as groups_chi2_table but by process""" return groups_chi2_table( - procs_data, pdf, groups_chi2_by_process, groups_each_dataset_chi2_by_process + groups_data_by_process, pdf, groups_chi2_by_process, groups_each_dataset_chi2_by_process ) diff --git a/validphys2/src/validphys/theorycovariance/construction.py b/validphys2/src/validphys/theorycovariance/construction.py index e4811c4478..e2d05513be 100644 --- a/validphys2/src/validphys/theorycovariance/construction.py +++ b/validphys2/src/validphys/theorycovariance/construction.py @@ -4,6 +4,7 @@ """ from collections import defaultdict, namedtuple +import dataclasses import logging import numpy as np @@ -13,7 +14,10 @@ from reportengine.table import table pass +from validphys.checks import check_pc_parameters +from validphys.core import PDF from validphys.results import results, results_central +from validphys.theorycovariance.higher_twist_functions import compute_deltas_pc from validphys.theorycovariance.theorycovarianceutils import ( check_correct_theory_combination, check_fit_dataset_order_matches_grouped, @@ -47,11 +51,18 @@ def theory_covmat_dataset(results, results_central_bytheoryids, point_prescripti return thcovmat -ProcessInfo = namedtuple("ProcessInfo", ("preds", "namelist", "sizes")) +@dataclasses.dataclass(frozen=True) +class ProcessInfo: + """Dataclass containing the information needed to construct the theory covariance matrix.""" + preds: dict + namelist: dict + sizes: dict + data_spec: dict -def combine_by_type(each_dataset_results_central_bytheory): - """Groups the datasets bu process and returns an instance of the ProcessInfo class + +def combine_by_type(each_dataset_results_central_bytheory, groups_data_by_process): + """Groups the datasets by process and returns an instance of the ProcessInfo class Parameters ---------- @@ -68,6 +79,7 @@ def combine_by_type(each_dataset_results_central_bytheory): dataset_size = defaultdict(list) theories_by_process = defaultdict(list) ordered_names = defaultdict(list) + data_spec = defaultdict(list) for dataset in each_dataset_results_central_bytheory: name = dataset[0][0].name theory_centrals = [x[1].central_value for x in dataset] @@ -77,8 +89,14 @@ def combine_by_type(each_dataset_results_central_bytheory): theories_by_process[proc_type].append(theory_centrals) for key, item in theories_by_process.items(): theories_by_process[key] = np.concatenate(item, axis=1) + + # Store DataGroupSpecs instances + for group_proc in groups_data_by_process: + for exp_set in group_proc.datasets: + data_spec[exp_set.name] = exp_set + process_info = ProcessInfo( - preds=theories_by_process, namelist=ordered_names, sizes=dataset_size + preds=theories_by_process, namelist=ordered_names, sizes=dataset_size, data_spec=data_spec ) return process_info @@ -206,6 +224,40 @@ def covmat_n3lo_ad(name1, name2, deltas1, deltas2): return 1 / norm * s +def covmat_power_corrections(deltas1, deltas2): + """Returns the the theory covariance sub-matrix for power + corrections. The two arguments `deltas1` and `deltas2` contain + the shifts for the firs and second experiment, respectively. + + The shifts are given in this form: + ``` + deltas1 = {shift1_label: array1_of_shifts1, + shift2_label: array1_of_shifts2, + shift3_label: array1_of_shifts3, + ...} + deltas2 = {shift1_label: array2_of_shifts1, + shift2_label: array2_of_shifts2, + shift3_label: array2_of_shifts3, + ...} + ``` + The sub-matrix is computed using the 5-point prescription, thus + + s = array1_of_shifts1 X array2_of_shifts1 + array1_of_shifts2 X array2_of_shifts2 + ... + + where `X` is the outer product. + """ + # Check that `deltas1` and `deltas2` have the same shifts + if deltas1.keys() != deltas2.keys(): + raise RuntimeError('The two dictionaries do not contain the same shifts.') + + size1 = next(iter(deltas1.values())).size + size2 = next(iter(deltas2.values())).size + s = np.zeros(shape=(size1, size2)) + for shift in deltas1.keys(): + s += np.outer(deltas1[shift], deltas2[shift]) + return s + + def compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2=None, deltas2=None): """Utility to compute the covariance matrix by prescription given the shifts with respect to the central value for a pair of processes. @@ -276,28 +328,30 @@ def compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2=None, del # alphas is correlated for all datapoints and the covmat construction is # therefore equivalent to that of the factorization scale variations s = covmat_3fpt(deltas1, deltas2) + elif point_prescription == 'power corrections': + # Shifts computed from power corrected predictions + s = covmat_power_corrections(deltas1, deltas2) return s @check_correct_theory_combination -def covs_pt_prescrip(combine_by_type, point_prescription): +def covs_pt_prescrip_mhou(combine_by_type, point_prescription): """Produces the sub-matrices of the theory covariance matrix according to a point prescription which matches the number of input theories. - If 5 theories are provided, a scheme 'bar' or 'nobar' must be chosen in the runcard in order to specify the prescription. Sub-matrices correspond to applying the scale variation prescription to each pair of processes in turn, using a different procedure for the case where the processes are the same relative to when they are different.""" - process_info = combine_by_type running_index = 0 + + covmats = defaultdict(list) start_proc = defaultdict(list) for name in process_info.preds: size = len(process_info.preds[name][0]) start_proc[name] = running_index running_index += size - covmats = defaultdict(list) for name1 in process_info.preds: for name2 in process_info.preds: central1, *others1 = process_info.preds[name1] @@ -307,6 +361,49 @@ def covs_pt_prescrip(combine_by_type, point_prescription): s = compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2, deltas2) start_locs = (start_proc[name1], start_proc[name2]) covmats[start_locs] = s + + return covmats + + +@check_pc_parameters +def covs_pt_prescrip_pc( + combine_by_type, + point_prescription, + pdf: PDF, + pc_parameters, + pc_included_procs, + pc_excluded_exps, +): + """Produces the sub-matrices of the theory covariance matrix for power + corrections. Sub-matrices correspond to applying power corrected shifts + to each pair of `datasets`.""" + process_info = combine_by_type + datagroup_spec = process_info.data_spec + running_index = 0 + + covmats = defaultdict(list) + start_proc_by_exp = defaultdict(list) + for exp_name, data_spec in datagroup_spec.items(): + start_proc_by_exp[exp_name] = running_index + running_index += data_spec.load_commondata().ndata + + for exp_name1, data_spec1 in datagroup_spec.items(): + for exp_name2, data_spec2 in datagroup_spec.items(): + process_type1 = process_lookup(exp_name1) + process_type2 = process_lookup(exp_name2) + + is_excluded_exp = any(name in pc_excluded_exps for name in [exp_name1, exp_name2]) + is_included_proc = any( + proc not in pc_included_procs for proc in [process_type1, process_type2] + ) + if not (is_excluded_exp or is_included_proc): + deltas1 = compute_deltas_pc(data_spec1, pdf, pc_parameters) + deltas2 = compute_deltas_pc(data_spec2, pdf, pc_parameters) + s = compute_covs_pt_prescrip( + point_prescription, exp_name1, deltas1, exp_name2, deltas2 + ) + start_locs = (start_proc_by_exp[exp_name1], start_proc_by_exp[exp_name2]) + covmats[start_locs] = s return covmats @@ -329,7 +426,7 @@ def theory_covmat_custom_per_prescription(covs_pt_prescrip, procs_index, combine covmat_index = pd.MultiIndex.from_tuples(indexlist, names=procs_index.names) # Put the covariance matrices between two process into a single covariance matrix - total_datapoints = sum(combine_by_type.sizes.values()) + total_datapoints = sum(process_info.sizes.values()) mat = np.zeros((total_datapoints, total_datapoints), dtype=np.float32) for locs, cov in covs_pt_prescrip.items(): xsize, ysize = cov.shape @@ -339,7 +436,7 @@ def theory_covmat_custom_per_prescription(covs_pt_prescrip, procs_index, combine @table -def fromfile_covmat(covmatpath, procs_data, procs_index): +def fromfile_covmat(covmatpath, groups_data_by_process, procs_index): """Reads a general theory covariance matrix from file. Then 1: Applies cuts to match experiment covariance matrix 2: Expands dimensions to match experiment covariance matrix @@ -353,7 +450,7 @@ def fromfile_covmat(covmatpath, procs_data, procs_index): # Reordering covmat to match exp order in runcard # Datasets in exp covmat dslist = [] - for group in procs_data: + for group in groups_data_by_process: for ds in group.datasets: dslist.append(ds.name) # Datasets in filecovmat in exp covmat order @@ -368,7 +465,7 @@ def fromfile_covmat(covmatpath, procs_data, procs_index): # ------------- # # Loading cuts to apply to covariance matrix indextuples = [] - for group in procs_data: + for group in groups_data_by_process: for ds in group.datasets: # Load cuts for each dataset in the covmat if ds.name in filecovmat.index.get_level_values(1): @@ -434,7 +531,7 @@ def fromfile_covmat(covmatpath, procs_data, procs_index): @table -def user_covmat(procs_data, procs_index, loaded_user_covmat_path): +def user_covmat(groups_data_by_process, procs_index, loaded_user_covmat_path, mult_factor_user_covmat): """ General theory covariance matrix provided by the user. Useful for testing the impact of externally produced @@ -444,7 +541,7 @@ def user_covmat(procs_data, procs_index, loaded_user_covmat_path): ``user_covmat_path`` in ``theorycovmatconfig`` in the runcard. For more information see documentation. """ - return fromfile_covmat(loaded_user_covmat_path, procs_data, procs_index) + return mult_factor_user_covmat * fromfile_covmat(loaded_user_covmat_path, groups_data_by_process, procs_index) @table diff --git a/validphys2/src/validphys/theorycovariance/higher_twist_functions.py b/validphys2/src/validphys/theorycovariance/higher_twist_functions.py new file mode 100644 index 0000000000..11876c3cee --- /dev/null +++ b/validphys2/src/validphys/theorycovariance/higher_twist_functions.py @@ -0,0 +1,479 @@ +""" +This module contains the utilities for the computation of the shifts in +the theoretical predictions due to the effects of power corrections. Contrary +to scale variations, in the case of power corrections the shifts are not +computed using theories. The shifts are computed at "runtime" during vp-setupfit. + +The aim is that, after shifts being computed, the covmat can be constructed using +the same functions implemented for scale variations (e.g. `covs_pt_prescrip`). +The way shifts are computed depends also on the point prescription. In the case of +scale variations, the point prescription specifies the theories whereby shifts are +computed. In the case of power corrections, shifts and covmat constructions are +computed using a 5-point prescription extended to every parameter used to define +the power correction. + +This module comprehends a bunch of ``factory`` functions such as `mult_dis_pc`. Each +of these functions returns another function that computes the shifts taking as arguments +the values of the parameters used to parametrise the power corrections. In +other words, these factory functions hard-code the dependence on the kinematic +and leave the dependence on the parameters free. In this way, the shifts +can be computed using different combinations of parameters (i.e. different prescriptions) +if needed. +""" + +from collections import defaultdict +import operator + +import numpy as np +import numpy.typing as npt + +from validphys.convolution import central_fk_predictions +from validphys.core import PDF, DataSetSpec + +F2P_exps = ['SLAC_NC_NOTFIXED_P_EM-F2', 'BCDMS_NC_NOTFIXED_P_EM-F2'] +F2D_exps = ['SLAC_NC_NOTFIXED_D_EM-F2', 'BCDMS_NC_NOTFIXED_D_EM-F2'] +NC_SIGMARED_P_EM = ['NMC_NC_NOTFIXED_P_EM-SIGMARED', 'HERA_NC_318GEV_EM-SIGMARED'] +NC_SIGMARED_P_EP = [ + 'HERA_NC_225GEV_EP-SIGMARED', + 'HERA_NC_251GEV_EP-SIGMARED', + 'HERA_NC_300GEV_EP-SIGMARED', + 'HERA_NC_318GEV_EP-SIGMARED', +] +NC_SIGMARED_P_EAVG = ['HERA_NC_318GEV_EAVG_CHARM-SIGMARED', 'HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED'] + + +def linear_bin_function( + a: npt.ArrayLike, y_shift: npt.ArrayLike, bin_edges: npt.ArrayLike +) -> np.ndarray: + """ + This function defines the linear bin function used to construct the prior. Specifically, + the prior is constructed using a triangular function whose value at the peak of the node + is linked to the right and left nodes using a straight line. + + Parameters + ---------- + a: ArrayLike of float + A one-dimensional array of points at which the function is evaluated. + y_shift: ArrayLike of float + A one-dimensional array whose elements represent the y-value of each bin + bin_nodes: ArrayLike of float + A one-dimensional array containing the edges of the bins. The bins are + constructed using pairs of consecutive points. + + Return + ------ + A one-dimensional array containing the function values evaluated at the points + specified in `a`. + """ + res = np.zeros_like(a) + for shift_pos, shift in enumerate(y_shift): + if shift_pos > 0 and shift_pos < len(y_shift) - 1: + bin_low = bin_edges[shift_pos - 1] + bin_high = bin_edges[shift_pos + 1] + bin_mid = bin_edges[shift_pos] + m1 = shift / (bin_mid - bin_low) + m2 = shift / (bin_high - bin_mid) + elif shift_pos == 0: # Left-most bin + bin_high = bin_edges[shift_pos + 1] + bin_mid = bin_edges[shift_pos] + bin_low = bin_mid + m1 = 0.0 + m2 = shift / (bin_high - bin_mid) + else: # Right-most bin + bin_low = bin_edges[shift_pos - 1] + bin_mid = bin_edges[shift_pos] + bin_high = bin_mid + m1 = shift / (bin_mid - bin_low) + m2 = 0.0 + cond_low = np.multiply( + a >= bin_low, a < bin_mid if shift_pos != len(y_shift) - 1 else a <= bin_mid + ) + cond_high = np.multiply( + a >= bin_mid, a < bin_high if shift_pos != len(y_shift) - 1 else a <= bin_high + ) + res = np.add(res, [m1 * (val - bin_low) if cond else 0.0 for val, cond in zip(a, cond_low)]) + res = np.add( + res, [-m2 * (val - bin_high) if cond else 0.0 for val, cond in zip(a, cond_high)] + ) + return res + + +def dis_pc_func( + delta_h: npt.ArrayLike, + nodes: npt.ArrayLike, + x: npt.ArrayLike, + Q2: npt.ArrayLike, +) -> npt.ArrayLike: + """ + This function builds the functional form of the power corrections for DIS-like processes. + Power corrections are modelled using a linear function, which interpolates between the nodes + of the parameterisation. The y-values for each node are given + by the array `delta_h`. The power corrections will be computed for the pairs (xb, Q2), + where `xb` is the Bjorken x. The power correction for DIS processes is divided by Q2. + + Parameters + ---------- + delta_h: ArrayLike + One-dimensional array containing the shifts for each bin. + nodes: ArrayLike + One-dimensional array containing the edges of the bins in x-Bjorken. + x: ArrayLike + List of x-Bjorken points at which the power correction function is evaluated. + Q2: ArrayLike + List of scales where the power correction function is evaluated. + + Returns + ------- + A one-dimensional array of power corrections for DIS-like processes where each point is + evaluated at the kinematic pair (x,Q2). + """ + PC = linear_bin_function(x, delta_h, nodes) / Q2 + return PC + + +def jets_pc_func( + delta_h: npt.ArrayLike, + nodes: npt.ArrayLike, + pT: npt.ArrayLike, + rap: npt.ArrayLike, +) -> npt.ArrayLike: + """ + Same as `dis_pc_func`, but for jet data. Here, the kinematic pair consists of the rapidity + `rap` and the transverse momentum `pT`. + + Parameters + ---------- + delta_h: ArrayLike + One-dimensional array containing the shifts for each bin. + nodes: ArrayLike + One-dimensional array containing the edges of the bins in rapidity. + rap: ArrayLike + List of rapidity points at which the power correction is evaluated. + pT: ArrayLike + List of pT points at which the power correction is evaluated. + + Returns + ------- + A one-dimensional array of power corrections for jet processes where each point is + evaluated at the kinematic pair (y, pT). + """ + PC = linear_bin_function(rap, delta_h, nodes) / pT + return PC + + +# def jet_single_par(delta_h: float, pT: npt.ArrayLike, rap: npt.ArrayLike) -> npt.ArrayLike: +# ret = [delta_h for _ in range(rap.size)] +# return np.array(ret) / pT + +# def mult_jet_pc_single_par(dataset_sp, pdf, pc_nodes, pT, rap, pc_func_type: str = "step"): +# """ +# As mult_jet_pc, but with one single shift for all rapidity bins. + +# This function is meant to be for development purposes only. It will either substitute +# mult_jet_pc or be deleted in the future.""" +# cuts = dataset_sp.cuts +# (fkspec,) = dataset_sp.fkspecs +# fk = fkspec.load_with_cuts(cuts) +# xsec = central_fk_predictions(fk, pdf) + +# def func(y_values): +# assert y_values.size == 1 +# ret = [y_values[0] for _ in range(rap.size)] +# ret = np.array(ret) / pT +# return np.multiply(ret, xsec.to_numpy()[:, 0]) + +# return func + + +def mult_dis_pc(nodes, x, q2, dataset_sp, pdf): + """ + Returns the function that computes the shift to observables due to + power corrections. Power corrections are treated as multiplicative + shifts. Hence if `O` is the observable, the prediction is shifted as + + O -> O * (1 + PC), + + and the shift is defined as + + Delta(O) = O * (1 + PC) - O = O * PC. + + This function returns a function that computes the shift + given the y-values of the nodes used to define the power corrections. + """ + cuts = dataset_sp.cuts + (fkspec,) = dataset_sp.fkspecs + fk = fkspec.load_with_cuts(cuts) + th_preds = central_fk_predictions(fk, pdf) + + def func(y_values): + result = dis_pc_func(y_values, nodes, x, q2) + return np.multiply(result, th_preds.to_numpy()[:, 0]) + + return func + + +def mult_dis_ratio_pc(p_nodes, d_nodes, x, q2, dataset_sp, pdf): + """ + Returns the function that computes the shift for the ratio of structure + functions F2_d / F2_p. For this observable, power corrections are defined + such that + + F2_d / F2_p -> F2_d * (1 + PC2_d) / F2_p * (1 + PC2_p) , + + and the shift is the defined as + + Delta(F2 ratio) = F2_d / F2_p * (PC2_d - PC2_p) / (1 + PC2_d). + + As for `mult_dis_pc`, this function returns a function that computes the shift + for the ratio of structure functions F2_d / F2_p given a set of y-values. + + Parameters + ---------- + p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for the proton (see `dis_pc_func`). + d_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for the deuteron (see `dis_pc_func`). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + dataset_sp: DataSetSpec + An instance of DataSetSpec used to extract information such as cuts + and fk tables. + pdf: PDF + An instance of the class PDF. This specifies the PDF to bo convoluted + with the FK tables. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2_d and P2_p. + """ + cuts = dataset_sp.cuts + fkspec_F2D, fkspec_F2P = dataset_sp.fkspecs + fk_F2D = fkspec_F2D.load_with_cuts(cuts) + fk_F2P = fkspec_F2P.load_with_cuts(cuts) + F2D = central_fk_predictions(fk_F2D, pdf) + F2P = central_fk_predictions(fk_F2P, pdf) + + F2D = np.concatenate(F2D.values) + F2P = np.concatenate(F2P.values) + F2_ratio = operator.truediv(F2D, F2P) + + def func(y_values_p, y_values_d): + h2d = dis_pc_func(y_values_d, d_nodes, x, q2) + h2p = dis_pc_func(y_values_p, p_nodes, x, q2) + num = np.sum([h2d, -h2p], axis=0) + denom = np.sum([np.ones_like(h2p), h2p], axis=0) + result = np.multiply(operator.truediv(num, denom), F2_ratio) + return result + + return func + +def mult_jet_pc(nodes, pT, rap, dataset_sp, pdf): + """ + As `mult_dis_pc`, but for jet data. The power corrections are defined as + + xsec -> xsec * ( 1 + PC ), + + and the shift is defined as + + Delta(xsec) = (xsec + xsec) - xsec = PC. + + The power correction is a function of the rapidity. + """ + cuts = dataset_sp.cuts + (fkspec,) = dataset_sp.fkspecs + fk = fkspec.load_with_cuts(cuts) + xsec = central_fk_predictions(fk, pdf) + + def func(y_values): + result = jets_pc_func(y_values, nodes, pT, rap) + return np.multiply(result, xsec.to_numpy()[:, 0]) + + return func + +def construct_pars_combs(parameters_dict): + """Construct the combination of parameters (the ones that parametrize the power + corrections) used to compute the shifts. + + Example + ------- + Given the following dictionary that specifies that power corrections + ``` + pc_dict = { + 'H1': {'list': [1,2], 'nodes': [0,1]} } + 'H2': {'list': [3,4,5], 'nodes': [0,1,2]} } + } + ``` + then this functions constructs a list as follows + ``` + pars_combs = [ + {'label': 'H1(1)', 'comb': {'H1': [1,0], 'H2': [0,0,0]}, + {'label': 'H1(2)', 'comb': {'H1': [0,1], 'H2': [0,0,0]}, + {'label': 'H2(1)', 'comb': {'H1': [0,0], 'H2': [3,0,0]}, + {'label': 'H2(2)', 'comb': {'H1': [0,0], 'H2': [0,4,0]}, + {'label': 'H2(3)', 'comb': {'H1': [0,0], 'H2': [0,0,5]}, + ] + ``` + + Parameters + ---------- + """ + combinations = [] + for key, values in parameters_dict.items(): + for i in range(len(values['yshift'])): + # Create a new dictionary with all keys and zeroed-out values + new_dict = {k: np.zeros_like(v['yshift']) for k, v in parameters_dict.items()} + # Set the specific value for the current index + label = key + f'({i})' + new_dict[key][i] = values['yshift'][i] + new_dict = {'label': label, 'comb': new_dict} + combinations.append(new_dict) + + return combinations + +def compute_deltas_pc(dataset_sp: DataSetSpec, pdf: PDF, pc_dict: dict): + """ + Computes the shifts due to power corrections for a single dataset given + the set of parameters that model the power corrections. The result is + a dictionary containing as many arrays of shifts as the number of + combinations of the parameters. For instance, the final dictionary + may look like: + ``` + deltas1 = {comb1_label: array_of_shifts1, + comb2_label: array_of_shifts2, + comb3_label: array_of_shifts3, + ...} + ``` + Note that, as of now, we don't need to specify different prescriptions. + For that reason, the prescription adopted to construct the shifts is hard + coded in the function `construct_pars_combs`, and the prescription used to + compute the sub-matrix is hard-coded in `covmat_power_corrections`. + """ + + exp_name = dataset_sp.name + process_type = dataset_sp.commondata.metadata.process_type.name + cuts = dataset_sp.cuts.load() + + pars_combs = construct_pars_combs(pc_dict) + deltas = defaultdict(list) + + pc_func = None + if process_type.startswith('DIS'): + f2_p_nodes = pc_dict["f2p"]['nodes'] + f2_d_nodes = pc_dict["f2d"]['nodes'] + dis_cc_nodes = pc_dict["dis_cc"]['nodes'] + + x = dataset_sp.commondata.metadata.load_kinematics()['x'].to_numpy().reshape(-1)[cuts] + q2 = dataset_sp.commondata.metadata.load_kinematics()['Q2'].to_numpy().reshape(-1)[cuts] + + # F2 ratio + if exp_name == "NMC_NC_NOTFIXED_EM-F2": + pc_func_ratio = mult_dis_ratio_pc(f2_p_nodes, f2_d_nodes, x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func_ratio(pars_pc['comb']['f2p'], pars_pc['comb']['f2d']) + + # F2 proton traget + elif exp_name in F2P_exps: + pc_func = mult_dis_pc(f2_p_nodes, x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['f2p']) + + # F2 deuteron traget + elif exp_name in F2D_exps: + pc_func = mult_dis_pc(f2_d_nodes, x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['f2d']) + + # EMC + elif exp_name.startswith('EMC_NC_250GEV'): + raise NotImplementedError( + f"The {process_type} observable for {exp_name} " + "has not been implemented." + ) + + # HERA NC xsec + elif exp_name in np.concatenate([NC_SIGMARED_P_EM, NC_SIGMARED_P_EP, NC_SIGMARED_P_EAVG]): + pc_func = mult_dis_pc(f2_p_nodes, x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['f2p']) + + # CHORUS + elif exp_name.startswith('CHORUS_CC'): + pc_func = mult_dis_pc(dis_cc_nodes, x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['dis_cc']) + + # NuTeV + elif exp_name.startswith('NUTEV_CC'): + pc_func = mult_dis_pc(dis_cc_nodes , x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['dis_cc']) + + # HERA_CC + elif exp_name.startswith('HERA_CC'): + pc_func = mult_dis_pc(dis_cc_nodes, x, q2, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['dis_cc']) + + else: + raise ValueError( + f"The {process_type} observable for {exp_name} " "has not been implemented." + ) + + elif process_type == 'JET': + pc_jet_nodes = pc_dict["Hj"]['nodes'] + eta = dataset_sp.commondata.metadata.load_kinematics()['y'].to_numpy().reshape(-1)[cuts] + pT = dataset_sp.commondata.metadata.load_kinematics()['pT'].to_numpy().reshape(-1)[cuts] + + pc_func = mult_jet_pc(pc_jet_nodes, pT, eta, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['Hj']) + + elif process_type == 'DIJET': + + + if dataset_sp.commondata.metadata.experiment == 'ATLAS': + pc_jet_nodes = pc_dict["H2j_ATLAS"]['nodes'] if pc_dict.get("H2j_ATLAS") else pc_dict["H2j"]['nodes'] + eta_star = ( + dataset_sp.commondata.metadata.load_kinematics()['ystar'] + .to_numpy() + .reshape(-1)[cuts] + ) + m_jj = ( + dataset_sp.commondata.metadata.load_kinematics()['m_jj'] + .to_numpy() + .reshape(-1)[cuts] + ) + pc_func = mult_jet_pc(pc_jet_nodes, m_jj, eta_star, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['H2j_ATLAS'] if pc_dict.get("H2j_ATLAS") else pars_pc['comb']['H2j']) + + elif dataset_sp.commondata.metadata.experiment == 'CMS': + pc_jet_nodes = pc_dict["H2j_CMS"]['nodes'] if pc_dict.get("H2j_CMS") else pc_dict["H2j"]['nodes'] + eta_diff = ( + dataset_sp.commondata.metadata.load_kinematics()['ydiff'] + .to_numpy() + .reshape(-1)[cuts] + ) + m_jj = ( + dataset_sp.commondata.metadata.load_kinematics()['m_jj'] + .to_numpy() + .reshape(-1)[cuts] + ) + pc_func = mult_jet_pc(pc_jet_nodes, m_jj, eta_diff, dataset_sp, pdf) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['H2j_CMS'] if pc_dict.get("H2j_CMS") else pars_pc['comb']['H2j']) + + else: + raise ValueError( + f"{dataset_sp.commondata.metadata.experiment} is not implemented for DIJET." + ) + + else: + raise RuntimeError(f"{process_type} has not been implemented.") + + return deltas