Merge pull request #638 from NNPDF/multiclosurev2

adding analysis tools for multiple closures/proxy fits

Author: Zaharid

validphys2/src/validphys/closuretest/__init__.py

@@ -5,3 +5,7 @@
 """
 from validphys.closuretest.closure_plots import *
 from validphys.closuretest.closure_results import *
+from validphys.closuretest.multiclosure import *
+from validphys.closuretest.multiclosure_output import *
+from validphys.closuretest.multiclosure_pdf import *
+from validphys.closuretest.multiclosure_pdf_output import *

validphys2/src/validphys/closuretest/closure_checks.py

@@ -1,8 +1,11 @@
 """
 closuretest/checks.py
 
-Module containing checks specfic to the closure tests
+Module containing checks specific to the closure tests.
+
 """
+from collections import defaultdict
+
 from reportengine.checks import make_argcheck, CheckError
 
 
@@ -51,3 +54,85 @@ def check_fits_areclosures(fits):
     for fit in fits:
         if not fit.as_input()["closuretest"]["fakedata"]:
             raise CheckError(f"Specified fit: {fit}, is not a closure test")
+
+
+@make_argcheck
+def check_t0pdfset_matches_law(t0pdfset, fit):
+    t0_from_fit = fit.as_input()["closuretest"]["fakepdf"]
+    if not str(t0pdfset) == t0_from_fit:
+        raise CheckError(
+            f"Underlying pdf: {t0_from_fit}, does not match t0pdfset: {t0pdfset}"
+        )
+
+
+@make_argcheck
+def check_at_least_10_fits(fits):
+    if len(fits) < 10:
+        raise CheckError(
+            "Multiclosure actions testing finite sampling effects require at least 10 fits"
+        )
+
+
+@make_argcheck
+def check_multifit_replicas(fits_pdf, _internal_max_reps, _internal_min_reps):
+    """Checks that all the fit pdfs have the same number of replicas N_rep. Then
+    check that N_rep is greater than the smallest number of replicas used in
+    actions which subsample the replicas of each fit.
+
+    This check also has the secondary
+    effect of filling in the namespace key _internal_max_reps
+    which can be used to override the number of replicas used at the level of
+    the runcard, but by default get filled in as the number of replicas in each
+    fit.
+
+    """
+    # we take off 1 here because we don't want to include replica 0
+    n_reps = {len(pdf) - 1 for pdf in fits_pdf}
+    if len(n_reps) != 1:
+        raise CheckError(
+            "all fits for multiclosure actions should have same number of replicas"
+        )
+    n_reps = n_reps.pop()
+    if _internal_max_reps is None:
+        _internal_max_reps = n_reps
+    elif _internal_max_reps > n_reps:
+        raise CheckError(
+            f"Specified _internal_max_reps to be {_internal_max_reps} "
+            f"however each fit only has {n_reps} replicas"
+        )
+
+    if _internal_max_reps < _internal_min_reps:
+        raise CheckError(
+            f"maximum replicas per fit, {_internal_max_reps}, is less than minimum replicas "
+            f", {_internal_min_reps}. If you have set _internal_max_reps and"
+            "_internal_min_reps then ensure that they take sensible values."
+        )
+    return {
+        "_internal_max_reps": _internal_max_reps,
+        "_internal_min_reps": _internal_min_reps,
+    }
+
+
+@make_argcheck
+def check_fits_different_filterseed(fits):
+    """Input fits should have different filter seeds if they are being
+    used for multiple closure test studies, because in high-level hand-waving
+    terms the different level 1 shifts represents different
+    'runs of the universe'!
+
+    """
+    seed_fits_dict = defaultdict(list)
+
+    for fit in fits:
+        pdf_name = fit.as_input()["pdf"]["id"]
+        seed = fit.as_input()["closuretest"]["filterseed"]
+        seed_fits_dict[seed].append(pdf_name)
+
+    bad_fits = [fits for _, fits in seed_fits_dict.items() if len(fits) > 1]
+
+    if bad_fits:
+        raise CheckError(
+        "Multiclosure actions require that fits have different level 1 "
+        "noise and therefore different filter seeds. The following groups "
+        f"of fits have the same seed: {bad_fits}."
+        )

validphys2/src/validphys/closuretest/closure_plots.py

@@ -24,6 +24,7 @@ def plot_biases(biases_table):
     ax.legend()
     return fig
 
+
 @figure
 def plot_delta_chi2(delta_chi2_bootstrap, fits):
     """Plots distributions of delta chi2 for each fit in `fits`.
@@ -41,6 +42,7 @@ def plot_delta_chi2(delta_chi2_bootstrap, fits):
     ax.set_title(r"Total $\Delta_{\chi^{2}}$ for each fit")
     return fig
 
+
 def errorbar_figure_from_table(df):
     """Given a table with even columns as central values as odd columns as errors
     plot an errorbar plot"""
@@ -49,9 +51,11 @@ def errorbar_figure_from_table(df):
         df.values[:, 1::2].T,
         df.index.values,
         df.columns.unique(0),
-        xlim=0)
+        xlim=0,
+    )
     return fig, ax
 
+
 @figure
 def plot_fits_bootstrap_variance(fits_bootstrap_variance_table):
     """Plot variance as error bars, with mean and central value calculated
@@ -61,10 +65,9 @@ def plot_fits_bootstrap_variance(fits_bootstrap_variance_table):
     ax.set_title("Variance by experiment for closure fits")
     return fig
 
+
 @figure
-def plot_fits_bootstrap_bias(
-    fits_bootstrap_bias_table
-):
+def plot_fits_bootstrap_bias(fits_bootstrap_bias_table):
     """Plot the bias for each experiment for all `fits` as a point with an error bar,
     where the error bar is given by bootstrapping the bias across replicas
 

validphys2/src/validphys/closuretest/closure_results.py

@@ -8,6 +8,9 @@
 import numpy as np
 import pandas as pd
 
+from reportengine import collect
+from reportengine.table import table
+
 from validphys.calcutils import calc_chi2, bootstrap_values
 from validphys.checks import check_pdf_is_montecarlo
 from validphys.closuretest.closure_checks import (
@@ -17,19 +20,16 @@
     check_fits_same_filterseed,
     check_fits_underlying_law_match,
 )
-from reportengine import collect
-from reportengine.table import table
 
 
 BiasData = namedtuple("BiasData", ("bias", "ndata"))
 
 underlying_results = collect("results", ("fitunderlyinglaw",))
 
+
 @check_fit_isclosure
 @check_use_fitcommondata
-def bias_dataset(
-    results, underlying_results, fit, use_fitcommondata, sqrt_covariance_matrix
-):
+def bias_dataset(results, underlying_results, fit, use_fitcommondata):
     """Calculate the bias for a given dataset and fit. The bias is defined as
     chi2 between the prediction from the underlying PDF (which was used to
     generate the closure pseudodata), also known as level zero closure data, and
@@ -39,11 +39,11 @@ def bias_dataset(
     is used to generate the multiplicative contributions to the covariance
     matrix
     """
-    _, th_ct = results
+    dt_ct, th_ct = results
     # does collect need to collect a list even with one element?
     (_, th_ul), = underlying_results
     central_diff = th_ct.central_value - th_ul.central_value
-    bias_out = calc_chi2(sqrt_covariance_matrix, central_diff)  # unnormalised
+    bias_out = calc_chi2(dt_ct.sqrtcovmat, central_diff)  # unnormalised
     return BiasData(bias_out, len(th_ct))
 
 
@@ -53,20 +53,12 @@ def bias_dataset(
 @check_fit_isclosure
 @check_use_fitcommondata
 def bias_experiment(
-    experiment_results,
-    underlying_experiment_results,
-    fit,
-    use_fitcommondata,
-    experiment_sqrt_covariance_matrix,
+    experiment_results, underlying_experiment_results, fit, use_fitcommondata
 ):
     """Like `bias_dataset` but for a whole experiment.
     """
     return bias_dataset(
-        experiment_results,
-        underlying_experiment_results,
-        fit,
-        use_fitcommondata,
-        experiment_sqrt_covariance_matrix,
+        experiment_results, underlying_experiment_results, fit, use_fitcommondata
     )
 
 
@@ -177,7 +169,7 @@ def fits_bootstrap_bias_table(
 
 @check_fit_isclosure
 @check_use_fitcommondata
-def variance_dataset(results, fit, use_fitcommondata, sqrt_covariance_matrix):
+def variance_dataset(results, fit, use_fitcommondata):
     """calculate the variance for a given dataset, which is the spread of
     replicas measured in the space of the covariance matrix. Given by:
 
@@ -188,25 +180,18 @@ def variance_dataset(results, fit, use_fitcommondata, sqrt_covariance_matrix):
     be made fully independent of the closure data. This is useful when checking
     the variance of data which was not included in the fit.
 
-    # TODO: here we require that use_fitcommondata is true, for the generic use
-    # case. we require another action which uses explicitly a t0pdf of the
-    # underlying law.
     """
-    _, th = results
+    dt, th = results
     diff = th.central_value[:, np.newaxis] - th._rawdata
-    var_unnorm = calc_chi2(sqrt_covariance_matrix, diff).mean()
+    var_unnorm = calc_chi2(dt.sqrtcovmat, diff).mean()
     return VarianceData(var_unnorm, len(th))
 
 
 @check_fit_isclosure
 @check_use_fitcommondata
-def variance_experiment(
-    experiment_results, fit, use_fitcommondata, experiment_sqrt_covariance_matrix
-):
+def variance_experiment(experiment_results, fit, use_fitcommondata):
     """Like variance_dataset but for a whole experiment"""
-    return variance_dataset(
-        experiment_results, fit, use_fitcommondata, experiment_sqrt_covariance_matrix
-    )
+    return variance_dataset(experiment_results, fit, use_fitcommondata)
 
 
 def bootstrap_variance_experiment(experiment_results, bootstrap_samples=500):