diff --git a/doc/sphinx/source/vp/theorycov/examples.rst b/doc/sphinx/source/vp/theorycov/examples.rst index 3e3aa7a341..e7bde9f987 100644 --- a/doc/sphinx/source/vp/theorycov/examples.rst +++ b/doc/sphinx/source/vp/theorycov/examples.rst @@ -154,39 +154,6 @@ a comprehensive set of plots and tables describing the covariance matrices. {@plot_diag_cov_comparison@} {@endwith@} - Experimental $\chi^2$ - --------------------- - {@with default_theory@} - {@total_experiments_chi2@} - - Total (exp. + th.) $\chi^2$ - --------------------------- - {@chi2_impact_custom@} - - Experimental $\chi^2$ by dataset - -------------------------------- - {@experiments_chi2_table@} - - Total (exp. + th.) $\chi^2$ by dataset - -------------------------------------- - {@experiments_chi2_table_theory@} - - $\chi^2$ including only diagonal theory elements - ------------------------------------------------ - {@chi2_diag_only@} - - Impact of theory covariance matrix on $\chi^2$s - ----------------------------------------------- - {@plot_datasets_chi2_theory@} - {@endwith@} - - Scale variations as a function of the kinematics - ------------------------------------------------ - {@with matched_datasets_from_dataspecs@} - [Plots for {@dataset_name@}]({@dataset_report report@}) - {@endwith@} - - Validation report ----------------- diff --git a/validphys2/examples/theory_covariance/template_matrix_plots.md b/validphys2/examples/theory_covariance/template_matrix_plots.md index 8e200fb46f..7b29c92016 100644 --- a/validphys2/examples/theory_covariance/template_matrix_plots.md +++ b/validphys2/examples/theory_covariance/template_matrix_plots.md @@ -18,35 +18,3 @@ Diagonal elements of covariance matrices {@with default_theory@} {@plot_diag_cov_comparison@} {@endwith@} - -Experimental $\chi^2$ ---------------------- -{@with default_theory@} - {@total_chi2_per_point_data@} - -Total (exp. + th.) $\chi^2$ ---------------------------- - {@chi2_impact_custom@} - -Experimental $\chi^2$ by dataset --------------------------------- - {@procs_chi2_table@} - -Total (exp. + th.) $\chi^2$ by dataset --------------------------------------- - {@procs_chi2_table_theory@} - -$\chi^2$ including only diagonal theory elements ------------------------------------------------- - {@chi2_diag_only@} - -Impact of theory covariance matrix on $\chi^2$s ------------------------------------------------ - {@plot_datasets_chi2_theory@} -{@endwith@} - -Scale variations as a function of the kinematics ------------------------------------------------- -{@with matched_datasets_from_dataspecs@} - [Plots for {@dataset_name@}]({@dataset_report report@}) -{@endwith@} diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py index f131ad3f6a..d9dc6b743b 100644 --- a/validphys2/src/validphys/pdfgrids.py +++ b/validphys2/src/validphys/pdfgrids.py @@ -346,7 +346,7 @@ def distance_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) = This method returns distance grids where the relative distance between both PDF set is computed. At least one grid will be identical to zero. """ - + gr2_stats = xplotting_grids[normalize_to].grid_values cv2 = gr2_stats.central_value() sg2 = gr2_stats.std_error() @@ -374,6 +374,46 @@ def distance_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) = return newgrids +@check_pdf_normalize_to +def pull_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) = None): + """Return an object containing the value of the distance PDF at the specified values + of x and flavour. + + The parameter ``normalize_to`` identifies the reference PDF set with respect to the + distance is computed. + + This method returns distance grids where the relative distance between both PDF + set is computed. At least one grid will be identical to zero. + """ + + gr2_stats = xplotting_grids[normalize_to].grid_values + cv2 = gr2_stats.central_value() + sg2 = gr2_stats.std_error() + N2 = pdfs[normalize_to].get_members() + + newgrids = [] + for grid, pdf in zip(xplotting_grids, pdfs): + if pdf == pdfs[normalize_to]: + # Zero the PDF we are normalizing against + pdf_zero = pdf.stats_class(np.zeros_like(gr2_stats.data[0:1])) + newgrid = grid.copy_grid(grid_values=pdf_zero) + newgrids.append(newgrid) + continue + + g_stats = grid.grid_values + cv1 = g_stats.central_value() + sg1 = g_stats.std_error() + N1 = pdf.get_members() + + # Wrap the distance into a Stats (1, flavours, points) + distance = Stats([np.sqrt((cv1 - cv2) ** 2 / (sg1**2 + sg2**2))]) + + newgrid = grid.copy_grid(grid_values=distance) + newgrids.append(newgrid) + + return newgrids + +pull_grids_list = collect(pull_grids, ('pdfs_list',)) @check_pdf_normalize_to def variance_distance_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) = None): diff --git a/validphys2/src/validphys/pdfplots.py b/validphys2/src/validphys/pdfplots.py index 80300060fa..12c134c724 100644 --- a/validphys2/src/validphys/pdfplots.py +++ b/validphys2/src/validphys/pdfplots.py @@ -273,7 +273,106 @@ def plot_pdf_uncertainties( PDF's central value is plotted. Otherwise it is the absolute values.""" yield from UncertaintyPDFPlotter(pdfs, xplotting_grids, xscale, normalize_to, ymin, ymax) +class PullPDFPlotter(metaclass=abc.ABCMeta): + """Auxiliary class which groups multiple pulls in one plot.""" + + def __init__(self, pdfs_list, pull_grids_list, xscale, normalize_to, ymin, ymax): + self.pdfs_list = pdfs_list + self.pull_grids_list = pull_grids_list + self._xscale = xscale + self.normalize_to = normalize_to + self.ymin = ymin + self.ymax = ymax + self.firstgrid = pull_grids_list[0][0] + + def legend(self, flstate): + return flstate.ax.legend() + + def get_ylabel(self): + return "Pull" + + @property + def Q(self): + return self.firstgrid.Q + + @property + def xscale(self): + if self._xscale is None: + return scale_from_grid(self.firstgrid) + return self._xscale + + def get_title(self, flstate): + return '$%s$' % flstate.parton_name + f', Q={self.Q : .1f} GeV' + + def draw(self, pdfs, grid, flstate): + ax = flstate.ax + flindex = flstate.flindex + color = ax._get_lines.get_next_color() + + # The grid for the distance is (1, flavours, points) + # take only the flavour we are interested in + gv = grid.select_flavour(flindex).grid_values.data.squeeze() + + ax.plot(grid.xgrid, gv, color=color, label = f'{pdfs[0].label}-{pdfs[1].label} pull') + + return gv + + def plot_call(self): + basis = self.firstgrid.basis + for flindex, fl in enumerate(self.firstgrid.flavours): + fig, ax = plotutils.subplots() + parton_name = basis.elementlabel(fl) + flstate = FlavourState(flindex=flindex, fl=fl, fig=fig, ax=ax, parton_name=parton_name) + ax.set_title(self.get_title(flstate)) + + all_vals = [] + for pdf, grids in zip(self.pdfs_list, self.pull_grids_list): + limits = self.draw([pdf['pdfs'][0],pdf['pdfs'][1]], grids[1], flstate) + if limits is not None: + all_vals.append(np.atleast_2d(limits)) + # Note these two lines do not conmute! + ax.set_xscale(self.xscale) + plotutils.frame_center(ax, self.firstgrid.xgrid, np.concatenate(all_vals)) + if self.ymin is not None: + ax.set_ylim(ymin=self.ymin) + if self.ymax is not None: + ax.set_ylim(ymax=self.ymax) + + ax.set_xlabel('$x$') + ax.set_xlim(self.firstgrid.xgrid[0]) + + ax.set_ylabel(self.get_ylabel()) + + ax.set_axisbelow(True) + + self.legend(flstate) + yield fig, parton_name + + def __call__(self): + for fig, parton_name in self.plot_call(): + ax = fig.get_axes()[0] + ymin, _ = ax.get_ylim() + ax.set_ylim(max(0, ymin), None) + yield fig, parton_name + + +@figuregen +@check_scale('xscale', allow_none=True) +def plot_pdf_pulls( + pdfs_list, + pull_grids_list, + xscale: (str, type(None)) = None, + normalize_to: (int, str, type(None)) = None, + ymin=None, + ymax=None, +): + """Plot the PDF standard deviations as a function of x. + If normalize_to is set, the ratio to that + PDF's central value is plotted. Otherwise it is the absolute values.""" + yield from PullPDFPlotter(pdfs_list, pull_grids_list, xscale, normalize_to, ymin, ymax)() + + class AllFlavoursPlotter(PDFPlotter): """Auxiliary class which groups multiple PDF flavours in one plot.""" diff --git a/validphys2/src/validphys/theorycovariance/construction.py b/validphys2/src/validphys/theorycovariance/construction.py index 147568daea..4f4ff2db59 100644 --- a/validphys2/src/validphys/theorycovariance/construction.py +++ b/validphys2/src/validphys/theorycovariance/construction.py @@ -31,43 +31,9 @@ log = logging.getLogger(__name__) -theoryids_procs_central_values = collect(procs_central_values, ("theoryids",)) - -theoryids_procs_central_values_no_table = collect(procs_central_values_no_table, ("theoryids",)) collected_theoryids = collect("theoryids", ["theoryconfig"]) - -def make_scale_var_covmat(predictions): - raise Exception( - "If you are seeing this error, please open an issue. This function should not be used" - ) - - -@check_correct_theory_combination -def theory_covmat_singleprocess_no_table( - theoryids_procs_central_values_no_table, procs_index, theoryids, fivetheories -): - """Calculates the theory covariance matrix for scale variations. - The matrix is a dataframe indexed by procs_index.""" - s = make_scale_var_covmat(theoryids_procs_central_values_no_table) - df = pd.DataFrame(s, index=procs_index, columns=procs_index) - return df - - -@table -@check_correct_theory_combination -def theory_covmat_singleprocess(theory_covmat_singleprocess_no_table, fivetheories): - """Duplicate of theory_covmat_singleprocess_no_table but with a table decorator.""" - return theory_covmat_singleprocess_no_table - - -results_central_bytheoryids = collect(results_central, ("theoryids",)) -each_dataset_results_central_bytheory = collect( - "results_central_bytheoryids", ("group_dataset_inputs_by_process", "data") -) - - @check_using_theory_covmat def theory_covmat_dataset( results, @@ -99,89 +65,9 @@ def theory_covmat_dataset( return thcovmat - -@check_correct_theory_combination -def theory_covmat_datasets(each_dataset_results_central_bytheory, fivetheories): - """Produces an array of theory covariance matrices. Each matrix corresponds - to a different dataset, which must be specified in the runcard.""" - dataset_covmats = [] - for dataset in each_dataset_results_central_bytheory: - theory_centrals = [x[1].central_value for x in dataset] - s = make_scale_var_covmat(theory_centrals) - dataset_covmats.append(s) - return dataset_covmats - - -@check_correct_theory_combination -def total_covmat_datasets(each_dataset_results_central_bytheory, fivetheories): - """Produces an array of total covariance matrices; the sum of experimental - and scale-varied theory covariance matrices. Each matrix corresponds - to a different dataset, which must be specified in the runcard. - These are needed for calculation of chi2 per dataset.""" - dataset_covmats = [] - for dataset in each_dataset_results_central_bytheory: - theory_centrals = [x[1].central_value for x in dataset] - s = make_scale_var_covmat(theory_centrals) - sigma = dataset[0][0].covmat - cov = s + sigma - dataset_covmats.append(cov) - return dataset_covmats - - -@check_correct_theory_combination -def total_covmat_diagtheory_datasets(each_dataset_results_central_bytheory, fivetheories): - """Same as total_covmat_theory_datasets but for diagonal theory only""" - dataset_covmats = [] - for dataset in each_dataset_results_central_bytheory: - theory_centrals = [x[1].central_value for x in dataset] - s = make_scale_var_covmat(theory_centrals) - # Initialise array of zeros and set precision to same as FK tables - s_diag = np.zeros((len(s), len(s)), dtype=np.float32) - np.fill_diagonal(s_diag, np.diag(s)) - sigma = dataset[0][0].covmat - cov = s_diag + sigma - dataset_covmats.append(cov) - return dataset_covmats - - -@table -def theory_block_diag_covmat(theory_covmat_datasets, procs_index): - """Takes the theory covariance matrices for individual datasets and - returns a data frame with a block diagonal theory covariance matrix - by dataset""" - s = la.block_diag(*theory_covmat_datasets) - df = pd.DataFrame(s, index=procs_index, columns=procs_index) - return df - - -@table -def theory_diagonal_covmat(theory_covmat_singleprocess_no_table, procs_index): - """Returns theory covmat with only diagonal values""" - s = theory_covmat_singleprocess_no_table.values - # Initialise array of zeros and set precision to same as FK tables - s_diag = np.zeros((len(s), len(s)), dtype=np.float32) - np.fill_diagonal(s_diag, np.diag(s)) - df = pd.DataFrame(s_diag, index=procs_index, columns=procs_index) - return df - - procs_results_theory = collect("procs_results", ("theoryids",)) -@check_correct_theory_combination -def total_covmat_procs(procs_results_theory, fivetheories): - """Same as total_covmat_datasets but per experiment rather than - per dataset. Needed for calculation of chi2 per experiment.""" - proc_result_covmats = [] - for proc_result in zip(*procs_results_theory): - theory_centrals = [x[1].central_value for x in proc_result] - s = make_scale_var_covmat(theory_centrals) - sigma = proc_result[0][0].covmat - cov = s + sigma - proc_result_covmats.append(cov) - return proc_result_covmats - - def dataset_names(data_input): """Returns a list of the names of the datasets, in the same order as they are inputted in the runcard""" @@ -189,7 +75,10 @@ def dataset_names(data_input): ProcessInfo = namedtuple("ProcessInfo", ("theory", "namelist", "sizes")) - +results_central_bytheoryids = collect(results_central, ("theoryids",)) +each_dataset_results_central_bytheory = collect( + "results_central_bytheoryids", ("group_dataset_inputs_by_process", "data") +) def combine_by_type(each_dataset_results_central_bytheory, dataset_names): """Groups the datasets according to processes and returns three objects: @@ -696,65 +585,14 @@ def procs_index_matched(groups_index, procs_index): return pd.MultiIndex.from_tuples(tups, names=("process", "dataset", "id")) - -@check_correct_theory_combination -def total_covmat_diagtheory_procs(procs_results_theory, fivetheories): - """Same as total_covmat_datasets but per proc rather than - per dataset. Needed for calculation of chi2 per proc.""" - exp_result_covmats = [] - for exp_result in zip(*procs_results_theory): - theory_centrals = [x[1].central_value for x in exp_result] - s = make_scale_var_covmat(theory_centrals) - # Initialise array of zeros and set precision to same as FK tables - s_diag = np.zeros((len(s), len(s)), dtype=np.float32) - np.fill_diagonal(s_diag, np.diag(s)) - sigma = exp_result[0][0].covmat - cov = s_diag + sigma - exp_result_covmats.append(cov) - return exp_result_covmats - - -@table -def theory_corrmat_singleprocess(theory_covmat_singleprocess): - """Calculates the theory correlation matrix for scale variations.""" - df = theory_covmat_singleprocess - covmat = df.values - diag_minus_half = (np.diagonal(covmat)) ** (-0.5) - mat = diag_minus_half[:, np.newaxis] * df * diag_minus_half - return mat - - -@table -def theory_blockcorrmat(theory_block_diag_covmat): - """Calculates the theory correlation matrix for scale variations - with block diagonal entries by dataset only""" - mat = theory_corrmat_singleprocess(theory_block_diag_covmat) - return mat - - @table def theory_corrmat_custom(theory_covmat_custom): """Calculates the theory correlation matrix for scale variations with variations by process type""" - mat = theory_corrmat_singleprocess(theory_covmat_custom) - return mat - - -@table -def theory_normcovmat_singleprocess(theory_covmat_singleprocess, procs_data_values): - """Calculates the theory covariance matrix for scale variations normalised - to data.""" - df = theory_covmat_singleprocess - mat = df / np.outer(procs_data_values, procs_data_values) - return mat - - -@table -def theory_normblockcovmat(theory_block_diag_covmat, procs_data_values): - """Calculates the theory covariance matrix for scale variations - normalised to data, block diagonal by dataset.""" - df = theory_block_diag_covmat - mat = df / np.outer(procs_data_values, procs_data_values) + df = theory_covmat_custom + covmat = df.values + diag_minus_half = (np.diagonal(covmat)) ** (-0.5) + mat = diag_minus_half[:, np.newaxis] * df * diag_minus_half return mat @@ -766,25 +604,6 @@ def theory_normcovmat_custom(theory_covmat_custom, procs_data_values): mat = df / np.outer(procs_data_values, procs_data_values) return mat - -@table -def experimentplustheory_covmat_singleprocess( - procs_covmat_no_table, theory_covmat_singleprocess_no_table -): - """Calculates the experiment + theory covariance matrix for - scale variations.""" - df = procs_covmat_no_table + theory_covmat_singleprocess_no_table - return df - - -@table -def experimentplusblocktheory_covmat(procs_covmat, theory_block_diag_covmat): - """Calculates the experiment + theory covariance - matrix for scale variations.""" - df = procs_covmat + theory_block_diag_covmat - return df - - @table def experimentplustheory_covmat_custom(procs_covmat, theory_covmat_custom): """Calculates the experiment + theory covariance matrix for @@ -793,28 +612,6 @@ def experimentplustheory_covmat_custom(procs_covmat, theory_covmat_custom): return df -@table -def experimentplustheory_normcovmat_singleprocess( - procs_covmat, theory_covmat_singleprocess, procs_data -): - """Calculates the experiment + theory covariance matrix for scale - variations normalised to data.""" - df = procs_covmat + theory_covmat_singleprocess - procs_data_array = np.array(procs_data) - mat = df / np.outer(procs_data_array, procs_data_array) - return mat - - -@table -def experimentplusblocktheory_normcovmat( - procs_covmat, theory_block_diag_covmat, procs_data_values, experimentplustheory_normcovmat -): - """Calculates the experiment + theory covariance matrix for scale - variations normalised to data, block diagonal by data set.""" - mat = experimentplustheory_normcovmat(procs_covmat, theory_block_diag_covmat, procs_data_values) - return mat - - @table def experimentplustheory_normcovmat_custom( procs_covmat, theory_covmat_custom, procs_data_values, experimentplustheory_normcovmat @@ -825,121 +622,14 @@ def experimentplustheory_normcovmat_custom( return mat - -@table -def experimentplustheory_corrmat_singleprocess(procs_covmat, theory_covmat_singleprocess): - """Calculates the correlation matrix for the experimental - plus theory covariance matrices.""" - total_df = procs_covmat + theory_covmat_singleprocess - total_cov = (procs_covmat + theory_covmat_singleprocess).values - diag_minus_half = (np.diagonal(total_cov)) ** (-0.5) - corrmat = diag_minus_half[:, np.newaxis] * total_df * diag_minus_half - return corrmat - - -@table -def experimentplusblocktheory_corrmat(procs_covmat, theory_block_diag_covmat): - """Calculates the correlation matrix for the experimental - plus theory covariance matrices, block diagonal by dataset.""" - corrmat = experimentplustheory_corrmat_singleprocess(procs_covmat, theory_block_diag_covmat) - return corrmat - - @table def experimentplustheory_corrmat_custom(procs_covmat, theory_covmat_custom): """Calculates the correlation matrix for the experimental plus theory covariance matrices, correlations by prescription.""" - corrmat = experimentplustheory_corrmat_singleprocess(procs_covmat, theory_covmat_custom) + total_df = procs_covmat + theory_covmat_custom + total_cov = (procs_covmat + theory_covmat_custom).values + diag_minus_half = (np.diagonal(total_cov)) ** (-0.5) + corrmat = diag_minus_half[:, np.newaxis] * total_df * diag_minus_half return corrmat - -def chi2_impact(theory_covmat_singleprocess, procs_covmat, procs_results): - """Returns total chi2 including theory cov mat""" - dataresults, theoryresults = zip(*procs_results) - dat_central_list = [x.central_value for x in dataresults] - th_central_list = [x.central_value for x in theoryresults] - dat_central = np.concatenate(dat_central_list) - th_central = np.concatenate([x for x in th_central_list]) - central_diff = dat_central - th_central - cov = theory_covmat_singleprocess.values + procs_covmat.values - return calc_chi2(la.cholesky(cov, lower=True), central_diff) / len(central_diff) - - -def data_theory_diff(procs_results): - """Returns (D-T) for central theory, for use in chi2 calculations""" - dataresults, theoryresults = zip(*procs_results) - dat_central_list = [x.central_value for x in dataresults] - th_central_list = [x.central_value for x in theoryresults] - dat_central = np.concatenate(dat_central_list) - th_central = np.concatenate(th_central_list) - central_diff = dat_central - th_central - return central_diff - - -def chi2_block_impact(theory_block_diag_covmat, procs_covmat, procs_results): - """Returns total chi2 including theory cov mat""" - chi2 = chi2_impact(theory_block_diag_covmat, procs_covmat, procs_results) - return chi2 - - -def chi2_impact_custom(theory_covmat_custom, procs_covmat, procs_results): - """Returns total chi2 including theory cov mat""" - chi2 = chi2_impact(theory_covmat_custom, procs_covmat, procs_results) - return chi2 - - -def theory_diagcovmat(theory_covmat_singleprocess): - """Returns theory covmat with only diagonal values""" - s = theory_covmat_singleprocess.values - # Initialise array of zeros and set precision to same as FK tables - s_diag = np.zeros((len(s), len(s)), dtype=np.float32) - np.fill_diagonal(s_diag, np.diag(s)) - return s_diag - - -def chi2_diag_only(theory_diagcovmat, procs_covmat, data_theory_diff): - """Returns total chi2 including only diags of theory cov mat""" - cov = theory_diagcovmat + procs_covmat.values - elements = np.dot(data_theory_diff.T, np.dot(la.inv(cov), data_theory_diff)) - chi2 = (1 / len(data_theory_diff)) * np.sum(elements) - return chi2 - - -each_dataset_results = collect(results, ("group_dataset_inputs_by_process", "data")) - - -def abs_chi2_data_theory_dataset(each_dataset_results, total_covmat_datasets): - """Returns an array of tuples (member_chi², central_chi², numpoints) - corresponding to each data set, where theory errors are included""" - chi2data_array = [] - for datresults, covmat in zip(each_dataset_results, total_covmat_datasets): - data_result, th_result = datresults - chi2s = all_chi2_theory(datresults, covmat) - central_result = central_chi2_theory(datresults, covmat) - chi2data_array.append( - Chi2Data(th_result.stats_class(chi2s[:, np.newaxis]), central_result, len(data_result)) - ) - return chi2data_array - - -def abs_chi2_data_theory_proc(procs_results, total_covmat_procs): - """Like abs_chi2_data_theory_dataset but for procs not datasets""" - chi2data_array = [] - for expresults, covmat in zip(procs_results, total_covmat_procs): - data_result, th_result = expresults - chi2s = all_chi2_theory(expresults, covmat) - central_result = central_chi2_theory(expresults, covmat) - chi2data_array.append( - Chi2Data(th_result.stats_class(chi2s[:, np.newaxis]), central_result, len(data_result)) - ) - return chi2data_array - - -def abs_chi2_data_diagtheory_proc(procs_results, total_covmat_diagtheory_procs): - """For a diagonal theory covmat""" - return abs_chi2_data_theory_proc(procs_results, total_covmat_diagtheory_procs) - - -def abs_chi2_data_diagtheory_dataset(each_dataset_results, total_covmat_diagtheory_datasets): - """For a diagonal theory covmat""" - return abs_chi2_data_theory_dataset(each_dataset_results, total_covmat_diagtheory_datasets) +each_dataset_results = collect(results, ("group_dataset_inputs_by_process", "data")) \ No newline at end of file diff --git a/validphys2/src/validphys/theorycovariance/output.py b/validphys2/src/validphys/theorycovariance/output.py index 0487d81215..158f31f8c1 100644 --- a/validphys2/src/validphys/theorycovariance/output.py +++ b/validphys2/src/validphys/theorycovariance/output.py @@ -22,33 +22,6 @@ log = logging.getLogger(__name__) -abs_chi2_data_theory_dataset_by_process = collect( - "abs_chi2_data_theory_dataset", ("group_dataset_inputs_by_process",) -) - - -@table -def procs_chi2_table_theory( - procs_data, pdf, abs_chi2_data_theory_proc, abs_chi2_data_theory_dataset_by_process -): - """Same as groups_chi2_table but including theory covariance matrix. - Note: we use groups_chi2_table here but provide data grouped by process.""" - return groups_chi2_table( - procs_data, pdf, abs_chi2_data_theory_proc, abs_chi2_data_theory_dataset_by_process - ) - - -@table -def procs_chi2_table_diagtheory( - procs_data, pdf, abs_chi2_data_diagtheory_proc, abs_chi2_data_theory_dataset_by_process -): - """Same as groups_chi2_table but including diagonal theory covariance matrix. - Note: we use groups_chi2_table here but provide data grouped by process.""" - return groups_chi2_table( - procs_data, pdf, abs_chi2_data_diagtheory_proc, abs_chi2_data_theory_dataset_by_process - ) - - def matrix_plot_labels(df): """Returns the tick locations and labels, and the starting point values for each category, based on a dataframe @@ -222,16 +195,6 @@ def plot_expcorrmat_heatmap(procs_corrmat): fig = plot_corrmat_heatmap(procs_corrmat, "Experimental Correlation Matrix") return fig - -@figure -def plot_normthblockcovmat_heatmap(theory_normblockcovmat): - """Matrix plot for block diagonal theory covariance matrix""" - fig = plot_covmat_heatmap( - theory_normblockcovmat, "Block diagonal theory covariance matrix by dataset" - ) - return fig - - @figure def plot_normthcovmat_heatmap_custom(theory_normcovmat_custom, theoryids, fivetheories): """Matrix plot for block diagonal theory covariance matrix by process type""" @@ -243,15 +206,6 @@ def plot_normthcovmat_heatmap_custom(theory_normcovmat_custom, theoryids, fiveth return fig -@figure -def plot_thblockcorrmat_heatmap(theory_blockcorrmat): - """Matrix plot of the theory correlation matrix""" - fig = plot_corrmat_heatmap( - theory_blockcorrmat, "Theory correlation matrix block diagonal by dataset" - ) - return fig - - @figure def plot_thcorrmat_heatmap_custom(theory_corrmat_custom, theoryids, fivetheories): """Matrix plot of the theory correlation matrix, correlations by process type""" @@ -262,42 +216,6 @@ def plot_thcorrmat_heatmap_custom(theory_corrmat_custom, theoryids, fivetheories fig = plot_corrmat_heatmap(theory_corrmat_custom, f"Theory Correlation matrix ({l} pt)") return fig - -@figure -def plot_normexpplusblockthcovmat_heatmap(experimentplusblocktheory_normcovmat): - """Matrix plot of the exp + theory covariance matrix normalised to data""" - fig = plot_covmat_heatmap( - experimentplusblocktheory_normcovmat, - "Experiment + theory (block diagonal by dataset) covariance matrix", - ) - return fig - - -@figure -def plot_normexpplusthcovmat_heatmap_custom( - experimentplustheory_normcovmat_custom, theoryids, fivetheories -): - """Matrix plot of the exp + theory covariance matrix normalised to data""" - l = len(theoryids) - if l == 5: - if fivetheories == "bar": - l = r"$\bar{5}$" - fig = plot_covmat_heatmap( - experimentplustheory_normcovmat_custom, f"Experimental + Theory Covariance Matrix ({l} pt)" - ) - return fig - - -@figure -def plot_expplusblockthcorrmat_heatmap(experimentplusblocktheory_corrmat): - """Matrix plot of the exp + theory correlation matrix""" - fig = plot_corrmat_heatmap( - experimentplusblocktheory_corrmat, - "Experiment + theory (block diagonal by dataset) correlation matrix", - ) - return fig - - @figure def plot_expplusthcorrmat_heatmap_custom( experimentplustheory_corrmat_custom, theoryids, fivetheories @@ -312,31 +230,6 @@ def plot_expplusthcorrmat_heatmap_custom( ) return fig - -@figure -def plot_blockcovdiff_heatmap(theory_block_diag_covmat, procs_covmat): - """Matrix plot (thcov + expcov)/expcov""" - df = (theory_block_diag_covmat.as_matrix() + procs_covmat.values) / np.mean(procs_covmat.values) - fig = plot_covmat_heatmap( - df, "(Theory + experiment)/mean(experiment)" + "for block diagonal theory covmat by dataset" - ) - return fig - - -@figure -def plot_covdiff_heatmap_custom(theory_covmat_custom, procs_covmat, theoryids, fivetheories): - """Matrix plot (thcov + expcov)/expcov""" - l = len(theoryids) - if l == 5: - if fivetheories == "bar": - l = r"$\bar{5}$" - df = (theory_covmat_custom + procs_covmat) / np.mean(procs_covmat.values) - fig = plot_covmat_heatmap( - df, f"(Theory + experiment)/mean(experiment) covariance matrices for {l} points" - ) - return fig - - @figure def plot_diag_cov_comparison( theory_covmat_custom, procs_covmat, procs_data_values, theoryids, fivetheories @@ -384,64 +277,3 @@ def plot_diag_cov_comparison( ax.legend(fontsize=20) ax.margins(x=0) return fig - - -@figure -def plot_diag_cov_impact( - theory_covmat_custom, procs_covmat, procs_data_values, theoryids, fivetheories -): - """Plot ((expcov)^-1_ii)^-0.5 versus ((expcov + thcov)^-1_ii)^-0.5""" - l = len(theoryids) - if l == 5: - if fivetheories == "bar": - l = r"$\bar{5}$" - matrix_theory = theory_covmat_custom.values - matrix_experiment = procs_covmat.values - inv_exp = (np.diag(la.inv(matrix_experiment))) ** (-0.5) / procs_data_values - inv_tot = (np.diag(la.inv(matrix_theory + matrix_experiment))) ** (-0.5) / procs_data_values - plot_index = theory_covmat_custom.index - df_inv_exp = pd.DataFrame(inv_exp, index=plot_index) - df_inv_exp.sort_index(axis=0, inplace=True) - oldindex = df_inv_exp.index.tolist() - newindex = sorted(oldindex, key=_get_key) - df_inv_exp = df_inv_exp.reindex(newindex) - df_inv_tot = pd.DataFrame(inv_tot, index=plot_index) - df_inv_tot.sort_index(axis=0, inplace=True) - df_inv_tot = df_inv_tot.reindex(newindex) - fig, ax = plotutils.subplots() - ax.plot(df_inv_exp.values, ".", label="Experiment", color="orange") - ax.plot(df_inv_tot.values, ".", label="Experiment + Theory", color="mediumseagreen") - ticklocs, ticklabels, startlocs = matrix_plot_labels(df_inv_exp) - ax.set_xticks(ticklocs) - ax.set_xticklabels(ticklabels, rotation="vertical", fontsize=20) - ax.vlines(startlocs, 0, len(matrix_theory), linestyles="dashed") - ax.set_ylabel(r"$\frac{1}{D_i}\frac{1}{\sqrt{[cov^{-1}_]{ii}}}$", fontsize=30) - ax.yaxis.set_tick_params(labelsize=20) - ax.set_title(f"Diagonal impact of adding theory covariance matrix for {l} points", fontsize=20) - ax.legend(fontsize=20) - ax.margins(x=0) - return fig - - -@figure -def plot_datasets_chi2_theory(procs_data, each_dataset_chi2, abs_chi2_data_theory_dataset): - """Plot the chi² of all datasets, before and after adding theory errors, with bars.""" - ds = iter(each_dataset_chi2) - dstheory = iter(abs_chi2_data_theory_dataset) - dschi2 = [] - dschi2theory = [] - xticks = [] - for proc in procs_data: - for dataset, dsres in zip(proc, ds): - dschi2.append(dsres.central_result / dsres.ndata) - xticks.append(dataset.name) - for proc in procs_data: - for dataset, dsres in zip(proc, dstheory): - dschi2theory.append(dsres.central_result / dsres.ndata) - plotvalues = np.stack((dschi2theory, dschi2)) - fig, ax = plotutils.barplot( - plotvalues, collabels=xticks, datalabels=["experiment + theory", "experiment"] - ) - ax.set_title(r"$\chi^2$ distribution for datasets") - ax.legend(fontsize=14) - return fig diff --git a/validphys2/src/validphys/theorycovariance/tests.py b/validphys2/src/validphys/theorycovariance/tests.py index b5b302ddd1..043d3b3a90 100644 --- a/validphys2/src/validphys/theorycovariance/tests.py +++ b/validphys2/src/validphys/theorycovariance/tests.py @@ -24,7 +24,7 @@ covmap, covs_pt_prescrip, process_starting_points, - theory_corrmat_singleprocess, + theory_corrmat_custom, theory_covmat_custom, ) from validphys.theorycovariance.output import _get_key, matrix_plot_labels @@ -248,7 +248,7 @@ def theory_matrix_threshold(theory_threshold: (int, float) = 0): def theory_corrmat_custom_dataspecs(theory_covmat_custom_dataspecs): """Calculates the theory correlation matrix for scale variations with variations by process type""" - mat = theory_corrmat_singleprocess(theory_covmat_custom_dataspecs) + mat = theory_corrmat_custom(theory_covmat_custom_dataspecs) return mat