diff --git a/validphys2/examples/data_theory_comparison_w_sv.yaml b/validphys2/examples/data_theory_comparison_w_sv.yaml new file mode 100644 index 0000000000..9894be049d --- /dev/null +++ b/validphys2/examples/data_theory_comparison_w_sv.yaml @@ -0,0 +1,90 @@ +meta: + title: "Example of comparison of data and theory including scale variations" + keywords: [data theory comparison, scale variations, example runcard] + author: "Validphys examples" + +use_cuts: "internal" + +theoryid: 717 # define the central theory +point_prescription: "3 point" +theoryids: + from_: scale_variation_theories +use_theorycovmat: true + +use_pdferr: true + +disable_pdferr: + use_pdferr: false + +dataset_inputs: + - { dataset: ATLASWZRAP36PB } + - { dataset: BCDMSP } + +dataspecs: + - pdf: NNPDF40_nnlo_lowprecision + speclabel: "pdf1" + + - pdf: Basic_runcard_3replicas_lowprec_221130 + speclabel: "pdf2" + +# for the chi2 comparison +pdfs: + - Basic_runcard_3replicas_lowprec_221130 + - NNPDF40_nnlo_lowprecision + +datanorm: + normalize_to: data + +template_text: | + Plots with scale variations + =========================== + + {@with matched_datasets_from_dataspecs@} + + {@dataset_name@} + ---------------- + {@ with disable_pdferr @} + {@ plot_fancy_sv_dataspecs @} + {@ datanorm plot_fancy_sv_dataspecs @} + {@ endwith @} + + ### chi2 distributions per PDF + {@ with dataspecs @} + {@ plot_chi2dist_sv @} + {@ endwith @} + + {@ endwith @} + + Plots without scale variations + ============================== + + {@with matched_datasets_from_dataspecs@} + + {@dataset_name@} + ---------------- + {@ with disable_pdferr @} + {@ plot_fancy_dataspecs @} + {@ datanorm plot_fancy_dataspecs @} + {@ endwith @} + + ### chi2 distributions per PDF + {@ with dataspecs @} + {@ plot_chi2dist @} + {@ endwith @} + + {@ endwith @} + + Chi2 Report per dataset + ======================= + + Without scale variations + ------------------------ + {@plot_datasets_pdfs_chi2@} + + With scale variations + --------------------- + {@plot_datasets_pdfs_chi2_sv@} + + +actions_: + - report(main=true) diff --git a/validphys2/src/validphys/checks.py b/validphys2/src/validphys/checks.py index 4ac9113f1d..1e4a33f103 100644 --- a/validphys2/src/validphys/checks.py +++ b/validphys2/src/validphys/checks.py @@ -43,6 +43,12 @@ def check_pdf_is_montecarlo_or_hessian(pdf, **kwargs): ) +@make_argcheck +def check_using_theory_covmat(use_theorycovmat): + """Check that the `use_theorycovmat` is set to True""" + check(use_theorycovmat, "Expecting `use_theorycovmat: true`") + + @make_argcheck def check_not_using_pdferr(use_pdferr=False, **kwargs): check( diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py index d9572cb203..09246b0f06 100644 --- a/validphys2/src/validphys/dataplots.py +++ b/validphys2/src/validphys/dataplots.py @@ -14,7 +14,6 @@ from matplotlib import colors as mcolors from matplotlib import ticker as mticker import numpy as np -import numpy.linalg as la import pandas as pd import scipy.stats as stats @@ -26,8 +25,8 @@ from validphys.checks import check_not_using_pdferr from validphys.core import CutsPolicy, MCStats, cut_mask from validphys.coredata import KIN_NAMES -from validphys.plotoptions import get_info, kitable, transform_result -from validphys.results import chi2_stat_labels +from validphys.plotoptions.core import get_info, kitable, transform_result +from validphys.results import chi2_stat_labels, chi2_stats from validphys.utils import sane_groupby_iter, scale_from_grid, split_ranges log = logging.getLogger(__name__) @@ -59,6 +58,13 @@ def plot_chi2dist(dataset, abs_chi2_data, chi2_stats, pdf): return fig +@figure +def plot_chi2dist_sv(dataset, abs_chi2_data_thcovmat, pdf): + """Same as ``plot_chi2dist`` considering also the theory covmat in the calculation""" + chi2_stats_thcovmat = chi2_stats(abs_chi2_data_thcovmat) + return plot_chi2dist(dataset, abs_chi2_data_thcovmat, chi2_stats_thcovmat, pdf) + + def _chi2_distribution_plots(chi2_data, stats, pdf, plot_type): fig, ax = plotutils.subplots() label = pdf.name @@ -421,7 +427,6 @@ def plot_fancy( See docs/plotting_format.md for details on the format of the PLOTTING files. """ - yield from _plot_fancy_impl( results=one_or_more_results, commondata=commondata, @@ -516,6 +521,32 @@ def plot_fancy_dataspecs( ) +@_check_same_dataset_name +@_check_dataspec_normalize_to +@figuregen +def plot_fancy_sv_dataspecs( + dataspecs_results_with_scale_variations, + dataspecs_commondata, + dataspecs_cuts, + dataspecs_speclabel, + normalize_to: (str, int, type(None)) = None, +): + """ + Exactly the same as ``plot_fancy_dataspecs`` but the theoretical results passed down + are modified so that the 1-sigma error bands correspond to a combination of the + PDF error and the scale variations collected over theoryids + + See: :py:func:`validphys.results.results_with_scale_variations` + """ + return plot_fancy_dataspecs( + dataspecs_results_with_scale_variations, + dataspecs_commondata, + dataspecs_cuts, + dataspecs_speclabel, + normalize_to=normalize_to, + ) + + def _scatter_marked(ax, x, y, marked_dict, *args, **kwargs): kwargs['s'] = kwargs.get('s', 30) + 10 x = np.array(x, copy=False) @@ -662,6 +693,16 @@ def plot_datasets_pdfs_chi2(data, each_dataset_chi2_pdfs, pdfs): return fig +each_dataset_chi2_sv = collect("abs_chi2_data_thcovmat", ("data",)) +each_dataset_chi2_pdfs_sv = collect("each_dataset_chi2_sv", ("pdfs",)) + + +@figure +def plot_datasets_pdfs_chi2_sv(data, each_dataset_chi2_pdfs_sv, pdfs): + """Same as ``plot_datasets_pdfs_chi2_sv`` with the chi²s computed including scale variations""" + return plot_datasets_pdfs_chi2(data, each_dataset_chi2_pdfs_sv, pdfs) + + @figure def plot_datasets_chi2_spider(groups_data, groups_chi2): """Plot the chi² of all datasets with bars.""" diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index b19eeff510..2fa314829d 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -19,10 +19,10 @@ from reportengine.table import table from validphys.calcutils import all_chi2, bootstrap_values, calc_chi2, calc_phi, central_chi2 from validphys.checks import ( - check_cuts_considered, check_pdf_is_montecarlo, check_speclabels_different, check_two_dataspecs, + check_using_theory_covmat, ) from validphys.convolution import PredictionsRequireCutsError, predictions from validphys.core import PDF, DataGroupSpec, DataSetSpec, Stats @@ -100,17 +100,21 @@ def sqrtcovmat(self): class ThPredictionsResult(StatsResult): - """Class holding theory prediction, inherits from StatsResult""" + """Class holding theory prediction, inherits from StatsResult + When created with `from_convolution`, it keeps tracks of the PDF for which it was computed + """ - def __init__(self, dataobj, stats_class, label=None): + def __init__(self, dataobj, stats_class, label=None, pdf=None, theoryid=None): self.stats_class = stats_class self.label = label statsobj = stats_class(dataobj.T) + self._pdf = pdf + self._theoryid = theoryid super().__init__(statsobj) @staticmethod def make_label(pdf, dataset): - """Deduce a reasonsble label for the result based on pdf and dataspec""" + """Deduce a reasonable label for the result based on pdf and dataspec""" th = dataset.thspec if hasattr(pdf, "label"): if hasattr(th, "label"): @@ -140,8 +144,44 @@ def from_convolution(cls, pdf, dataset): ) from e label = cls.make_label(pdf, dataset) + thid = dataset.thspec.id + + return cls(th_predictions, pdf.stats_class, label, pdf=pdf, theoryid=thid) + + +class ThUncertaintiesResult(StatsResult): + """Class holding central theory predictions and the error bar corresponding to + the theory uncertainties considered. + The error members of this class correspond to central +- error_bar + """ + + def __init__(self, central, std_err, label=None): + # All should be (ndata, 1) + self._central = central + self._std_err = std_err + self.stats = None + self.label = label + + @property + def rawdata(self): + return self._central + + @property + def error_members(self): + upper = self._central + self._std_err + lower = self._central - self._std_err + return np.concatenate([self._central, lower, upper], axis=1) + + @property + def central_value(self): + return self._central + + @property + def std_error(self): + return self._std_err - return cls(th_predictions, pdf.stats_class, label) + def __len__(self): + return self._central.shape[0] class PositivityResult(StatsResult): @@ -156,11 +196,11 @@ 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 ---------- @@ -209,13 +249,7 @@ def groups_index(groups_data): data_id = np.arange(dataset.commondata.ndata, dtype=int) for idat in data_id: records.append( - dict( - [ - ("group", str(group.name)), - ("dataset", str(dataset.name)), - ("id", idat), - ] - ) + dict([("group", str(group.name)), ("dataset", str(dataset.name)), ("id", idat)]) ) columns = ["group", "dataset", "id"] @@ -262,13 +296,7 @@ def group_result_table_no_table(groups_results, groups_index): ) result_records.append( - dict( - [ - ("data_central", dt_central), - ("theory_central", th_central), - *replicas, - ] - ) + dict([("data_central", dt_central), ("theory_central", th_central), *replicas]) ) if not result_records: log.warning("Empty records for group results") @@ -469,12 +497,55 @@ def results(dataset: (DataSetSpec), pdf: PDF, covariance_matrix, sqrt_covmat): The theory is specified as part of the dataset (a remnant of the old C++ layout) A group of datasets is also allowed. """ + # TODO: is the message about the usage of the theory covariance matrix here true? + # probably not in most cases... return ( DataResult(dataset, covariance_matrix, sqrt_covmat), ThPredictionsResult.from_convolution(pdf, dataset), ) +def results_with_theory_covmat(dataset, results, theory_covmat_dataset): + """Returns results with a modfy ``DataResult`` such that the covariance matrix includes + also the theory covmat. + This can be used to make use of results that consider scale variations without including + the theory covmat as part of the covariance matrix used by other validphys function. + Most notably, this can be used to compute the chi2 including theory errors while plotting + data theory covariance in which the experimental uncertainties are not stained by the thcovmat + """ + # TODO: in principle this function could be removed, and `results` could automagically + # include the theory covmat when `use_theorycovmat: true` by changing the nodes in `config.py` + # however at the moment config.py _loads_ theory covmats and we need to compute it on the fly + from .covmats import sqrt_covmat + + data_result, central_th_result = results + total_covmat = theory_covmat_dataset + data_result.covmat + + data_result = DataResult(dataset, total_covmat, sqrt_covmat(total_covmat)) + return (data_result, central_th_result) + + +def results_with_scale_variations(results, theory_covmat_dataset): + """Use the theory covariance matrix to generate a ThPredictionsResult-compatible object + modified so that its uncertainties correspond to a combination of + the PDF and theory (scale variations) errors added in quadrature. + This allows to plot results including scale variations + + By doing this we lose all information about prediction for the individual replicas or theories + """ + data_result, central_th_result = results + + # Use the central value and PDF error from the central theory + cv = central_th_result.central_value + pdf_error = central_th_result.std_error + sv_error_sq = np.diag(theory_covmat_dataset) + + total_error = np.sqrt(pdf_error**2 + sv_error_sq) + + theory_error_result = ThUncertaintiesResult(cv, total_error, label=central_th_result.label) + return (data_result, theory_error_result) + + def dataset_inputs_results( data, pdf: PDF, dataset_inputs_covariance_matrix, dataset_inputs_sqrt_covmat ): @@ -486,10 +557,7 @@ def dataset_inputs_results( # ``results`` to support this. # TODO: The above comment doesn't make sense after adding T0. Deprecate this def pdf_results( - dataset: (DataSetSpec, DataGroupSpec), - pdfs: Sequence, - covariance_matrix, - sqrt_covmat, + dataset: (DataSetSpec, DataGroupSpec), pdfs: Sequence, covariance_matrix, sqrt_covmat ): """Return a list of results, the first for the data and the rest for each of the PDFs.""" @@ -512,11 +580,9 @@ def one_or_more_results( and the next is either the prediction for pdf or for each of the pdfs. Which of the two is selected intelligently depending on the namespace, when executing as an action.""" - if pdf: + if pdf is not None: return results(dataset, pdf, covariance_matrix, sqrt_covmat) - else: - return pdf_results(dataset, pdfs, covariance_matrix, sqrt_covmat) - raise ValueError("Either 'pdf' or 'pdfs' is required") + return pdf_results(dataset, pdfs, covariance_matrix, sqrt_covmat) Chi2Data = namedtuple("Chi2Data", ("replica_result", "central_result", "ndata")) @@ -534,6 +600,11 @@ def abs_chi2_data(results): return Chi2Data(th_result.stats_class(chi2s[:, np.newaxis]), central_result, len(data_result)) +def abs_chi2_data_thcovmat(results_with_theory_covmat): + """The same as ``abs_chi2_data`` but considering as well the theory uncertainties""" + return abs_chi2_data(results_with_theory_covmat) + + def dataset_inputs_abs_chi2_data(dataset_inputs_results): """Like `abs_chi2_data` but for a group of inputs""" return abs_chi2_data(dataset_inputs_results) @@ -595,10 +666,7 @@ def dataset_inputs_bootstrap_phi_data(dataset_inputs_results, bootstrap_samples= dt, th = dataset_inputs_results diff = np.array(th.error_members - dt.central_value[:, np.newaxis]) phi_resample = bootstrap_values( - diff, - bootstrap_samples, - apply_func=(lambda x, y: calc_phi(y, x)), - args=[dt.sqrtcovmat], + diff, bootstrap_samples, apply_func=(lambda x, y: calc_phi(y, x)), args=[dt.sqrtcovmat] ) return phi_resample @@ -616,11 +684,7 @@ def dataset_inputs_bootstrap_chi2_central( diff = np.array(th.error_members - dt.central_value[:, np.newaxis]) cchi2 = lambda x, y: calc_chi2(y, x.mean(axis=1)) chi2_central_resample = bootstrap_values( - diff, - bootstrap_samples, - boot_seed=boot_seed, - apply_func=(cchi2), - args=[dt.sqrtcovmat], + diff, bootstrap_samples, boot_seed=boot_seed, apply_func=(cchi2), args=[dt.sqrtcovmat] ) return chi2_central_resample @@ -666,10 +730,7 @@ def groups_chi2_table(groups_data, pdf, groups_chi2, groups_each_dataset_chi2): def procs_chi2_table(procs_data, 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, + procs_data, pdf, groups_chi2_by_process, groups_each_dataset_chi2_by_process ) @@ -764,23 +825,14 @@ def dataset_chi2_table(chi2_stats, dataset): procs_chi2 = collect("dataset_inputs_abs_chi2_data", ("group_dataset_inputs_by_process",)) fits_groups_chi2_data = collect("groups_chi2", ("fits", "fitcontext")) -fits_groups = collect( - "groups_data", - ( - "fits", - "fitcontext", - ), -) +fits_groups = collect("groups_data", ("fits", "fitcontext")) # TODO: Possibly get rid of the per_point_data parameter and have separate # actions for absolute and relative tables. @table def fits_groups_chi2_table( - fits_name_with_covmat_label, - fits_groups, - fits_groups_chi2_data, - per_point_data: bool = True, + fits_name_with_covmat_label, fits_groups, fits_groups_chi2_data, per_point_data: bool = True ): """A table with the chi2 computed with the theory corresponding to each fit for all datasets in the fit, grouped according to a key in the metadata, the @@ -840,10 +892,7 @@ def fits_groups_phi_table(fits_name_with_covmat_label, fits_groups, fits_groups_ @table @check_speclabels_different def dataspecs_groups_chi2_table( - dataspecs_speclabel, - dataspecs_groups, - dataspecs_groups_chi2_data, - per_point_data: bool = True, + dataspecs_speclabel, dataspecs_groups, dataspecs_groups_chi2_data, per_point_data: bool = True ): """Same as fits_groups_chi2_table but for an arbitrary list of dataspecs.""" return fits_groups_chi2_table( @@ -862,10 +911,7 @@ def dataspecs_groups_chi2_table( @table def fits_datasets_chi2_table( - fits_name_with_covmat_label, - fits_groups, - fits_datasets_chi2_data, - per_point_data: bool = True, + fits_name_with_covmat_label, fits_groups, fits_datasets_chi2_data, per_point_data: bool = True ): """A table with the chi2 for each included dataset in the fits, computed with the theory corresponding to the fit. The result are indexed in two @@ -911,10 +957,7 @@ def fits_datasets_chi2_table( @table @check_speclabels_different def dataspecs_datasets_chi2_table( - dataspecs_speclabel, - dataspecs_groups, - dataspecs_datasets_chi2_data, - per_point_data: bool = True, + dataspecs_speclabel, dataspecs_groups, dataspecs_datasets_chi2_data, per_point_data: bool = True ): """Same as fits_datasets_chi2_table but for arbitrary dataspecs.""" return fits_datasets_chi2_table( @@ -932,10 +975,7 @@ def dataspecs_datasets_chi2_table( # TODO: Decide what to do with the horrible totals code. @table def fits_chi2_table( - fits_total_chi2_data, - fits_datasets_chi2_table, - fits_groups_chi2_table, - show_total: bool = False, + fits_total_chi2_data, fits_datasets_chi2_table, fits_groups_chi2_table, show_total: bool = False ): """Show the chi² of each and number of points of each dataset and experiment of each fit, where experiment is a group of datasets according to the `experiment` key in @@ -1016,7 +1056,7 @@ def total_chi2_data_from_experiments(experiments_chi2_data, pdf): """ central_result = np.sum( - [exp_chi2_data.central_result for exp_chi2_data in experiments_chi2_data], + [exp_chi2_data.central_result for exp_chi2_data in experiments_chi2_data] ) # we sum data, not error_members here because we feed it back into the stats @@ -1025,9 +1065,7 @@ def total_chi2_data_from_experiments(experiments_chi2_data, pdf): [exp_chi2_data.replica_result.data for exp_chi2_data in experiments_chi2_data], axis=0 ) - ndata = np.sum( - [exp_chi2_data.ndata for exp_chi2_data in experiments_chi2_data], - ) + ndata = np.sum([exp_chi2_data.ndata for exp_chi2_data in experiments_chi2_data]) return Chi2Data(pdf.stats_class(data_sum), central_result, ndata) @@ -1064,8 +1102,7 @@ def perreplica_chi2_table(groups_data, groups_chi2, total_chi2_data): total_chis[1:, :] /= np.array(ls)[:, np.newaxis] columns = pd.MultiIndex.from_arrays( - (["Total", *[str(exp) for exp in groups_data]], [total_n, *ls]), - names=["name", "npoints"], + (["Total", *[str(exp) for exp in groups_data]], [total_n, *ls]), names=["name", "npoints"] ) return pd.DataFrame(total_chis.T, columns=columns) @@ -1172,6 +1209,7 @@ def dataspecs_dataset_chi2_difference_table( dataspecs_groups_chi2_data = collect("groups_chi2", ("dataspecs",)) dataspecs_groups_bootstrap_phi = collect("groups_bootstrap_phi", ("dataspecs",)) dataspecs_results = collect("results", ("dataspecs",)) +dataspecs_results_with_scale_variations = collect("results_with_scale_variations", ("dataspecs",)) dataspecs_total_chi2 = collect("total_chi2_per_point_data", ("dataspecs",)) dataspecs_speclabel = collect("speclabel", ("dataspecs",), element_default=None) diff --git a/validphys2/src/validphys/theorycovariance/construction.py b/validphys2/src/validphys/theorycovariance/construction.py index 04d9cbfe4a..0bc03f798e 100644 --- a/validphys2/src/validphys/theorycovariance/construction.py +++ b/validphys2/src/validphys/theorycovariance/construction.py @@ -15,6 +15,7 @@ from reportengine import collect from reportengine.table import table from validphys.calcutils import all_chi2_theory, calc_chi2, central_chi2_theory +from validphys.checks import check_using_theory_covmat from validphys.results import Chi2Data, procs_central_values, procs_central_values_no_table, results from validphys.theorycovariance.theorycovarianceutils import ( check_correct_theory_combination, @@ -32,21 +33,9 @@ def make_scale_var_covmat(predictions): - """Takes N theory predictions at different scales and applies N-pt scale - variations to produce a covariance matrix.""" - l = len(predictions) - central, *others = predictions - deltas = (other - central for other in others) - if l == 3: - norm = 0.5 - elif l == 5: - norm = 0.5 - elif l == 7: - norm = 1 / 3 - elif l == 9: - norm = 0.25 - s = norm * sum(np.outer(d, d) for d in deltas) - return s + raise Exception( + "If you are seeing this error, please open an issue. This function should not be used" + ) @check_correct_theory_combination @@ -73,6 +62,38 @@ def theory_covmat_singleprocess(theory_covmat_singleprocess_no_table, fivetheori ) +@check_using_theory_covmat +def theory_covmat_dataset( + results, + results_bytheoryids, + use_theorycovmat, # for the check + point_prescription, + fivetheories=None, + seventheories=None, +): + """ + Compute the theory covmat for a collection of theoryids for a single dataset. + + In general this will come from some point prescription and it could be guessed from the input + however, it takes as input all relevant variables for generality + """ + _, theory_results = zip(*results_bytheoryids) + _, central_th_result = results + l = len(results_bytheoryids) + + # Remove the central theory from the list if it was included + theory_results = [i for i in theory_results if i._theoryid != central_th_result._theoryid] + cv = central_th_result.central_value + + # Compute the theory contribution to the covmats + deltas = list((t.central_value - cv for t in theory_results)) + thcovmat = compute_covs_pt_prescrip( + point_prescription, l, "A", deltas, fivetheories=fivetheories, seventheories=seventheories + ) + + return thcovmat + + @check_correct_theory_combination def theory_covmat_datasets(each_dataset_results_bytheory, fivetheories): """Produces an array of theory covariance matrices. Each matrix corresponds @@ -359,6 +380,109 @@ def covmat_n3lo_ad(name1, name2, deltas1, deltas2): return 1 / norm * s +def compute_covs_pt_prescrip( + point_prescription, + l, + name1, + deltas1, + name2=None, + deltas2=None, + fivetheories=None, + seventheories=None, +): + """Utility to compute the covariance matrix by prescription given the + shifts with respect to the central value for a pair of processes. + + The processes are defined by the variables ``name1`` and ``name2`` with + ``deltas1`` and ``deltas2`` the associated shifts wrt the central prediction. + + This utility also allows for the computation of the theory covmat for a single process + or dataset if ``names2=deltas2=None``. + + Parameters + --------- + point_prescription: str + defines the point prescription to be utilized + l: int + Number of theory variations (counting the central theory) + name1: str + Process name of the first set of shifts + deltas1: list(np.ndarray) + list of shifts for each of the non-central theories + name2: str + Process name of the second set of shifts + deltas2: list(np.ndarray) + list of shifts for each of the non-central theories + fivetheories: str + 5-point prescription variation + seventheories: str + 7-point prescription variation + """ + if name2 is None and deltas2 is not None: + raise ValueError( + f"Error building theory covmat: predictions have been given with no associated process/dataset name" + ) + elif deltas2 is None and name2 is not None: + raise ValueError( + f"Error building theory covmat: a process/dataset name has been given {name2} with no predictions" + ) + + if name2 is None: + name2 = name1 + deltas2 = deltas1 + + if l == 3: + if point_prescription == "3f point": + s = covmat_3fpt(name1, name2, deltas1, deltas2) + elif point_prescription == "3r point": + s = covmat_3rpt(name1, name2, deltas1, deltas2) + else: + s = covmat_3pt(name1, name2, deltas1, deltas2) + elif l == 5: + # 5 point -------------------------------------------------------------- + if fivetheories == "nobar": + s = covmat_5pt(name1, name2, deltas1, deltas2) + # 5bar-point ----------------------------------------------------------- + else: + s = covmat_5barpt(name1, name2, deltas1, deltas2) + # --------------------------------------------------------------------- + elif l == 7: + # Outdated 7pts implementation: left for posterity --------------------- + if seventheories == "original": + s = covmat_7pt_orig(name1, name2, deltas1, deltas2) + # 7pt (Gavin) ---------------------------------------------------------- + else: + s = covmat_7pt(name1, name2, deltas1, deltas2) + elif l == 9: + s = covmat_9pt(name1, name2, deltas1, deltas2) + # n3lo ad variation prescriprion + elif l == 62: + s = covmat_n3lo_singlet(name1, name2, deltas1, deltas2) + # n3lo ihou prescriprion + elif l == 64: + s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-2], deltas2[:-2]) + s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:]) + s = s_ad + s_cf + # n3lo 3 pt MHOU see also + # see https://github.com/NNPDF/papers/blob/e2ac1832cf4a36dab83a696564eaa75a4e55f5d2/minutes/minutes-2023-08-18.txt#L148-L157 + elif l == 66: + s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-4], deltas2[:-4]) + s_mhou = covmat_3pt(name1, name2, deltas1[-4:-2], deltas2[-4:-2]) + s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:]) + s = s_ad + s_cf + s_mhou + # n3lo full covmat prescriprion + elif l == 70: + # spit deltas and compose thcovmat + # splittin functions variatons + s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-8], deltas2[:-8]) + # scale variations + s_mhou = covmat_7pt(name1, name2, deltas1[-8:-2], deltas2[-8:-2]) + # massive coefficient function variations + s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:]) + s = s_ad + s_cf + s_mhou + return s + + @check_correct_theory_combination def covs_pt_prescrip( combine_by_type, @@ -385,55 +509,9 @@ def covs_pt_prescrip( deltas1 = list((other - central1 for other in others1)) central2, *others2 = process_info.theory[name2] deltas2 = list((other - central2 for other in others2)) - if l == 3: - if point_prescription == "3f point": - s = covmat_3fpt(name1, name2, deltas1, deltas2) - elif point_prescription == "3r point": - s = covmat_3rpt(name1, name2, deltas1, deltas2) - else: - s = covmat_3pt(name1, name2, deltas1, deltas2) - elif l == 5: - # 5 point -------------------------------------------------------------- - if fivetheories == "nobar": - s = covmat_5pt(name1, name2, deltas1, deltas2) - # 5bar-point ----------------------------------------------------------- - else: - s = covmat_5barpt(name1, name2, deltas1, deltas2) - # --------------------------------------------------------------------- - elif l == 7: - # Outdated 7pts implementation: left for posterity --------------------- - if seventheories == "original": - s = covmat_7pt_orig(name1, name2, deltas1, deltas2) - # 7pt (Gavin) ---------------------------------------------------------- - else: - s = covmat_7pt(name1, name2, deltas1, deltas2) - elif l == 9: - s = covmat_9pt(name1, name2, deltas1, deltas2) - # n3lo ad variation prescriprion - elif l == 62: - s = covmat_n3lo_singlet(name1, name2, deltas1, deltas2) - # n3lo ihou prescriprion - elif l == 64: - s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-2], deltas2[:-2]) - s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:]) - s = s_ad + s_cf - # n3lo 3 pt MHOU see also - # see https://github.com/NNPDF/papers/blob/e2ac1832cf4a36dab83a696564eaa75a4e55f5d2/minutes/minutes-2023-08-18.txt#L148-L157 - elif l == 66: - s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-4], deltas2[:-4]) - s_mhou = covmat_3pt(name1, name2, deltas1[-4:-2], deltas2[-4:-2]) - s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:]) - s = s_ad + s_cf + s_mhou - # n3lo full covmat prescriprion - elif l == 70: - # spit deltas and compose thcovmat - # splittin functions variatons - s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-8], deltas2[:-8]) - # scale variations - s_mhou = covmat_7pt(name1, name2, deltas1[-8:-2], deltas2[-8:-2]) - # massive coefficient function variations - s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:]) - s = s_ad + s_cf + s_mhou + s = compute_covs_pt_prescrip( + point_prescription, l, name1, deltas1, name2, deltas2, fivetheories, seventheories + ) start_locs = (start_proc[name1], start_proc[name2]) covmats[start_locs] = s return covmats