From 24093bd4b82240f65b2ca03497dd8804c642c157 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 9 Feb 2022 19:04:04 +0100 Subject: [PATCH 1/7] create a copy of xplotting grid instead of _replace --- validphys2/src/validphys/correlations.py | 2 +- validphys2/src/validphys/deltachi2.py | 4 +-- validphys2/src/validphys/eff_exponents.py | 4 +-- validphys2/src/validphys/pdfgrids.py | 41 ++++++++++++++++++----- 4 files changed, 37 insertions(+), 14 deletions(-) diff --git a/validphys2/src/validphys/correlations.py b/validphys2/src/validphys/correlations.py index 0691a71f96..72dd307cb3 100644 --- a/validphys2/src/validphys/correlations.py +++ b/validphys2/src/validphys/correlations.py @@ -67,7 +67,7 @@ def obs_pdf_correlations(pdf, results, xplotting_grid): values on a grid of (x,f) points in a format similar to `xplotting_grid`.""" _ , th = results corrs = _basic_obs_pdf_correlation(xplotting_grid.grid_values, th._rawdata) - return xplotting_grid._replace(grid_values=corrs) + return xplotting_grid.copy_grid(grid_values=corrs) corrpair_results = collect('results', ['corrpair']) diff --git a/validphys2/src/validphys/deltachi2.py b/validphys2/src/validphys/deltachi2.py index 65f0a9f8a8..503a514df1 100644 --- a/validphys2/src/validphys/deltachi2.py +++ b/validphys2/src/validphys/deltachi2.py @@ -152,8 +152,8 @@ def pos_neg_xplotting_grids(delta_chi2_hessian, xplotting_grid): pos_grid = xplotting_grid.grid_values[pos_mask] neg_grid = xplotting_grid.grid_values[neg_mask] - pos_xplotting_grid = xplotting_grid._replace(grid_values=pos_grid) - neg_xplotting_grid = xplotting_grid._replace(grid_values=neg_grid) + pos_xplotting_grid = xplotting_grid.copy_grid(grid_values=pos_grid) + neg_xplotting_grid = xplotting_grid.copy_grid(grid_values=neg_grid) return [xplotting_grid, pos_xplotting_grid, neg_xplotting_grid] diff --git a/validphys2/src/validphys/eff_exponents.py b/validphys2/src/validphys/eff_exponents.py index a7b07c2ecc..97361e5931 100644 --- a/validphys2/src/validphys/eff_exponents.py +++ b/validphys2/src/validphys/eff_exponents.py @@ -73,7 +73,7 @@ def alpha_eff(pdf: PDF, *, warnings.simplefilter('ignore', RuntimeWarning) alphaGrid_values = -np.log(abs(pdfGrid_values/xGrid))/np.log(xGrid) alphaGrid_values[alphaGrid_values == - np.inf] = np.nan # when PDF_i =0 - alphaGrid = pdfGrid._replace(grid_values=alphaGrid_values) + alphaGrid = pdfGrid.copy_grid(grid_values=alphaGrid_values) return alphaGrid @check_positive('Q') @@ -117,7 +117,7 @@ def beta_eff(pdf, *, warnings.simplefilter('ignore', RuntimeWarning) betaGrid_values = np.log(abs(pdfGrid_values/xGrid))/np.log(1-xGrid) betaGrid_values[betaGrid_values == -np.inf] = np.nan # when PDF_i =0 - betaGrid = pdfGrid._replace(grid_values=betaGrid_values) + betaGrid = pdfGrid.copy_grid(grid_values=betaGrid_values) return betaGrid # .grid_values diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py index 0f9061c943..67fc9f97b2 100644 --- a/validphys2/src/validphys/pdfgrids.py +++ b/validphys2/src/validphys/pdfgrids.py @@ -3,6 +3,7 @@ to facilitate plotting and analysis. """ from collections import namedtuple +from dataclasses import dataclass import numbers import numpy as np @@ -11,7 +12,7 @@ from reportengine import collect from reportengine.checks import make_argcheck, CheckError, check_positive, check -from validphys.core import PDF +from validphys.core import PDF, Stats from validphys.gridvalues import (evaluate_luminosity) from validphys.pdfbases import (Basis, check_basis) from validphys.checks import check_pdf_normalize_to, check_xlimits @@ -37,9 +38,30 @@ def xgrid(xmin:numbers.Real=1e-5, xmax:numbers.Real=1, return (scale, arr) - -XPlottingGrid = namedtuple('XPlottingGrid', ('Q', 'basis', 'flavours', 'xgrid', - 'grid_values', 'scale')) +@dataclass +class XPlottingGrid: + """DataClass holding the value of the PDF at the specified + values of x, Q and flavour. + It also exposes a `stats_gv` attribute with a `Stats` instance + of the raw `grid_values` object in order to compute statistical + estimators in a sensible manner. + """ + Q: float + basis: (str, Basis) + flavours: (list, tuple, type(None)) + xgrid: np.ndarray + grid_values: np.ndarray + scale: str + stats_gv: Stats # TODO: transitional variable during development, it should substitute grid_values by the end of this PR!!!!!!!! + + def copy_grid(self, grid_values=None): + """Create a copy of the grid with potentially a different set of values""" + new_values = {} + if grid_values is not None: + new_values["grid_values"] = grid_values + new_values["stats_gv"] = self.stats_gv.__class__(grid_values) + new_variables = {**self.__dict__, **new_values} + return self.__class__(**new_variables) @make_argcheck(check_basis) @@ -73,8 +95,9 @@ def xplotting_grid(pdf:PDF, Q:(float,int), xgrid=None, basis:(str, Basis)='flavo #Eliminante Q axis #TODO: wrap this in pdf.stats_class? gv = gv.reshape(gv.shape[:-1]) + gv_stats = pdf.stats_class(gv) - res = XPlottingGrid(Q, basis, flavours, xgrid, gv, scale) + res = XPlottingGrid(Q, basis, flavours, xgrid, gv, scale, gv_stats) return res xplotting_grids = collect(xplotting_grid, ('pdfs',)) @@ -242,7 +265,7 @@ def distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(None))=None for grid, pdf in zip(xplotting_grids, pdfs): if pdf == pdfs[normalize_to]: - newgrid = grid._replace(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2]))) + newgrid = grid.copy_grid(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2]))) newgrids.append(newgrid) continue @@ -253,7 +276,7 @@ def distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(None))=None # the distance definition distance = np.sqrt((cv1-cv2)**2/(sg1**2/N1+sg2**2/N2)) - newgrid = grid._replace(grid_values=distance) + newgrid = grid.copy_grid(grid_values=distance) newgrids.append(newgrid) return newgrids @@ -281,7 +304,7 @@ def variance_distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(No for grid, pdf in zip(xplotting_grids, pdfs): if pdf == pdfs[normalize_to]: - newgrid = grid._replace(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2]))) + newgrid = grid.copy_grid(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2]))) newgrids.append(newgrid) continue @@ -293,7 +316,7 @@ def variance_distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(No # the distance definition variance_distance = np.sqrt((sg1**2-sg2**2)**2/(s1+s2)) - newgrid = grid._replace(grid_values=variance_distance) + newgrid = grid.copy_grid(grid_values=variance_distance) newgrids.append(newgrid) return newgrids From cdba34db1bd3c97eafe8b52ef6f29692ad7fdab3 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 10 Feb 2022 12:31:12 +0100 Subject: [PATCH 2/7] continue using stats wherever it can be used --- validphys2/src/validphys/eff_exponents.py | 11 ++++------- validphys2/src/validphys/mc2hessian.py | 6 +++--- validphys2/src/validphys/pdfplots.py | 5 +---- 3 files changed, 8 insertions(+), 14 deletions(-) diff --git a/validphys2/src/validphys/eff_exponents.py b/validphys2/src/validphys/eff_exponents.py index 97361e5931..fcb40aa466 100644 --- a/validphys2/src/validphys/eff_exponents.py +++ b/validphys2/src/validphys/eff_exponents.py @@ -348,17 +348,14 @@ def next_effective_exponents_table( eff_exp_data = [] - alphagrid = alpha_effs.grid_values - betagrid = beta_effs.grid_values - - alphastats = pdf.stats_class(alphagrid) - betastats = pdf.stats_class(betagrid) + alphastats = alpha_effs.stats_gv + betastats = beta_effs.stats_gv with warnings.catch_warnings(): warnings.simplefilter('ignore', RuntimeWarning) - alpha_cv = np.nanmean(alphagrid, axis=0) - beta_cv = np.nanmean(betagrid, axis=0) + alpha_cv = np.nanmean(alphastats.error_members(), axis=0) + beta_cv = np.nanmean(betastats.error_members(), axis=0) #tuple of low and high values repectively alpha68 = alphastats.errorbar68() beta68 = betastats.errorbar68() diff --git a/validphys2/src/validphys/mc2hessian.py b/validphys2/src/validphys/mc2hessian.py index 4e045ced64..33eb2aab69 100644 --- a/validphys2/src/validphys/mc2hessian.py +++ b/validphys2/src/validphys/mc2hessian.py @@ -108,9 +108,9 @@ def _create_mc2hessian(pdf, Q, xgrid, Neig, output_path, name=None): def _get_X(pdf, Q, xgrid, reshape=False): pdf_grid = xplotting_grid(pdf, Q, xgrid=xgrid) - pdf_grid_values = pdf_grid.grid_values - replicas = pdf_grid_values - mean = pdf_grid_values.mean(axis=0) + pdf_grid_values = pdf_grid.stats_gv + replicas = pdf_grid_values.error_members() + mean = pdf_grid_values.central_value() Xt = replicas - mean if reshape: Xt = Xt.reshape(Xt.shape[0], Xt.shape[1] * Xt.shape[2]) diff --git a/validphys2/src/validphys/pdfplots.py b/validphys2/src/validphys/pdfplots.py index 74397e2b83..fc5571b812 100644 --- a/validphys2/src/validphys/pdfplots.py +++ b/validphys2/src/validphys/pdfplots.py @@ -77,10 +77,7 @@ def fp_error(tp, flag): np.seterrcall(fp_error) for grid in self._xplotting_grids: newvalues = grid.grid_values/normvals - #newgrid is like the old grid but with updated values - newgrid = type(grid)(**{**grid._asdict(), - 'grid_values':newvalues}) - newgrids.append(newgrid) + newgrids.append(grid.copy_grid(grid_values=newvalues)) return newgrids return self._xplotting_grids From da3069c7ade11551032396c8c99d48f3631c53ce Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 10 Feb 2022 16:32:00 +0100 Subject: [PATCH 3/7] dont use internal _rawdata, move thpredictions to the StatResults class --- validphys2/src/validphys/calcutils.py | 4 +-- .../validphys/closuretest/closure_results.py | 6 ++-- .../src/validphys/closuretest/multiclosure.py | 10 +++--- .../closuretest/multiclosure_pseudodata.py | 2 +- validphys2/src/validphys/correlations.py | 6 ++-- validphys2/src/validphys/covmats.py | 2 +- validphys2/src/validphys/mc2hessian.py | 3 ++ validphys2/src/validphys/pdfgrids.py | 12 +++---- validphys2/src/validphys/results.py | 33 ++++++++++++------- .../src/validphys/tests/test_closuretest.py | 2 +- .../src/validphys/tests/test_pyfkdata.py | 6 ++-- 11 files changed, 47 insertions(+), 39 deletions(-) diff --git a/validphys2/src/validphys/calcutils.py b/validphys2/src/validphys/calcutils.py index af4e721735..2601e084d6 100644 --- a/validphys2/src/validphys/calcutils.py +++ b/validphys2/src/validphys/calcutils.py @@ -70,7 +70,7 @@ def all_chi2(results): """Return the chi² for all elements in the result. Note that the interpretation of the result will depend on the PDF error type""" data_result, th_result = results - diffs = th_result._rawdata - data_result.central_value[:,np.newaxis] + diffs = th_result.rawdata - data_result.central_value[:,np.newaxis] return calc_chi2(sqrtcov=data_result.sqrtcovmat, diffs=diffs) def central_chi2(results): @@ -85,7 +85,7 @@ def all_chi2_theory(results, totcov): """Like all_chi2 but here the chi² are calculated using a covariance matrix that is the sum of the experimental covmat and the theory covmat.""" data_result, th_result = results - diffs = th_result._rawdata - data_result.central_value[:,np.newaxis] + diffs = th_result.rawdata - data_result.central_value[:,np.newaxis] total_covmat = np.array(totcov) return calc_chi2(sqrtcov=la.cholesky(total_covmat, lower=True), diffs=diffs) diff --git a/validphys2/src/validphys/closuretest/closure_results.py b/validphys2/src/validphys/closuretest/closure_results.py index 9f30cf15ee..857720fe4a 100644 --- a/validphys2/src/validphys/closuretest/closure_results.py +++ b/validphys2/src/validphys/closuretest/closure_results.py @@ -122,7 +122,7 @@ def bootstrap_bias_experiment( """ dt_ct, th_ct = dataset_inputs_results ((_, th_ul),) = underlying_dataset_inputs_results - th_ct_boot_cv = bootstrap_values(th_ct._rawdata, bootstrap_samples) + th_ct_boot_cv = bootstrap_values(th_ct.rawdata, bootstrap_samples) boot_diffs = th_ct_boot_cv - th_ul.central_value[:, np.newaxis] boot_bias = calc_chi2(dt_ct.sqrtcovmat, boot_diffs) / len(dt_ct) return boot_bias @@ -191,7 +191,7 @@ def variance_dataset(results, fit, use_fitcommondata): """ dt, th = results - diff = th.central_value[:, np.newaxis] - th._rawdata + diff = th.central_value[:, np.newaxis] - th.rawdata var_unnorm = calc_chi2(dt.sqrtcovmat, diff).mean() return VarianceData(var_unnorm, len(th)) @@ -210,7 +210,7 @@ def bootstrap_variance_experiment(dataset_inputs_results, bootstrap_samples=500) normalised to the number of data in the experiment. """ dt_ct, th_ct = dataset_inputs_results - diff = th_ct.central_value[:, np.newaxis] - th_ct._rawdata + diff = th_ct.central_value[:, np.newaxis] - th_ct.rawdata var_unnorm_boot = bootstrap_values( diff, bootstrap_samples, diff --git a/validphys2/src/validphys/closuretest/multiclosure.py b/validphys2/src/validphys/closuretest/multiclosure.py index c784c8996a..c48388b238 100644 --- a/validphys2/src/validphys/closuretest/multiclosure.py +++ b/validphys2/src/validphys/closuretest/multiclosure.py @@ -130,7 +130,7 @@ def fits_dataset_bias_variance( """ closures_th, law_th, _, sqrtcov = internal_multiclosure_dataset_loader # The dimentions here are (fit, data point, replica) - reps = np.asarray([th._rawdata[:, :_internal_max_reps] for th in closures_th]) + reps = np.asarray([th.rawdata[:, :_internal_max_reps] for th in closures_th]) # take mean across replicas - since we might have changed no. of reps centrals = reps.mean(axis=2) # place bins on first axis @@ -226,7 +226,7 @@ def dataset_xi(internal_multiclosure_dataset_loader): """ closures_th, law_th, covmat, _ = internal_multiclosure_dataset_loader - replicas = np.asarray([th._rawdata for th in closures_th]) + replicas = np.asarray([th.rawdata for th in closures_th]) centrals = np.mean(replicas, axis=-1) underlying = law_th.central_value @@ -297,7 +297,7 @@ class BootstrappedTheoryResult: """ def __init__(self, data): - self._rawdata = data + self.rawdata = data self.central_value = data.mean(axis=1) @@ -341,7 +341,7 @@ def _bootstrap_multiclosure_fits( # construct proxy fits theory predictions for fit_th in fit_boot_th: rep_boot_index = rng.choice(n_rep_max, size=n_rep, replace=use_repeats) - boot_ths.append(BootstrappedTheoryResult(fit_th._rawdata[:, rep_boot_index])) + boot_ths.append(BootstrappedTheoryResult(fit_th.rawdata[:, rep_boot_index])) return (boot_ths, *input_tuple) @@ -741,7 +741,7 @@ def dataset_fits_bias_replicas_variance_samples( """ closures_th, law_th, _, sqrtcov = internal_multiclosure_dataset_loader # The dimentions here are (fit, data point, replica) - reps = np.asarray([th._rawdata[:, :_internal_max_reps] for th in closures_th]) + reps = np.asarray([th.rawdata[:, :_internal_max_reps] for th in closures_th]) # take mean across replicas - since we might have changed no. of reps centrals = reps.mean(axis=2) # place bins on first axis diff --git a/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py b/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py index ff32d86fd1..22b2c78a1a 100644 --- a/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py +++ b/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py @@ -126,7 +126,7 @@ def experiments_closure_pseudodata_estimators_table( fits_erep_delta_chi2 = [] for i_fit, fit_th in enumerate(closures_th): # some of these could be done outside of loop, but easier to do here. - th_replicas = fit_th._rawdata + th_replicas = fit_th.rawdata th_central = np.mean(th_replicas, axis=-1) dt_replicas = fits_reps_pseudo[i_fit].xs(exp_name, axis=0, level=0).to_numpy() dt_central = fits_cv.xs(exp_name, axis=0, level=0).iloc[:, i_fit].to_numpy() diff --git a/validphys2/src/validphys/correlations.py b/validphys2/src/validphys/correlations.py index 72dd307cb3..7a254467e5 100644 --- a/validphys2/src/validphys/correlations.py +++ b/validphys2/src/validphys/correlations.py @@ -60,13 +60,13 @@ def _basic_obs_obs_correlation(obsarr1, obsarr2): return x@y/np.outer(la.norm(x,axis=1),la.norm(y,axis=0)) -#TODO: Implement for other error types. Do not use the _rawdata. +#TODO: Implement for other error types. Do not use the rawdata. @check_pdf_is_montecarlo def obs_pdf_correlations(pdf, results, xplotting_grid): """Return the correlations between each point in a dataset and the PDF values on a grid of (x,f) points in a format similar to `xplotting_grid`.""" _ , th = results - corrs = _basic_obs_pdf_correlation(xplotting_grid.grid_values, th._rawdata) + corrs = _basic_obs_pdf_correlation(xplotting_grid.grid_values, th.rawdata) return xplotting_grid.copy_grid(grid_values=corrs) @@ -77,4 +77,4 @@ def obs_pdf_correlations(pdf, results, xplotting_grid): def obs_obs_correlations(pdf, corrpair_results): """Return the theoretical correlation matrix between a pair of observables.""" (_,th1), (_,th2) = corrpair_results - return _basic_obs_obs_correlation(th1._rawdata, th2._rawdata) + return _basic_obs_obs_correlation(th1.rawdata, th2.rawdata) diff --git a/validphys2/src/validphys/covmats.py b/validphys2/src/validphys/covmats.py index bd09c15d5b..d4d389ea6b 100644 --- a/validphys2/src/validphys/covmats.py +++ b/validphys2/src/validphys/covmats.py @@ -531,7 +531,7 @@ def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered): True """ th = ThPredictionsResult.from_convolution(pdf, dataset) - pdf_cov = np.cov(th._rawdata, rowvar=True) + pdf_cov = np.cov(th.rawdata, rowvar=True) return pdf_cov + covmat_t0_considered diff --git a/validphys2/src/validphys/mc2hessian.py b/validphys2/src/validphys/mc2hessian.py index 33eb2aab69..b8d1aeef02 100644 --- a/validphys2/src/validphys/mc2hessian.py +++ b/validphys2/src/validphys/mc2hessian.py @@ -17,6 +17,8 @@ from validphys.lhio import hessian_from_lincomb from validphys.pdfgrids import xplotting_grid +from validphys.checks import check_pdf_is_montecarlo + log = logging.getLogger(__name__) @@ -58,6 +60,7 @@ def mc2hessian_xgrid( ) +@check_pdf_is_montecarlo def mc2hessian( pdf, Q, Neig: int, mc2hessian_xgrid, output_path, gridname, installgrid: bool = False ): diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py index 67fc9f97b2..7c9f677641 100644 --- a/validphys2/src/validphys/pdfgrids.py +++ b/validphys2/src/validphys/pdfgrids.py @@ -3,7 +3,7 @@ to facilitate plotting and analysis. """ from collections import namedtuple -from dataclasses import dataclass +import dataclasses import numbers import numpy as np @@ -38,7 +38,7 @@ def xgrid(xmin:numbers.Real=1e-5, xmax:numbers.Real=1, return (scale, arr) -@dataclass +@dataclasses.dataclass class XPlottingGrid: """DataClass holding the value of the PDF at the specified values of x, Q and flavour. @@ -56,12 +56,8 @@ class XPlottingGrid: def copy_grid(self, grid_values=None): """Create a copy of the grid with potentially a different set of values""" - new_values = {} - if grid_values is not None: - new_values["grid_values"] = grid_values - new_values["stats_gv"] = self.stats_gv.__class__(grid_values) - new_variables = {**self.__dict__, **new_values} - return self.__class__(**new_variables) + new_stats_gv = self.stats_gv.__class__(grid_values) + return dataclasses.replace(self, grid_values=grid_values, stats_gv=new_stats_gv) @make_argcheck(check_basis) diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index d989da6caa..9d03014566 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -73,6 +73,10 @@ class StatsResult(Result): def __init__(self, stats): self.stats = stats + @property + def rawdata(self): + return self.stats.data.T + @property def central_value(self): return self.stats.central_value() @@ -106,7 +110,7 @@ def sqrtcovmat(self): return self._sqrtcovmat -class ThPredictionsResult(NNPDFDataResult): +class ThPredictionsResult(StatsResult): """Class holding theory prediction For legacy purposes it still accepts libNNPDF datatypes, but prefers python-pure stuff """ @@ -122,11 +126,20 @@ def __init__(self, dataobj, stats_class, label=None, central_value=None): except AttributeError: self._std_error = dataobj.get_error() self._rawdata = dataobj.get_data() - super().__init__(dataobj, central_value=central_value) + statsobj = stats_class(dataobj.T) + super().__init__(statsobj) + self._central = central_value + + # TODO: transitional method so that the test pass by themselves for the next commit before + # offloading this to the stats class (and updating the regressions) @property - def std_error(self): - return self._std_error + def central_value(self): + if hasattr(self, '_central') and self._central is not None: + return np.array(self._central).reshape(-1) + return self.stats.central_value() + + @staticmethod def make_label(pdf, dataset): @@ -175,10 +188,6 @@ def from_convolution(cls, pdf, posset): stats = pdf.stats_class(data.T) return cls(stats) - @property - def rawdata(self): - return self.stats.data - # TODO: finish deprecating all dependencies on this index largely in theorycovmat module groups_data = collect("data", ("group_dataset_inputs_by_metadata",)) @@ -257,7 +266,7 @@ def group_result_table_no_table(groups_results, groups_index): ): replicas = ( ("rep_%05d" % (i + 1), th_rep) - for i, th_rep in enumerate(th._rawdata[index, :]) + for i, th_rep in enumerate(th.rawdata[index, :]) ) result_records.append( @@ -612,7 +621,7 @@ def dataset_inputs_bootstrap_phi_data(dataset_inputs_results, bootstrap_samples= For more information on how phi is calculated see `phi_data` """ dt, th = dataset_inputs_results - diff = np.array(th._rawdata - dt.central_value[:, np.newaxis]) + diff = np.array(th.rawdata - dt.central_value[:, np.newaxis]) phi_resample = bootstrap_values( diff, bootstrap_samples, @@ -632,7 +641,7 @@ def dataset_inputs_bootstrap_chi2_central( a different value can be specified in the runcard. """ dt, th = dataset_inputs_results - diff = np.array(th._rawdata - dt.central_value[:, np.newaxis]) + diff = np.array(th.rawdata - dt.central_value[:, np.newaxis]) cchi2 = lambda x, y: calc_chi2(y, x.mean(axis=1)) chi2_central_resample = bootstrap_values( diff, @@ -727,7 +736,7 @@ def positivity_predictions_data_result(pdf, posdataset): def count_negative_points(possets_predictions): """Return the number of replicas with negative predictions for each bin in the positivity observable.""" - return np.sum([(r.rawdata < 0).sum(axis=1) for r in possets_predictions], axis=0) + return np.sum([(r.rawdata < 0).sum(axis=0) for r in possets_predictions], axis=0) chi2_stat_labels = { diff --git a/validphys2/src/validphys/tests/test_closuretest.py b/validphys2/src/validphys/tests/test_closuretest.py index a55eb4063a..68fc91dc6e 100644 --- a/validphys2/src/validphys/tests/test_closuretest.py +++ b/validphys2/src/validphys/tests/test_closuretest.py @@ -13,7 +13,7 @@ class TestResult: """class for testing base level estimators which expect a results object""" def __init__(self, central_value, rawdata=None): self.central_value = central_value - self._rawdata = rawdata + self.rawdata = rawdata self.ndata = len(central_value) self.sqrtcovmat = np.identity(self.ndata) diff --git a/validphys2/src/validphys/tests/test_pyfkdata.py b/validphys2/src/validphys/tests/test_pyfkdata.py index 0f6b90b6a2..37c45720c9 100644 --- a/validphys2/src/validphys/tests/test_pyfkdata.py +++ b/validphys2/src/validphys/tests/test_pyfkdata.py @@ -68,7 +68,7 @@ def test_predictions(pdf_name): ds = l.check_dataset(**daset, theoryid=THEORYID) preds = predictions(ds, pdf) core_predictions = ThPredictionsResult.from_convolution(pdf, ds) - assert_allclose(preds.values, core_predictions._rawdata, atol=1e-8, rtol=1e-3) + assert_allclose(preds.values, core_predictions.rawdata, atol=1e-8, rtol=1e-3) # Now check that the stats class does the right thing with the data cv_predictions = central_predictions(ds, pdf).squeeze() stats_predictions = pdf.stats_class(preds.T) @@ -90,12 +90,12 @@ def test_positivity(pdf_name): ps = l.check_posset(setname=posset, theoryID=THEORYID, postlambda=1e6) preds = predictions(ps, pdf) core_predictions = PositivityResult.from_convolution(pdf, ps) - assert_allclose(preds.values, core_predictions.rawdata.T) + assert_allclose(preds.values, core_predictions.rawdata) # Now do the same with the API api_predictions = API.positivity_predictions_data_result( theoryid=THEORYID, pdf=pdf_name, posdataset={"dataset": posset, "maxlambda": 1e6} ) - assert_allclose(preds.values, api_predictions.rawdata.T) + assert_allclose(preds.values, api_predictions.rawdata) # And now check that the results are correct for any kind of PDF cv_predictions = central_predictions(ps, pdf).squeeze() assert_allclose(cv_predictions, api_predictions.central_value, atol=1e-3) From 0c70f5c9205e073653eb3637d29cfd9480b80ac2 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 10 Feb 2022 17:16:31 +0100 Subject: [PATCH 4/7] update regression tests to the StatsResult predictions --- validphys2/src/validphys/results.py | 36 +--- .../tests/baseline/test_dataspecschi2.png | Bin 23215 -> 23213 bytes .../tests/regressions/test_datasetchi2.csv | 6 +- .../regressions/test_replicachi2data.csv | 202 +++++++++--------- 4 files changed, 108 insertions(+), 136 deletions(-) diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index 9d03014566..052eb1f6bf 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -111,35 +111,13 @@ def sqrtcovmat(self): class ThPredictionsResult(StatsResult): - """Class holding theory prediction - For legacy purposes it still accepts libNNPDF datatypes, but prefers python-pure stuff + """Class holding theory prediction, inherits from StatsResult """ - def __init__(self, dataobj, stats_class, label=None, central_value=None): + def __init__(self, dataobj, stats_class, label=None): self.stats_class = stats_class self.label = label - # Ducktype the input into numpy arrays - try: - self._rawdata = dataobj.to_numpy() - # If the numpy conversion worked then we don't have a libNNPDF in our hands - stats = stats_class(self._rawdata.T) - self._std_error = stats.std_error() - except AttributeError: - self._std_error = dataobj.get_error() - self._rawdata = dataobj.get_data() - statsobj = stats_class(dataobj.T) super().__init__(statsobj) - self._central = central_value - - # TODO: transitional method so that the test pass by themselves for the next commit before - # offloading this to the stats class (and updating the regressions) - @property - def central_value(self): - if hasattr(self, '_central') and self._central is not None: - return np.array(self._central).reshape(-1) - return self.stats.central_value() - - @staticmethod def make_label(pdf, dataset): @@ -165,20 +143,14 @@ def from_convolution(cls, pdf, dataset): datasets = (dataset,) try: - all_preds = [] - all_centrals = [] - for d in datasets: - all_preds.append(predictions(d, pdf)) - all_centrals.append(central_predictions(d, pdf)) + th_predictions = pd.concat([predictions(d, pdf) for d in datasets]) except PredictionsRequireCutsError as e: raise PredictionsRequireCutsError("Predictions from FKTables always require cuts, " "if you want to use the fktable intrinsic cuts set `use_cuts: 'internal'`") from e - th_predictions = pd.concat(all_preds) - central_values = pd.concat(all_centrals) label = cls.make_label(pdf, dataset) - return cls(th_predictions, pdf.stats_class, label, central_value=central_values) + return cls(th_predictions, pdf.stats_class, label) class PositivityResult(StatsResult): diff --git a/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png b/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png index 7472b78f635b44f1ad57150a23ef4b4b9f843581..0115b2daa34ff91bd168412dc1607937cb407fbd 100644 GIT binary patch delta 13668 zcmZv@cQ}^)8$W)b_2fZY$!KVqPfsb4(bTQ%6_J)?Zz5cpDn(g6WoF*oRz@~GNs@;O zH)W3`goI@OUhh8N-*J5Z_;wr}N8=vX^}f#Yb-uuhgr=WK0pPRh>FT>9hLAH@_!{b7I1(21h>&XV68Or=B%iZT*-MklSI@xN;i z`B2I-Crdo@7x4-**GL+@xtP`0Y*#v4oxWqojt^Py6j>P=f4O^ji0$8RT$iZfrf90K zpHMo}BW=k}4k;=ssxR@nB6%{S<8N<%8NOk8s(yd?wS^SB<8x(@WqV2SEk2Pn`Qf50 zRnuJ;-YgT6b3P_Q&peBXiE-|DL38XAoS)zHqBn{@JwMZ{G~VeelI$|jFf~#k;xav$ zYc<}_?*I6#yJt$@a#fof#dNn0H8V^T=AX@%;kMme=HNBnl1*9hwA>gHXZdgAx{Q#Nm)IyiZK5A6w4 zi}L#UIIW;xtkrY0`pMsBcP~~y>6zI*mzvwA?d+Uwlxr)=$H#Yn-L9i$w|IiW!ovI# z)AKL4bG#T|vW&978Xdi1;o`-tvt_qxD!RJ6d+cR{LqfJNENkQL7?;fcJQcCuU3=ZS zb*c=8*7Zd^tClbCOnzpyAE#}vz`1}jj=nN6)Kp5BE?v4bJRm+U?u~7(xTt7c?eo*8 zPcP-=)pvK#|MK7C*k;?p>ao0>MJRwgVGd4;9UZhpNS@wkXJ^N#EpQni_t)#fy|>*n zOH52`w{44TcYAxtUEHdrFuo@FU+L1B7=IyoxyCzdNi*kfvLdRWHum%WTt!-?7HiBjWs)%#^v_zWs5DY9Rw`Uc*t>^}c)RWJynO>uBP6(&*Z(GB)Ye!#w zDQEe?Gn5avcB1sb`Juw$`W!0*wa0%S*|>4z-8*+a=x8X!qvx}1?F>VKF>#+iUl{-P zsqTAr4)?GEMJc!HviyCFnc3O9Q;qYrQCd%Snfp6BI_|#s-&!FB4_yL0wr(@O92PF7 z7qjBx4%^z=KG|a<5cR@zQ`;4a@+nACUwpf)^Fj5+I`b&6gj+o8RRa@FJdW9T=+@Wb ziMAj#sc6!H>x)!RoX|aR;3~O>XXh)#{{H*NJ)IkxE=3cl1{a+O-Ut6{7kQc^H15~0 z%lia*@wj%#$r%M+$ZcM|di53t-T(gS^}xWu-4?Z?-hY%V`2XKpvMBEV*(#sIP0bBM z@A~tNZnR831+hkGJAAHtLc4lXjkey)?2c7oFuK7W*yWiEJ*9vo49(5#8T{u>K z>7Red3Xs0=w?2r_N^9) zQoYMcd-Td~D))kdbnC_)@7Y=PJNIdigrDpA{CU|1d=!856N7^REB!|P33l!`{qyR% zjpUx{ml)sR&APgAW@%Q&hJQH~x=$BNKhiCh_HBRQ-EX=tEh0v}OQDo=WVAxUdoFLm z@WkrdQ3_x1*Dc8-v(F7f&Gg!<@BpY6`3_6h1aWMVpz|(lQxCTvUqbl^T%$NB2y?u*rrVfQ8V$QlXdFh zvQP2Ha=q9ti6svNzVtjY)yLc$Dw(#$WQ%YK8*I)U$nTA;>n!u-(O4WUWmM{2j6VN7 zJS@#{baD!+zUZWgmKWrueI=?kouG96Q1HcLwa4!U1T;>x)9KPBzYacq`cy!{BYTa+ zNzFd<_b+d9%b@G|)&1sQ?>H&*23IuJ9iY@78=$Bid-Q%yzDe_Q8yWi!NYK;Gmjx!4vEqdG7Xs{3htHMugFQ!rdIJ@z#=z;g0O4|+1tx&KgL&^C<%N8=fO&bURBxMg%cJO22iQCkBe zqvW~SX{(X0@&!wLs%%l~O@(ghukWwl#<2XD@zm(Kb<&$R$8s#$s>_!zU&_s`Qxz(? zdhJ?fM+bf8+JYBbK4qM}F6TD%?cSeNne~}w$BV}sP62j(UI&PzlkV;5>1ldy6Xx;r zTVL0m5vv8qcrtC8_Z&HV_)T5iXi;;1VIi60O(kB1Pb;ZU`fe*N^*skjEK(J$2? z0}*>{n&L1Cz7`Ib1c`pK-jj3b!$BTeVE6ml_~Yzs3mtw9-sr93Vq(e8-F|fpoa|P2wg*jf+t~|o~KmMQR{wh-L?{s>7ws1*S%Bt)p31lvNwu6b#kqELj6rSEbfF?`S!~_cxuy4mGF2i`fB3? z=jg$xl8n(YlV`3&ElI(_t4!QlUD5b=uI*cI0QXo%n_g%BaXv)#Jz632N*4 zPuz8SX4n!K_j7)6;8+E~ErB!kZ41NC3g>pFZnP~N+=GJK0da$4u&hhCkK$#= z9O0I+{-}w`Q}Y)GwQ`O0>6AyOJs1-$P62IR6DICgu2kXItIt-cD62VO-xcUq_nj3<;$1#1uj-wL`AnNDxLH1>6ej)={`l+HxwXJC=JAYw?WE>DfZncExG0qO{6JdMQ0`~u)Z-ay<%*=S2UA&l# zMiv$pwrnfP2UZ6^co6;Z<4M9Lz!#gIo%dx`yrw(tt_egJ6&c5>1T?hdJFzp)9wL2; z0;~lf2L}fe3TM+^V#2*f;yswu=ic5#RJK6?Pju+1e8+fPUEe)EMu<93j>K_12U(m7H z+J-|~$vAyEF*!9w9y;b?qhxMN%)i$c7yfQg`p(L2tAvD70nEn&oq9&Ic8H19;))H$ zF%|{@0Y~o(#DD*8&v2Y;dDcAmQ?n|<{}<}V{#Ky2Yq>kv>A*$sBUw|{kqH1 zlpn20moVhU8Rs=MA1MCn;JPm$qe_G+22Uc&iU5aH{_M}@9`xdRyCN+S&+k`%4z;qn zy3U1!hN{wXE~>8px^AlXn^V97S4>_r=|6znc^45Xv;2^f;*fi^tLjf1T z1QsJTG(V~|ADg6;n7M-9>2v%?!zsHR+qS8!+@^j(GF!l(%73|K>sCwfX)$qe0ya>6 zphdv2bgT4Yeb5rZmWXEV?d|P%;HY0d(2%Nw-=rIaWbKtjyv81S&y1+V$Hyn3gE8%N zt*uicJ%2EG*8P2>E*+Dut@sKH!_kB8jc^MOfY?`63)`ay+0MPK?B*urH2d`8pHH@) z_$T`PdrcD;wq5~_(4|T;ciDQ%H~~)?=$uMMw$q-7^xeC4Q8R2%r=k84dzvGUXL$9{i>q9LakaJDHXkz% z;701=7(wA~mv7&`t*NUU3@`|1#c;jo$p-1zvTd8eU2Y8x4X}^xn>M|?8N|U%1#J5Z z?u*IG%Tsl0sj95hGBys!jg)UtPCFB7)>pl^9T=|=kLtI=r#Iz0CAto`rGgb)qCaJt zLA*aoOiV<-ZmP=MZPg&LW{T6^knnGJ>`O`{6-^d!sN5J#;CTW2&lm zQNQD_*V~b&1k|>q8&Vz%$XLcrjSp1S)dj_Fppw~<3Rj$ed}g$|wR`=Tm@q);;OR=^ z$h~IA8wqyOby!PO=?2_A<@@#RYb{+sDO`lV(-sVUnc z7NQkb^#`M&JAl=|YDhiB(qrj^Txfw=uM!Y= z`(ZitNEF~mz!wKc6l!(YUaR&JFJg{>S%IyjLn4!}zz7ad`1wzPb6;2gGY8%UZ&s(F zV5cXB33<@g*3PzSltO1%;LUQ5^yk~(Tnovgq36Bt3o9567f(d8R2FIaf%EPF^*a7Y zO3!bNB)qA3!X*f`Mgf3tb?=vMe)moT;}N`9e(xJfbMWrnyTM^$Eb=-$VIt|ZwF7Z8 ziBaBD40KB|OXbx*`=t%5YpSc&2q&!yr6t*Te5&tKPZt~MB&xIF3~pdR1h|ynrurX1 zluS-fi*MX`0Pv}gR`8e@O1*n`X-IJJcAUn(ef!?Nd83IhsjaKK!$&=e!;EgH=ZhEp zxmL5UUJNEzuS+J;*ovn$H8+Teh#1_vD+XHRQ#VoIDl7Q5riN|+9oUTl0mSM5HPTH= z#?!NCD6gB>AYdtGCMX$XhLK8@r%*@P$=Z1H;HBu+m%2?9D zX!{sURrHm9pFA$&^9Q`XC;OZdp4+x21ul%{FhLhcv;ynK1E_*a)7@nB&zMq{)Q0q60fa@%-g=^Q#KwKmhN ziXa&<8{S?=- z4F*i5nNB{{W5+aIU7tJUj6#}=@7-$v9~jDODIcE!>YWWzeqV53VBr~#ycB{5L^%@W zLI^ygHpMV&4E}i)%7{z^3^a7R4DCI0xz)amhwx#)#6_Ax#ED1i5!<{zSx+3Z_B8AW zebA;SGFDp|DZ6yJp7Fj}T?6U&=n?q-u?l}79ew>!Fw}SmFCzXF6imS6*q{j*Tn4{| z$j4pp_QxLaThnJcKlKO(omgMJRujQ1T50ZAEWj@jB57c1Al)r8vD#%LB~@PHH9dwI zFOHH;N=gDi)PgBt2IxX-f2j@=-@5f^gq(}^^v^*yEX5-#Duo`_GmALmZe>aw2fB!D z-%bYp>Xj?Yw!z-(oQN>ToPc(j+`Mb!9JYbdFEQcNiTEpf-~?d{`tDdJ0M=dbTKnwgu6@7TdW zXGw*x&?{fQjCx_`C@V;=QKr#27Dthcl1Ec(U$_tlA7h-JpP*|B&vGcc7(TatOq>m1 z^kdH=i&{Mj^wzkdy9 zo@l7qzSvm6eU0CHd6)h;ULm>dIO3~UuS!L(;J*%MYG!LyCOmPZ^mN}j*3dh*^|<2&b2c*b0!+BxTk-F zs$T(CDCaT8F* z4kefy(xDKPsN}?N->fkSG*BZl`!`B0A}QTvU_7%*h6|S` zn+_V+N7THdBk(hGtV+6cyNnEb(i3rwbJwdsVAP!CTFhI|HZV>7=N}kne@5>)nRW5L z$Jp1WIF{LW^D_ZY=sh;cfPE`)8>0T<&AC*|G=9Z2(1u}yDqpxknWB2fKd6XUKxz?R z5pY(YVz>(}$I^&Wv}{a^xcDwIRYPf37Ya8JWyY$2)V^S5M&5uEs)KXa)7`DpXY>H0 zlbHYX+-THzh0>f>T>`D<)R~>9nR0^{H-tuGB-WIb0DH!ImF5kHhK2|?gdaf+M+7dt zeSE!zmsU~9=ivH~@a)x{m4yDoxWPyFjzY>1&4sHRLuA-5O;2)LU8~;^F0mM&PJ+KI z^LP4>q)f(Ma!X(O-tli!4sG?QrB%c0R*r3}97c~im~-G^Pb%+%((G6e8pLHVr%_JC z3BTr?rWQnc&rT(Q*cm~|vW}J?`F%D(lXwy*M7Xtmz@yHAnL;E7#Pw}b<^V;P{`$Uf zw>L0)$k(_GwHQL}20MRQygrehZ-J@6PB%4$5Pfa_s(dq!Tdy|X$$WNrelERaZpw<`LvV?5Br(+goayn@>G{{1{`P^xLkdQ?Rdn#2jP=0fZQi zbqEp1;G(g>!7wi1@TWMNOHg?yG1bI&?4UnDs1c6)fAl}EZ}&Z$0N?~HW*8B#7Y22m zbKsW%Ajx@+^!BZOZ?h*1z75_a<2B{%H8Da&S=LjqoZ+rvfhezJ%GiMp}6xb2YngfhIhi5M6LsP-{c(1l@6vvLAyHA%poVY=1fKA+uaXbd z15>sw8aawAOouovh^@#lWFbhJLSgT|^nL>pC`pq7U4#dkmX=~}Wt{d5Sa(Qyp^Zoh z648hmQgc&%{uFKorrn;b>Fen`qCCsi3zsjE-Da!G9WVli>U40 zm%pDGY%4xRtck;KKqTP*ZsJ$Fd;EY8Mb-S~XZ1-z;F;T?ih9oz#sroN6*^inqv9_J zTv0RGC*@;w1{9R6BM|;cYL%1=@yCF8ga?(DIMN$1cp(_vfI=sKUpyHccX32!_kf z&_-pAaxJglpTY&uK7?mM#;`?Qf(4GtIBT(B9XMj zWVfW$LK4o`*w?A;-SO97n$UE8cD7cAp=et&1GtD0LGoMNk$QURT*3bR`@;iHg~BBu z{w9y0tmd6N#Ez{4qsTsA!OO3H&##<5@|E)06$hgR(^Rji;bPWu5lZS4pxxCS05-s&SkUE0fmwY!Kj4St4hBR6y8$W7TZyd92U)N|_+ zO&Nh@otniSaxgF++I4UM7?zRreote(rMIp{A$%O?=5#-t_Ve{|%>4@ZTOS1XCr zhM3(McB<8n&y}QyxD(6cyQEH1KB`!CkTNwjB@GNi`iZ<-8c~8EQB`GShv2kn!m&V_ zm4copIZVqnek5Gk7C~ez|GQX3Rzw^l18YFzZC8eA!u^PQoZIfrA~8RlFsNLq@U@Rf z9IBfM9!ZjJSiKGK^N0-qmH^X5e!lX}oA;PQAuxg+x?fw=#xc@_*DLcpcka_h!j}xh zt*O{!+rlJej=!);{*M(I3%DU3Umlxh15$VdY1#KTS9WWYnW;*`d*`p3WT#0~^f26+ zf-g_3!r#H<8+iZp5D^4Gf_k{i{Oa7%dp;svCzc?oQ_NOC6#}4iLcA2_es$OQ)!6}b z)S?tk7z_q7mZPZMa3B;;b(M&SY2Qu;M8WI~JwqIFnBC(;t<6V*IGJj|ZmnjEwMJv7`@)wyI2eiQTsh73psJc08wO1>2iKy!%SOgx z*!t0QTo^;5@=#E&EnVh&a}0B z?As_ASxj{@)C1KcF7^z0e=%K2-p{_HF@FxJz7>pNG%7K>U_^W#pa;Y;LKw+oyDmbr z5!L%EUxzzmg&d_ENK80PDrMI~6Yw=%Fb^C}qI!fJoMUwa)du)16NRllMnLO5gfsva zSE1REnqZb;8|pCrorE=@KpXh@0nf^~oNk{}^%!v}Hf;KdOn}|^&B`qyYY+fYG*r?Hhze=zj z2OhC3sS^XhN3;)+U1wnG)w5ry3z}CB8F{?k`w5uvUeol?9sBkfl2`)V1uV=*ASEyF zBpC{Glkezc2ttp~4^>DO1VTuxQ+d~K#P9+oh{XjX?MM9kf_=Vw1)z{()DF}0SMZO0 z*kan(E9t6%MG-hB;&|@kNHekz%2CGm0s~iIQ|O>w@p$MS8lf~H&+q8DEGLB+SSMuf zj06pE9MXz^8u4q~JAMScBg#OsP+O;=l!`vquT3gnq%j}nqC^?n!c*~>pPi6#|8W6O z_Z3QonA>T9hC+|sy1RywmQW8wTb2&3-@esFKswu} zL2H9_5s%&rpeL~N+>KBP!x+G!LBJi$SEMm?kZfY6)JztYb#`1r&i@4%2Jef+T)WNQ zt}33J8X#&2b7`B5O;bI)`{@&VZVt+agjWHIQ^(*qI;SamHo=$BJ9!cYqSC?z^nznk z^ND60FXC0z)f>QQh)|ACNWe~wiz@O~j0N~u@wTl64g&g;*vX`_?K=-v|FiQD1Xp(4 zeCgD8^y(NEfCe6!dOvKnZ3JSBW4s{9e_Otii`hpsC#u0j?aA(ks70`=`sf2(t**^B zBh41KMUx<1`SY)CU<5vl?^p{9bQM1von{m#RO=z;+$ zMlD&}BU^lMzc2tC*?E-ZRd{hAkZ&;C*6*=?4;6#OyCmZLAizK_Kq6V zd=>mVfX4QnJGHRdO7IBA4_Va%t4l&Rz}$KJ?p-6S;pl<%(l++ySONP53#k4)QrTFF z0j((L2|e}MUUn>+Np$nus9#!Z9(M4E+5pv8Y37R)903eqn5Q*Un6i;vy2HObQ5N3VA+F ze-6@sg4zmOpqz|9Nzy&wG!|Gqg8AV-Y<=$hc`^@Yk*5pD+DBI|cUT9g(X0S7!Z+Xa z^|c!N`XMATOX4`O3y_XO1HsO%M$1032E*)cD=q}o*G$dyZN^OBNXP~pwG;lnRre?k znnETf|E@GN4y}%vD^LIN9W8z1kCn#zHu)IErKFqya2XgHJ|?Z7PnLQyL5bCRy{^-- zTVNT*?1twpTDeXn@bi#&|0rVm-Nfi3S*q}rZ##Pr-UNHYOX{E7oWefzi*+K<*eT5y zFI@NvmHK;Y@ZUci@*PM;$?n~|w|Dj@H#avCJctB{p#{WcI#;kj(<*P`gm|$1<%}fo z($rWVvAUrR_4~LUn2OwwK;47Il`nr_bc^Efk+a5zF=2Huhq6J6u-PU%tXK!Pl<;;) z338WoJ>ndv$aEfxKgRrBWWmVm_2bA2#gKe+E&5Xv%VtQ5%yOGuuwz*af-{B3fPND- z2*A^o!*lQ@D+tg~Xp(}TO;+sGQ+vn1emDd>j3`Uu>>?oVN0z0~Yhp=toP!sr0x)Ef zd=x~4)u=+wbJ02l9T?lO!tKyZ}3seNNj0B}n2S&rC^Cf{n zK{Pwp&QxLgTULKvVgSsl9)Qe3?loU;t+%xwXHu5=B<%?aA&w`S8R3h=__4eb1k=FQ z$w7EQ>MaGs0r-*V3@*rFA~XE`UB{6it%cv1%c(MW5XoSU{bS=B@!Sf!{pDA!Uft7GNl`f*D&!l?vOWRb$gZ~6R9|A6 z@M22yI4HBG3OcY?%ioFoL;xN$mT=D4|4k{|fb1#@p=$6xu~zomObNUWTkOXI(R?KS za~W=H)jN-M{dII22P*XvkQ|4oiKqz95Xm-CI8n<*D*{9e;1-i5cW`cssg}nhP1Wy& zi}PcAJVX@P9LuaK^zOZTx-ccNqW%;MWMruwNkbD-9b%ynHi*|qMG8`y>&Z1QMnrK5 zxIk#C_axoJPZk0IVBm{S$RxyvCt(+Y@Eo{Te%lsN^=n`Z)bB+08E9sRqi2m>#{al} z`}h>J@DfNMSvMlC7l}zjmWsoM=87f?}e~dlaTde;qG3n zFv2-V52+`AQkwj~r_uw^GE%V{H+H5Z8qYSoyg^aMV*ccE0>i^4kg&sifwlKLKJle_ zg~%QeR^(>F8?RE7rU@#Kc}PFyB!-nV3jDBDL&^uN=f?ul8QM5EQsr1HBeQ~VVB`q! zmW=nyP|W{+V^O1Xn z!A6E1q)9k`c(e7Mv|0QROBw0;&1h9NvMj{-a@}RNN~0HVgQAky;4S#ys`cv!t}!^M zoOLMmXkJ;bxUx=r-Mjv{bVe<#`5f>qBD`@AL=|I{1V4P(nfws%SoQlI8K11-62Z4f zYDhOthnyvbKk*9p*tMmAAd$TbT%3c0!~ezRDdx@J)17`rxDS3}g_tjL4qZ3!`9wp& zF(qf)5;lYvf4aTmq>S}Ry>pK%>YZn>Pm^$^;O?JJG>xrtIJ zuXwDTUc!SRJ%S~H;UcRw6ve*38J$>h+1;HcgMJEYmIzGkNHj|3CUF4iInOTa<{JK< z0?|NWv?xO^=X(c{k`Q3zyga(@Y#YP9S5;J0lvQiBs?AuaO@PeuN~mf`h6*i2KH`^La~zg8tp;wq~BMb?UpB Tbwvt`{4?W-`r)*LXRiH!$?``e delta 13651 zcmZX*2~>@18$bM1$lN&xMMx+#sN+y7O~&kKo|Fb88Z=L&dcvV&DjI1P)!vPoH0%^Y zNTu2d%__A^g$m98-`)HDzxA#4z3ZH{&Z6CZp8LM8-*jDXMm|qgKF_TKiaTVsON-ii zc(}ML$;vt&_i#Dm>~Y5Kgp9M>Y1y}N`)VkPTHtcf)SaRPj+1|Rm})6@6lJz1n`Q|b zB?wDfm74SMFEzX_xb4K4n+Dtdvx} z)$|XceESAXOG``pw$jIf(nnlYm0gbt{I!;{7YPcW)a4XH*R2+(D`wRAe~tCTH>%By z>G=Kb2()i0&i9}EuF>K@F<|!CCtu0y>)E2xaN*0dpNzIKjYtcm0zN%Im1{THTqL~1 z_MKvZ>BD0e28;V;l2THnL`B=vo96KI_uoAoS6eQi_SN@m^{%0?f`S6osge3V5%$YY zzd%N`yrX_kb*%Ui$1u+8hEC45c!sd#&A(EP-OaPgQDw%RXQBe} zC8uQ;Ywm3G%}qXhH_Ax<8{}`P#J^?=y~8mY46fetDzz;qmgWFJHcV<~qo< zdQ&bcymKgX(JfmJK7IPs{=>t2!`Z)o(`OaeEVZ$*VT45OuuF6YdfUe|ab2gaV=hDV_lPZRb?w|3;7O&5?@g#gXzV;7_3Rv`9Ni-Ax$~;vPMN?4r8$vtQuc2MNM=8pM z_0RV0;h~|OecSqCLPGvJ=i^iCb3j^rUsnp`Z$&IhQ`Lpy}h<7V_l)e zGt(1Gd^P0whg8hPH*fAy`~IEvK*1!~AaL?%PF&>Rt>wLarhT{8u+C~{yScf||AUe< z&RBRWbhA#X-^5R0Ro}wFlHqDezp7j7_j>fcTkkjdz2e%nYjqd+_`6&bVruy3&~s)R z-P}B;uOTgEsBFfvJW#+S_f!?9jGn0+&EP1M{Bn7$kiGH_gOM`5Vv8UzA7y^aY=tP_ zyv6SVc=%FFrY_E1l;^?K|6h-)hb$Knv9ybBeSF?zeAqfzB^iZo_u|@$$-gG0N`B$o znb(y3gwp?S1JeJ$4pi#@8DV2ro48+7TiduMPGyPCQTDf1K8oUK<@e2hby@S!q53!B zQr+F%NvPQ^n>TBauBpARA?G*l(c@JSyw1eM<)NmgCb}ruu0Hwuf=ThukNNS%(S05P z3Tu{rJR|S-^W$mKOk#HQn=lFCe@^~~GFCZv{`{d308hg6=lfPHUv7(M5Z>jLzG36W z8ioDyjVjIyP?YY@1}?Q=;X-|U$FR^b#Fpc*?XuzsnbcgcX-Ubt^hG&--g-1jr1Gen?IIQumP@Nqfy}F%kZQ(b5 ziQr(3^wPASj(0~|aZ+^F^~6^Mt(9W5!0|Q}dC% zR~FYwWtDtx{?he2RCxQT8#y^SU4TA7Qa?Jm!NkJ5#b$7I8qe7w9=!hw`N|w=69IsV zLw&MVUvp9JgPSMT?+H~6-C-y!F1TTbUD<<~aCxU@rgrS_p_cI9A3HTCCta^jg^cbV zJ-IhSY3?FC>yz-Cp~umEsRQyQckkYP^ZGR_Gji{W@UXB~_4V}))*Qx#3l|nHS>idm z5g;ax}Xv*hK1BPqEXvlA^nTYhB zq+Mcji+v)LHfc!48KlR=tO1Zx?YxG=n!TXirHk)|Gya$_fy^aQz)muB+(6 z#kr#Eb9q#n(|&b_r2gvuwz^JdN~$8?fWs*1F%k8jzVo7HU1a-IDew2(3T6N4!hLe3 zZ8(thzVnKBPD){?_EX-_@1oLl| zSr`5wkMj5NUxO_idJ;lBpFfFQKAh z0n+WOPYL$*^)*)!Obys5D)$zEvBRzF;MC-Bjii*+u6}tC(3KRLS2r@ls>~9f^PG2Z z|NhI;m2Te+Df{h}IXmcS-z6*mRxg{G8j^EqR|i@20)RbK=}%XeZ63jY7TO!e%;nK9 zFPZ#NDJdiKzQwz3D;P1?Dc5hpOXXLyyOQtkhaFssv%kMF%f!$yY;tN!oZ;l;g!icK z*;*M=-a+58G5vQzfc@<|cbJiuK8)|*zmHE(42`}Jvf8p^hv9vV;P-$H$*IL=g)aAA z-dNT7%Z1Tr-bY!sy*R6E`NYhY8(_1Y|?od1~ z9|^sk|25c9=-Mg9z@r#eMr>WNY8A^k!%SkshCRoRAGiOQlf=u%hlimL-=Ye|R?DB! zg&2vwecUv{?9mxkL95Rpf}6j;caFZ^9M3E%F>iY0XaJD2P~PT$lVa_DyuH2U&V0~qDsbMpV^u(LYGUFZ^8Em6JcG4kYl&(9`wTQ;bl#~dVL!ZCh*)%0 zvfx%T=5G4ZJhNgCR+PMB-{;pnoRjgMVn+@31(x3n3tN6UFfe-EU+S?rj!p9C&Yfe~ zR7NaZy42XiBk#)MRiAyQ0U)1>+*9Pwe4rCX`x|j{j?ei~`+|bWW5YfDAa#Y#BK~#0 zx@t3%xw!TFcve#BfHcpiVG2Srd$M+qkM^hFdbe5>8`jNaP1Or|#g$6Eua4an9ukt+ zpUJK^y3fnpXq;u4@uezySJojZSy_{F=Zc24kKh|k8Z#~APQ6`OMaL)_tlsJ5G4O=4 zSFX$%nP{7_jNRoO(c+}mU$wQ>f5LI@oH@OyNE16d)~7tX&Wj9J8Tso+LpUpwl9Gx$ zuB^`U98Q#op4<(L2U=l0cZuWqd;)DxrV9wT&$=Slb9{*!0XslC@_OQG-FWq5z~ zFRG#agrwOm{_pAdqm;{f=*h6i$jxeMY7J>7L+MS@raAkj#zzfMIN`T%A3wV%MgVP= zZu!{kR~ubcS6?r8^3@-sy|sJSuV4S;+c(|ayDvU@0um!znuI%*)2 z&7M%Cga0MR3Q^eroT+99Y=7G+Un2M1v7luWAOVFeg=291l&6Q{HOsTTcV)@yH(;Js zo2X;YKEJp=T6aj68+U%NmOWG|2q2!3VN=tuHjBKizp#$G+^=c zcz;8-wZ`SkmsvpEYd3GwMrY3!5^)Yqs-ROa+ZB$V_W7D$P~g?)TP`U0rm@i^gp+Bp z3zRO2S3}kP`1Q*KhRdwLDdXRN=e($>kP;Kyhob%B>p(>|yt=*7^IQAw$D{9$+Bc+z ziy0L+g~L`5-Gc*Sq+Gl#z@W^pLGW;IYlzT#1x_tt@$p%2-ySB73~x~ty^~<|Uw{3z zeD&(I%Zr306%|c)d3|*m{qS3MQcIT)rRp%!TN|~@J5yNx43qRJ5K;+%91;=|788>` z&{W{?^x|9?w>y)QllbwTzHAmFQv_Usw{ zFcaYG`SnBX<@(R&c>iZP(nf}c%>CB|qm(>TiNm|EA)qq+b~|}YVl#ToL@AddRDrsR z&&WF=E8lyuyJ+t>RaFD&jj`d8k^N21Wkx4YKIm)A8oJ2f;cRkjdei`Q0$;+%=&`iw zeRprKq=bYX+iT!(P@v=nb^H0zRisN8DXMmh7cQ(rB_!zGKfukcQLP90p{#}oaEKSQ z`LofBo}D*$_FWQCeFx2NWNw~bSEtXQQ?#%C2Iywe)5G!7K?@Vi3fXKEZl*E%)L36X zWO$QmVUCMXx^NEGJj)n)tn`PKBx z#^ryCniN9V6T1cNK;V_ZxGG(#4-E4h$ZwQ{EJ5!Y;bKwT<|;A)N+(*1JxOXox&aYM z@BLuax?r_wW?tTo6@r3&U^{erL#U1_YyI9UF5@c#4!vi4JLfw#84uAV!$HH@Mrp<@ zG~~6QAg_*-g}QTItr*O2w5-KAUA}U~{_BTq9~e9f z6%iSlD(HCpj2O@mLo>5fs1;kh3;-!)c&WgjfBuO^(FYi~8+lNLF+_B9b#-aFc5t~Y zs1Fh!xco!qhnIJ}yS-5Z4%>tqbY4+W@!Iw4?s$=vc#*d`(vQc!MEmryQu3e*9b1Zx zc6kqS23y$qHkE=~>7$Ph!=e2c88Jk^o-~cX6TjGFMaiTAp`|u$dRa8@71wIsnx%nf z7Ql2JD}H}Zv6Nt@G`{%j$jnt3CX(8q7hjK= zo0{JJ`Ss&{0#48ypYk0d-)C9TbnE`7RDx1EG%#MbI$kejrb_E?%3LOTr*j0G%_aik z-8)^JsvX#sG#^x?c1k3rxZQ3Xc~;r98BvdT2&hNWmV?FKmpeC%NPPe;jXdEKIH=?qWQ}rdZf+*| zj+}FA^SLPvWL?YGt;@nwfu4%a+Ee)!$Brc^D=QnC)8^^mcoHd*n#LQPu(oa}_H-c0 z9mEBc_3b@-_bMv!wuf}}^f&-~$-yD#$v-!I3P-5}Qzi*z1*ifpKEPLig8Q>8L`VW- zt|DK3fhu|iShAg&o`QHVfRpta_>_dsARpEUqSxGRaWWxSgVYHe)&KZctNZ9H%09cK zL>bb{mRpx(Uz>OtUNa9`f=Fu=z&n_km>!#nv9Z^nN9f+3d8?>YIO^MO?(X^cIUtu7 z02=1LuGq7G1oo7RZt))dPk~?iY4I}<0nwz8;dGMEXrB?c2-3oVnB|VZbz%$><*fWu zeQn5{J6mX)78es^x^})KtuWO3*ho+B9#{ZWA%#d`)v8ta5`7;ZrSClLAz={_eG@}v zjmSGBW6in;1_nyBAL+!aR#z-qxUj3Uv-6ss-hRj;mvPI+rGa}te*8#(%C-&!L9l>0 zYt{VCNy!8UH?H+=?&#A4XF`(c>*~H+MHjJ)ie6V#XsD{H+P9VzL5kjwb}ls~DPq)e zb;s-zFBYnNd%_2@^qzMR2eEMBLJoL!`=ei@qjD-6Wd7jM=OZ^*9is^05DQJ5vQZ;w z3yr{#7(CC>29q>;W_scL`7Wb(sJ;OFCDXaB)TsDL{i)p7EwN>%DDx;dO-R+$x=(wD zzf{#%tfSM{`TyQ60)C2FDs-uARVL{ZlG(P>GkzJSxd(s+3k3ukm+Wj#ATj}60#b+{ zj&N_D5Dr*?!l(7o_6V}P4jlXDKL1JgL$3idoHsG|=Zs%LU>Rjus(ikGWrKRA?mBuh zdwL{AjX-5FH#yp)a%%z8s?7iN8m_{WkSfjP`+L2rs^#|6)f{9JPA7$!t7-lMjK2LU9K zyLKH7@%a7icm{-`d^D-6>U_0}7cX`V4b?b$^Q<&Da)b#s)iX#wtaoYN5<^fMF>pGn zV-Kds8_S@=j3Air0@;$%)1`4^#E(KvQTBy?^u$jVG!AlBKwy|yQgU)ClE4t{mTe;= z{5f9iSOFUcinnXwXrS8U?ga}Lys57C%Jfo&HShlNfyX1 znK@6AJ`^~|Tl$aLkT#Q(Gcq+bRf1AFBOgt?cF8&S7tfzF4N|npWTDt=AQdsFN$|}! z^#G18;<3Cu!?|7Pc)DxQq!nY2M+K9j(iIMjk92^u@KKt+|4cKum9XkL4 zTiq+FN*P8%m8?~C+lSu##P-gD*9yOTmpl2ZqdX9~?|mXPV3oW^3x6+(xoET(86j{R zj=*3xg;7A~;0TtiRpN~G)VpR|S1cu2#JTJ-ZpBprs^JBI7Wq8PaHhvg4cIms0pNqL zU(X=)m5`KtbJ-v(wc$jEiC(uaBR8K&%eV4n63mu;$0`aQ*kr3Qx=tQcu^!0c1>Y z|53EOXTnG3JGGdtRr1iq7kkvG{3`T0W5~n9OBLGU!{plW?BwiGnsyLj z3l}fuKoy&-hy<7sRj3gx(tz6`rpM;=X{n7H_n!axRS&^D+z7h-L08vFli+KTqBt9) zQ>GEagU043m$tj@YHJyhuB7#0OW`-dVq-JF1S~MzmCesz%>J;xo}O#x^XLRw8ENCWpo)F&0wD^564AqY9y4zP6!Mj9}VkRM|~@vKezqX{#01a zzXT#W_IiBcDe`zdkpA0A!UgRZg?sQ!4w!xvkN2gmShlS5 z*;O9SJcsw`dch(p13^ZddA=JW6pj;HBL^MnMGfO{lEsOQB-*$HD^B6 zUZOhQyAf7vc9;iK2{0j(hf!|5#=xii`mY~zN4Y8S_GoNemHTKRWq$!kM?L_-L6SQJ zk&+gRaG>`j43c|m=;_?)-xX^1pYmusnHWK#YvJ#Bs2(($1ebH}1Hqo@$sV=JmKFhKl1$*XQ8jnoHq9i`x!2>CBtnPaG@*Wue*Ri#Bdm!)N#>g2b* zQ{A%kajjdz|AuzsxJKV@ud#&y?5j(d^F@g z0V?uqO1wV-$jrw-&3dd6Syt29@lWcrtkj6GgOl}WU>44wJ)4q}k{yU+9uHUV5>|A3 zK*AIVKxHHfIvkZy!jPaf+WLH>hK0a;cWllkiRjl#3| zl?j?k-h=5gQ*AR1jaGC9Nh5Kx22fC*BfWY+wjv*8D*qq`Y)T%z-=do^j*!~E-2gNY z>9@&n9pY8#Anv1&6wkPd{Bj> zuQmo3o;9Q7KH>m8>a})qVS5lX>4Fi7i9rNK3=b#_F2@)%E&u7U(<+nnw{j8_q*c6} z+b(kiDC6;{e-~gL!D^z6MRBip*+lGto5Dme3LWN8kiUB41_v2f*0RmoTco58vAu_k zAt_;Y^L)qL)~fgzhH&!i49I`rn`Giv3SHX>>>4In1?kj>5C3lQ9j>Z&&p;QXB2uEs zbwcP84SW3fdWQ=X(-=r5DJ^XXV1SzRY=6Q-a&`bI9Kn0Q-e`r|bZvw6XI_Ge?nq9f zNP@)(fsRF%vctJmndClJ^i)1v%jeIZM>qyLEJ!JtI6-yCZ^(;ShazQ)5zo}HEnI2W zY_v>wXXi#F1bgf6DUn~AOY%c9Hipc#N7C}?;qgm^3_XWhEz~Bz;Y5YS5M?DDitK|u z*VNMLb~-}ITt#%}?j8znL`smD#-)MfMen$@l(cjhLc0zLS61O7L8?97IMZALvjP(N z(Xhh}ceYv2js|{r29JL$7a)gAnV(_w7yi2T_QfFW{BLyk@5@@32J-@5e&jM>2u2;{I^w=mqC$Fh9vlWJ8~ zN(Qh_EyWUwiuA`9%lMdFQh@-Jsm^sX#BI+L5J_)gdsiZwxG~=pLIU}3S0ckz0)xvZ zrbM5`DtYQZ_8IL$Ch40F#q?uz6fV@IHF5?C6L5tkRWVn0Q92}OhFN|c(VwJNQ7-TC z_XP2X4r(-6N{fk$vk;aO6uM3x-U!4&OvK68e-RUb79k4V-SY&yMF#&==H#8O5%C;s z{+Q3l$NY{Kc>q#Og`yyNz1sA50U~<^)FvMSg4ZJ+VW6_S{Y5WE1$qLsPj8Wy=Agag z96#--y0dL>6^*cM-@bpBiYPPjb=>}}GGi=QkQxTbO2S}R6nnVNN=2|I6ATD6jI=RD z(r^cLGm0DQC%Y4AS*4OlA7d3<_kjFr!BG|}|H$xyf)p+xjs54Jf5L9xzPI)0BU=b} z5;DTXjt;Vi`e@8-jKEl&9GgQVUPZ>D5CFU zZpxyRxw@;*e9B|c&&Y*Yj$u6i zWR+7fchlT?^ZMZFrm7sbR_Y%=eh(brd<;QEL*oAHBX}e3XronRU|`^|Yh?~10+Me6 ziY2#gBXKNI6kcCHbUeH2nq5fGUrV`M0lI)qvo)JOoss{{_e7ix*);;bFuW$f<<#5( ziI0tqaS`C??b&k?sX-Qsx5!rCToPF8as@C#^cf%*FO$k-{=xs2$6OnVTCH zN{%e2WZn}dKuh$qQsUg^RdDOGyoOr0qKP4&#(w`MSO~fe$;!&Qp7L^5zJyO=(%08k zRtBd+!a_7J3ysqEA^ScIU=`PxN^QxurYi_mzXsES{@qkDH_Y4|cXkg++O)NUGE2wp zwc?ab(L-(>d4-0}n91tC{?T38-QS-^8W@2zerE%5g0KVU&YmUlBN7WEOtX+pC14H{ zbo3i}&pcS>_b%Qk@^J#>PMJgtI`Oa;vV;@kZ zq|8$RzvEwXYdGnq+mUl<)hG|&0Hi#)z(?~JyPPmzNH~*$GQF~Ntt1@a<$wR(f^??Q z-2LY8*yN-y#OXm|2!I5@Ma&=}rI3P>NRZSi;Z{fwE5FfjBD`=K&F5Ul3Fu%ON=icF zD#dyVUxwBJBWz*|(*+ud3bUZ>t5_TgKW!vgg`B#D8OV_DNW2T3+js$aR z)4%WITMf7o?TCX2TN1VS%noS2M3#DDhoc54WO+y(X15HMG&R_z4E^h=q)oN!H@jA> zC0TjtBp*`JF7o?`^Kc=)4OVh2I$p-D--2o&Qy4Pm z+HvCfe3I>e0ZD@oBB;@9+G>_}5D$g6zIlDqeovxHwc-5?@GnU}`t&HydLAAgvVkLK z{cG-1J(tW#}X(Roi8C`6^89SNz|vFg*CTWTg9 zkQ6>EQu@xdT1>T{aQyIaHyRa15RM?i*N^r?uqMF>_kX!~t&$;>r;wVrJ&*Aoe(GQd z3E{w0ubc{qUL)QMcg|rMY$mvbK{Z(nNvoKnTbYOMhiD+MNSX~(6L1-7h>4(r3~NBx zT^${+?*;bj4x_3mMqR2=8o>dTAM~q@5OKGN&4dGKZZ<~~C=CV8MsSczqyovvbk<5? zfIw(Y3{a&9Z>Y&&)PhQ&kdqKCkje74ZQDo;4Q%(sfiK%2w*~}oON>A^Oh8sl%HB1~ zXLb28BhOnnj)Fsv(5hG$O2b%!ghFW*^K>t*z@8K!v_TE)cyjifgTcx_{zI%820=)o zQ*4WnR1L}mmI>2-S;T%YY~~x#LZXcEsUaGFChPa(_2!!`$n=*<%)>O*>(9WfE2l)dkGt$!!V;!kV$zDo#Eg}$8+Vn0Nkg82_b%^L989io5ni~(?co!mm$k5Vq zuy32*&Dq8vlg+@u3dMMDg*9uIV9p<-d`=(P7wJf{s0W0$AZ`e}^Z>K=VjpG0WlrlU z>K%;OXsyONk~mg31nTTUF0pXYA|x9r0Rc3{n3>6wrV0&TzkW@UGh&oE97AxWHna_) zg)Ga*oQ$J8waZg&$Uz9om@`*xKh+Dv6)`envd$9GQk?&|-PDgtS;R=f3N8tj^pq!p zxT3QsD7hBA9;t>j>==<+gt{V|OHEE}Ei&h(%gTuP%$%JNU$Z}e?i)c~eH250dW4Yc zxtzXoOj8>+1}i7x=lihvR}V8rRPKHa4Fby$+3H#aZrz&1=b~Xzgu+T5%<*4CL*_zr z@`X@>?~oZ2sUn3qVeA8$5+MjA^YC{Tw0Da*qZmtGc`hA$Ncsf=$P^Y*--=DCWC%!b zn{D}(CB(&G!IUIF7+^?l+H`=7WJ#|5=b!gBX+^n=pWga(e=jzG96uRuKlS$A)`I2@ zT#lRw8D>IX3CdbZuUof{QP*7LP6mQRq#^FlTtlzIS0TRxXozpxBq=VwA9zH9AEL+* zL8BWm5i|h(lW5oj3wvyzlZJd`e|4=lDjOq^Eszxj&2HDB_NPa)Xf6aM@#@=E#a<(v zBWq?D@%GP~rz_Rf-94MV4Hor?w7A!(R98sRl9G8)qj+;kSUB9cS&>^3QY&+Xt*Og2 z)Tx9ZcZD@RU$ZM_7yn3l1{NhN8H_rjvj6()b2CBWdSt;Yl#NNy;)E-eH8qX+?I|y= z3prqBfk^jOd53~^%RRsgNxsoYnX3gGf-cYJNuYyrn+A$8E}hj|+`9d!;>SVoMyKP> z9e^*Yp8JGE@AeIn*zWTD-(ACLk*ydcH8Er?8d$y_Gr-fnG1I=2n1a_oxxfcP?0Vy> zSET^G(5FzDL+Tg1w>nX3bSW84V51C*DY?!op7#B5m56A>0)>SBUXbpArXh<*b*x=O zv07)&q=OIBSANlmiF>UT5L4{7(lo-z@6dZA>G^EAf3hQNd&Q??p_IXR%{Gn)x zI>AT)ffhmB@WHd5kZDCOIH3#J{vwk=t$0-n5_Mw@SkjOB(SJW1b*~E*hsA+Pu#9dX z*)k@f4#FWKim=(%({h@mr6k*jl^}0Pqcy#&hx2*b_l$h~Xt1;O?0qsT#N^xlNdOO) z&2%v-!lF~qI<-t}xESFYiqLOeeYo8t{J4jQSV4eEdW|M&v^<}IX-l$*iq9*=VQ`ix zswoZzM1$2gfhQD&IVxEK<1+mB6KRKnBWrL*=w+wMJxsErllA`@3w z0%U8Xseme0+%a6K1Ga3iGGs5bVnvdG>hIb0S-7rQV1YuhM+WBeYcS*9-oU>1N1Z)K z?-VX3Any4AqqI@DgZ`;sZ6K&3xA^Dz-b>(E+q4!H#f^1K0)f0N+36tm;V zg0)cW+?piRd{&L>#Nq2gvc{O=5hNATnMcvgB}7Fggw&@Ti7i!qj(hWNb`6APGFM(j zCCIm9Grc_>feg?M4{;yu$!X}sQuG4Sgj)j3JCxtJ$`F_Q$T1pAIE%ICQOqjxy7y5v zXXKA2vCZIQ(5cwyX6PgFh7U=?^p|aam{Q+bPRpTE`H3dSA!1Zyf@BD@O;R}9FLq)7 zjikj!CKlb$+*JO8f{wWXO34U)OG`^b#4wx7Kl4%l_U+qb(`Cbk4cM1SZf`%9Z(F;E zOm#?vLd@Bjk999X`6zon^2`ay=Ls%AqSnin(UXa&Qwi*r)v*k62@@hB!Y%~qI;2-F z-=E|EYH$(D=Srt=YB|z}aP;9KvRI&k?I%ksL{JifMCL+-A2!m2v=&iJM{)$#kZy?q z0Td~JtUTh=aI@i=Q@`^OUonw~VE=Ls#W5zfs_zaCfp|nHI}+ot+2@m9@JmTb5)D$a zOz=Ao#nk@a4C z{>0<#4GBed7eT@oKrR`d^u43+T_Sr#7zf4=D`ERRbrCAh^VPK#gshTK;I`O`9RCvS zjHLrH1_@TkfVkA}yeGH<8knp#BU`$jGLPbHnti-mz)Kcb+-IyilB9@GT_kHH3yBbo zBs~GRz=c|qjnUsUHXI>Mi+#?Z8lKHQC4qtSLyt}H#boA$$@P#gY_2+Pq#jm~gg2-+ z?#tEA9vykOo&S*V=>`@S6_Eiwu_nV8-J`2s=C7kD>)D>)hgK!CER69_`+gWW-MV<2 ztgH^k1{f>Rckws*pnD`eHSSg#CC?jh~X%%xRj{X`R)FW9s!4Vo`4@HEYbj6TP zi`aD48R)^I(OUZ@L4t3>vX#K^!583cwG z9YM%oE3c7ye{w39`;>njhs@8L-D|i|7YZ~c3(HM*rK@neKC9%3g&Pzx6`@iJfArcGEYEtTs`-m*_8n-2r zH4^tCDYhc{1m+bU7-_JPo2zfd4^c4?;Bq+i7NK)Ju^9ijjVj Date: Fri, 11 Feb 2022 13:14:19 +0100 Subject: [PATCH 5/7] make grid_values to be a Stats class, make the corresponding changes everywhere --- validphys2/src/validphys/arclength.py | 4 +- .../validphys/closuretest/multiclosure_pdf.py | 2 +- validphys2/src/validphys/correlations.py | 32 ++++----- validphys2/src/validphys/dataplots.py | 2 +- validphys2/src/validphys/deltachi2.py | 11 +-- validphys2/src/validphys/eff_exponents.py | 12 ++-- validphys2/src/validphys/mc2hessian.py | 2 +- validphys2/src/validphys/pdfgrids.py | 69 ++++++++++--------- validphys2/src/validphys/pdfplots.py | 15 ++-- validphys2/src/validphys/replica_selector.py | 4 +- 10 files changed, 76 insertions(+), 77 deletions(-) diff --git a/validphys2/src/validphys/arclength.py b/validphys2/src/validphys/arclength.py index 5ae170fc25..3484c5742e 100644 --- a/validphys2/src/validphys/arclength.py +++ b/validphys2/src/validphys/arclength.py @@ -50,7 +50,7 @@ def arc_lengths( eps = (b - a) / npoints ixgrid = xgrid(a, b, "linear", npoints) # PDFs evaluated on grid - xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values * ixgrid[1] + xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values.data * ixgrid[1] fdiff = np.diff(xfgrid) / eps # Compute forward differences res += integrate.simps(1 + np.square(fdiff), ixgrid[1][1:]) stats = pdf.stats_class(res) @@ -127,6 +127,6 @@ def integrability_number( checked = check_basis(basis, flavours) basis, flavours = checked["basis"], checked["flavours"] ixgrid = xgrid(1e-9, 1e-6, "log", 3) - xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values + xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values.data res = np.sum(np.abs(xfgrid), axis=2) return res.squeeze() diff --git a/validphys2/src/validphys/closuretest/multiclosure_pdf.py b/validphys2/src/validphys/closuretest/multiclosure_pdf.py index 7881c1e557..a8d0f5f3f8 100644 --- a/validphys2/src/validphys/closuretest/multiclosure_pdf.py +++ b/validphys2/src/validphys/closuretest/multiclosure_pdf.py @@ -90,7 +90,7 @@ def xi_grid_values(xi_pdfgrids): glu_sin_grid, nonsin_grid = xi_pdfgrids # grid values have shape: replica, flavour, x # concatenate along flavour - return np.concatenate((glu_sin_grid.grid_values, nonsin_grid.grid_values), axis=1) + return np.concatenate((glu_sin_grid.grid_values.data, nonsin_grid.grid_values.data), axis=1) def underlying_xi_grid_values( diff --git a/validphys2/src/validphys/correlations.py b/validphys2/src/validphys/correlations.py index 7a254467e5..2e6132c3cf 100644 --- a/validphys2/src/validphys/correlations.py +++ b/validphys2/src/validphys/correlations.py @@ -14,12 +14,12 @@ #This would be a good candidate to be optimized to calculate everything in one #pass over x, -def _basic_obs_pdf_correlation(pdfarr, obsarr): +def _basic_obs_pdf_correlation(pdf_stats, obs_stats): """Calculate the correlation between pdfs and observables. - The expected format is: + The expected format is are Stats instances of: - obsarr: (nbins x nreplicas), as returned from thresult. - pdfarr: (nreplicas x nf x nx), as returned from xplotting_grid.grid_values + obs_stats: (nbins x nreplicas), as returned from thresult. + pdf_stats: (nreplicas x nf x nx), as returned from xplotting_grid.grid_values The returned array contains the PDF correlation between the value of the obsevable and the PDF at the corresponding point in (fl,x) @@ -27,11 +27,8 @@ def _basic_obs_pdf_correlation(pdfarr, obsarr): (nbins x nf x nx), compatible with grid_values, upon changing replicas->bins. """ - - #Remove mean - #TODO: This should be done at the Result level - x = pdfarr - np.mean(pdfarr, axis=0) - y = obsarr.T - np.mean(obsarr, axis=1) + x = pdf_stats.data - pdf_stats.central_value() + y = obs_stats.data - obs_stats.central_value() #We want to compute: #sum(x*y)/(norm(x)*norm(y)) @@ -46,17 +43,18 @@ def _basic_obs_pdf_correlation(pdfarr, obsarr): return num/den -def _basic_obs_obs_correlation(obsarr1, obsarr2): +def _basic_obs_obs_correlation(obs1, obs2): """Calculate the correlation between two observables. The expected format is - obsarr1: (nbins1, nreplicas) - obsarr2: (nbins2, nreplicas) + Stats instances of: + + obs1: (nreplicas, nbins1) + obs2: (nreplicas, nbins2) The result is (nbins1 , nbins2), a mareix containing the correlation coefficients between the two sets. """ - #TODO: Do this at Result level taking into account error type - x = (obsarr1.T - np.mean(obsarr1, axis=1)).T - y = (obsarr2.T - np.mean(obsarr2, axis=1)) + x = (obs1.data - obs1.central_value()).T + y = (obs2.data - obs2.central_value()) return x@y/np.outer(la.norm(x,axis=1),la.norm(y,axis=0)) @@ -66,7 +64,7 @@ def obs_pdf_correlations(pdf, results, xplotting_grid): """Return the correlations between each point in a dataset and the PDF values on a grid of (x,f) points in a format similar to `xplotting_grid`.""" _ , th = results - corrs = _basic_obs_pdf_correlation(xplotting_grid.grid_values, th.rawdata) + corrs = pdf.stats_class(_basic_obs_pdf_correlation(xplotting_grid.grid_values, th.stats)) return xplotting_grid.copy_grid(grid_values=corrs) @@ -77,4 +75,4 @@ def obs_pdf_correlations(pdf, results, xplotting_grid): def obs_obs_correlations(pdf, corrpair_results): """Return the theoretical correlation matrix between a pair of observables.""" (_,th1), (_,th2) = corrpair_results - return _basic_obs_obs_correlation(th1.rawdata, th2.rawdata) + return _basic_obs_obs_correlation(th1.stats, th2.stats) diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py index bbb1279f50..b9bfd1def3 100644 --- a/validphys2/src/validphys/dataplots.py +++ b/validphys2/src/validphys/dataplots.py @@ -871,7 +871,7 @@ def plot_smpdf(pdf, dataset, obs_pdf_correlations, mark_threshold:float=0.9): basis = obs_pdf_correlations.basis - fullgrid = obs_pdf_correlations.grid_values + fullgrid = obs_pdf_correlations.grid_values.data fls = obs_pdf_correlations.flavours x = obs_pdf_correlations.xgrid diff --git a/validphys2/src/validphys/deltachi2.py b/validphys2/src/validphys/deltachi2.py index 503a514df1..bba69a87e5 100644 --- a/validphys2/src/validphys/deltachi2.py +++ b/validphys2/src/validphys/deltachi2.py @@ -134,11 +134,6 @@ def plot_delta_chi2_hessian_distribution(delta_chi2_hessian, pdf, total_chi2_dat return fig -XPlottingGrid = namedtuple( - "XPlottingGrid", ("Q", "basis", "flavours", "xgrid", "grid_values", "scale") -) - - def pos_neg_xplotting_grids(delta_chi2_hessian, xplotting_grid): """ Generates xplotting_grids correspodning to positive and negative delta chi2s. @@ -149,8 +144,8 @@ def pos_neg_xplotting_grids(delta_chi2_hessian, xplotting_grid): pos_mask = np.append(True, positive_eigenvalue_mask) neg_mask = np.append(True, ~positive_eigenvalue_mask) - pos_grid = xplotting_grid.grid_values[pos_mask] - neg_grid = xplotting_grid.grid_values[neg_mask] + pos_grid = xplotting_grid.grid_values.data[pos_mask] + neg_grid = xplotting_grid.grid_values.data[neg_mask] pos_xplotting_grid = xplotting_grid.copy_grid(grid_values=pos_grid) neg_xplotting_grid = xplotting_grid.copy_grid(grid_values=neg_grid) @@ -220,7 +215,7 @@ def draw(self, pdf, grid, flstate): # Basically stats is an object which says what is the type class of the replicas: # MCStats, HessianStats, SymmHessianStats. In this way validphys use the right methods # to compute statistical values. - stats = pdf.stats_class(grid.grid_values[:, flindex, :]) + stats = pdf.stats_class(grid.grid_values.data[:, flindex, :]) # Ignore spurious normalization warnings with warnings.catch_warnings(): diff --git a/validphys2/src/validphys/eff_exponents.py b/validphys2/src/validphys/eff_exponents.py index fcb40aa466..dce6a49d09 100644 --- a/validphys2/src/validphys/eff_exponents.py +++ b/validphys2/src/validphys/eff_exponents.py @@ -66,14 +66,14 @@ def alpha_eff(pdf: PDF, *, pdfGrid = pdfgrids.xplotting_grid( pdf, Q, xgrid=xGrid, basis=basis, flavours=flavours) - pdfGrid_values = pdfGrid.grid_values + pdfGrid_values = pdfGrid.grid_values.data # NOTE: without this I get "setting an array element with a sequence" xGrid = pdfGrid.xgrid with warnings.catch_warnings(): warnings.simplefilter('ignore', RuntimeWarning) alphaGrid_values = -np.log(abs(pdfGrid_values/xGrid))/np.log(xGrid) alphaGrid_values[alphaGrid_values == - np.inf] = np.nan # when PDF_i =0 - alphaGrid = pdfGrid.copy_grid(grid_values=alphaGrid_values) + alphaGrid = pdfGrid.copy_grid(grid_values=pdf.stats_class(alphaGrid_values)) return alphaGrid @check_positive('Q') @@ -110,14 +110,14 @@ def beta_eff(pdf, *, pdfGrid = pdfgrids.xplotting_grid( pdf, Q, xgrid=xGrid, basis=basis, flavours=flavours) - pdfGrid_values = pdfGrid.grid_values + pdfGrid_values = pdfGrid.grid_values.data # NOTE: without this I get "setting an array element with a sequence" xGrid = pdfGrid.xgrid with warnings.catch_warnings(): warnings.simplefilter('ignore', RuntimeWarning) betaGrid_values = np.log(abs(pdfGrid_values/xGrid))/np.log(1-xGrid) betaGrid_values[betaGrid_values == -np.inf] = np.nan # when PDF_i =0 - betaGrid = pdfGrid.copy_grid(grid_values=betaGrid_values) + betaGrid = pdfGrid.copy_grid(grid_values=pdf.stats_class(betaGrid_values)) return betaGrid # .grid_values @@ -348,8 +348,8 @@ def next_effective_exponents_table( eff_exp_data = [] - alphastats = alpha_effs.stats_gv - betastats = beta_effs.stats_gv + alphastats = alpha_effs.grid_values + betastats = beta_effs.grid_values with warnings.catch_warnings(): warnings.simplefilter('ignore', RuntimeWarning) diff --git a/validphys2/src/validphys/mc2hessian.py b/validphys2/src/validphys/mc2hessian.py index b8d1aeef02..2f9ee758fc 100644 --- a/validphys2/src/validphys/mc2hessian.py +++ b/validphys2/src/validphys/mc2hessian.py @@ -111,7 +111,7 @@ def _create_mc2hessian(pdf, Q, xgrid, Neig, output_path, name=None): def _get_X(pdf, Q, xgrid, reshape=False): pdf_grid = xplotting_grid(pdf, Q, xgrid=xgrid) - pdf_grid_values = pdf_grid.stats_gv + pdf_grid_values = pdf_grid.grid_values replicas = pdf_grid_values.error_members() mean = pdf_grid_values.central_value() Xt = replicas - mean diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py index 7c9f677641..18b4ea2158 100644 --- a/validphys2/src/validphys/pdfgrids.py +++ b/validphys2/src/validphys/pdfgrids.py @@ -5,6 +5,7 @@ from collections import namedtuple import dataclasses import numbers +import logging import numpy as np import scipy.integrate as integrate @@ -17,6 +18,8 @@ from validphys.pdfbases import (Basis, check_basis) from validphys.checks import check_pdf_normalize_to, check_xlimits +log = logging.getLogger(__name__) + @make_argcheck def _check_scale(scale): scales = ('linear', 'log') @@ -42,22 +45,22 @@ def xgrid(xmin:numbers.Real=1e-5, xmax:numbers.Real=1, class XPlottingGrid: """DataClass holding the value of the PDF at the specified values of x, Q and flavour. - It also exposes a `stats_gv` attribute with a `Stats` instance - of the raw `grid_values` object in order to compute statistical - estimators in a sensible manner. + The `grid_values` attribute corresponds to a `Stats` instance + in order to compute statistical estimators in a sensible manner. """ Q: float basis: (str, Basis) flavours: (list, tuple, type(None)) xgrid: np.ndarray - grid_values: np.ndarray + grid_values: Stats scale: str - stats_gv: Stats # TODO: transitional variable during development, it should substitute grid_values by the end of this PR!!!!!!!! def copy_grid(self, grid_values=None): """Create a copy of the grid with potentially a different set of values""" - new_stats_gv = self.stats_gv.__class__(grid_values) - return dataclasses.replace(self, grid_values=grid_values, stats_gv=new_stats_gv) + if not isinstance(grid_values, Stats): + log.warning("XPlottingGrid being called with a numpy grid, should be using Stats instead!") + grid_values = self.grid_values.__class__(grid_values) + return dataclasses.replace(self, grid_values=grid_values) @make_argcheck(check_basis) @@ -89,11 +92,9 @@ def xplotting_grid(pdf:PDF, Q:(float,int), xgrid=None, basis:(str, Basis)='flavo raise TypeError(f"Invalid xgrid {xgrid!r}") gv = basis.grid_values(pdf, flavours, xgrid, Q) #Eliminante Q axis - #TODO: wrap this in pdf.stats_class? - gv = gv.reshape(gv.shape[:-1]) - gv_stats = pdf.stats_class(gv) + gv = pdf.stats_class(gv.reshape(gv.shape[:-1])) - res = XPlottingGrid(Q, basis, flavours, xgrid, gv, scale, gv_stats) + res = XPlottingGrid(Q, basis, flavours, xgrid, gv, scale) return res xplotting_grids = collect(xplotting_grid, ('pdfs',)) @@ -252,25 +253,28 @@ def distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(None))=None set is computed. At least one grid will be identical to zero. """ - gr2 = xplotting_grids[normalize_to] - cv2 = pdfs[normalize_to].stats_class(gr2.grid_values).central_value() - sg2 = pdfs[normalize_to].stats_class(gr2.grid_values).std_error() - N2 = gr2.grid_values.shape[0] + gr2_stats = xplotting_grids[normalize_to].grid_values + cv2 = gr2_stats.central_value() + sg2 = gr2_stats.std_error() + N2 = gr2_stats.get_members() - newgrids = list() + newgrids = [] for grid, pdf in zip(xplotting_grids, pdfs): if pdf == pdfs[normalize_to]: - newgrid = grid.copy_grid(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2]))) + # Zero the PDF we are normalizing against + pdf_zero = pdf.stats_class(np.zeros_like(gr2_stats.data[0])) + newgrid = grid.copy_grid(grid_values=pdf_zero) newgrids.append(newgrid) continue - cv1 = pdf.stats_class(grid.grid_values).central_value() - sg1 = pdf.stats_class(grid.grid_values).std_error() - N1 = grid.grid_values.shape[0] + g_stats = grid.grid_values + cv1 = g_stats.central_value() + sg1 = g_stats.std_error() + N1 = g_stats.get_members() # the distance definition - distance = np.sqrt((cv1-cv2)**2/(sg1**2/N1+sg2**2/N2)) + distance = pdf.stats_class(np.sqrt((cv1-cv2)**2/(sg1**2/N1+sg2**2/N2))) newgrid = grid.copy_grid(grid_values=distance) newgrids.append(newgrid) @@ -290,27 +294,30 @@ def variance_distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(No set is computed. At least one grid will be identical to zero. """ - gr2 = xplotting_grids[normalize_to] - sg2 = pdfs[normalize_to].stats_class(gr2.grid_values).std_error() - mo2 = pdfs[normalize_to].stats_class(gr2.grid_values).moment(4) - N2 = gr2.grid_values.shape[0] + gr2_stats = xplotting_grids[normalize_to].grid_values + sg2 = gr2_stats.std_error() + mo2 = gr2_stats.moment(4) + N2 = gr2_stats.get_members() s2 = (mo2-(N2-3)/(N2-1)*sg2**4)/N2 - newgrids = list() + newgrids = [] for grid, pdf in zip(xplotting_grids, pdfs): if pdf == pdfs[normalize_to]: - newgrid = grid.copy_grid(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2]))) + # Zero the PDF we are normalizing against + pdf_zero = pdf.stats_class(np.zeros_like(gr2_stats.data[0])) + newgrid = grid.copy_grid(grid_values=pdf_zero) newgrids.append(newgrid) continue - sg1 = pdf.stats_class(grid.grid_values).std_error() - mo1 = pdf.stats_class(grid.grid_values).moment(4) - N1 = grid.grid_values.shape[0] + g_stats = grid.grid_values + sg1 = g_stats.std_error() + mo1 = g_stats.moment(4) + N1 = g_stats.get_members() s1 = (mo1-(N1-3)/(N1-1)*sg1**4)/N1 # the distance definition - variance_distance = np.sqrt((sg1**2-sg2**2)**2/(s1+s2)) + variance_distance = pdf.stats_class(np.sqrt((sg1**2-sg2**2)**2/(s1+s2))) newgrid = grid.copy_grid(grid_values=variance_distance) newgrids.append(newgrid) diff --git a/validphys2/src/validphys/pdfplots.py b/validphys2/src/validphys/pdfplots.py index fc5571b812..3b668dba45 100644 --- a/validphys2/src/validphys/pdfplots.py +++ b/validphys2/src/validphys/pdfplots.py @@ -61,8 +61,7 @@ def normalize(self): if normalize_to is not None: normalize_pdf = self.normalize_pdf normalize_grid = self._xplotting_grids[normalize_to] - normvals = normalize_pdf.stats_class( - normalize_grid.grid_values).central_value() + normvals = normalize_grid.grid_values.central_value() #Handle division by zero more quietly def fp_error(tp, flag): @@ -75,8 +74,8 @@ def fp_error(tp, flag): newgrids = [] with np.errstate(all='call'): np.seterrcall(fp_error) - for grid in self._xplotting_grids: - newvalues = grid.grid_values/normvals + for pdf, grid in zip(self.pdfs, self._xplotting_grids): + newvalues = pdf.stats_class(grid.grid_values.data/normvals) newgrids.append(grid.copy_grid(grid_values=newvalues)) return newgrids @@ -187,7 +186,7 @@ def draw(self, pdf, grid, flstate): ax = flstate.ax next_prop = next(ax._get_lines.prop_cycler) color = next_prop['color'] - gv = grid.grid_values[:,flstate.flindex,:] + gv = grid.grid_values.data[:,flstate.flindex,:] ax.plot(grid.xgrid, gv.T, alpha=0.2, linewidth=0.5, color=color, zorder=1) stats = pdf.stats_class(gv) @@ -224,7 +223,7 @@ def get_ylabel(self, parton_name): def draw(self, pdf, grid, flstate): ax = flstate.ax flindex = flstate.flindex - gv = grid.grid_values[:,flindex,:] + gv = grid.grid_values.data[:,flindex,:] stats = pdf.stats_class(gv) res = stats.std_error() @@ -331,7 +330,7 @@ def draw(self, pdf, grid, flstate): next_prop = next(pcycler) color = next_prop['color'] - gv = grid.grid_values[flindex,:] + gv = grid.grid_values.data[flindex,:] ax.plot(grid.xgrid, gv, color=color, label='$%s$' % flstate.parton_name) return gv @@ -411,7 +410,7 @@ def draw(self, pdf, grid, flstate): hatchit = flstate.hatchit labels = flstate.labels handles = flstate.handles - stats = pdf.stats_class(grid.grid_values[:,flindex,:]) + stats = pdf.stats_class(grid.grid_values.data[:,flindex,:]) pcycler = ax._get_lines.prop_cycler #This is ugly but can't think of anything better diff --git a/validphys2/src/validphys/replica_selector.py b/validphys2/src/validphys/replica_selector.py index 200c2df36b..08a82bdf1c 100644 --- a/validphys2/src/validphys/replica_selector.py +++ b/validphys2/src/validphys/replica_selector.py @@ -169,7 +169,7 @@ def not_outlier_mask( """Return a boolean mask with the replicas that are never outside of the given percentile in the constrained region""" lim = unconstrained_region_index - gv = xplotting_grid.grid_values[:, :, lim:] + gv = xplotting_grid.grid_values.data[:, :, lim:] delta = nsigma * np.std(gv, axis=0) mean = np.mean(gv, axis=0) return ((gv >= mean - delta) & (gv <= mean+delta)).all(axis=2).all(axis=1) @@ -309,7 +309,7 @@ def draw(self, pdf, grid, flstate): #PDFs, and instead would do something different next_prop = next(ax._get_lines.prop_cycler) color = next_prop['color'] - gv = grid.grid_values[filt, flstate.flindex, :] + gv = grid.grid_values.data[filt, flstate.flindex, :] ax.plot( grid.xgrid, gv.T, From 64d09b7240ea71355b38ba0bd8a5b7d24a17e7a9 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 11 Feb 2022 18:35:01 +0100 Subject: [PATCH 6/7] same changes for n3fit --- n3fit/src/n3fit/tests/test_vpinterface.py | 10 +++++----- n3fit/src/n3fit/vpinterface.py | 2 +- validphys2/src/validphys/pdfgrids.py | 8 ++++---- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/n3fit/src/n3fit/tests/test_vpinterface.py b/n3fit/src/n3fit/tests/test_vpinterface.py index 5a90d6cbf9..26e77ce2cd 100644 --- a/n3fit/src/n3fit/tests/test_vpinterface.py +++ b/n3fit/src/n3fit/tests/test_vpinterface.py @@ -56,7 +56,7 @@ def test_N3PDF(members, layers): assert n3pdf.compute_arclength().shape == (5,) # Try to get a plotting grid res = xplotting_grid(n3pdf, 1.6, xx) - assert res.grid_values.shape == (members, 8, xsize) + assert res.grid_values.data.shape == (members, 8, xsize) def test_vpinterface(): @@ -68,7 +68,7 @@ def test_vpinterface(): res_2 = xplotting_grid(fit_2, 1.6, xgrid) distances = distance_grids([fit_1, fit_2], [res_1, res_2], 0) assert len(distances) == 2 - assert distances[0].grid_values.shape == (8, 40) - assert distances[1].grid_values.shape == (8, 40) - np.testing.assert_allclose(distances[0].grid_values, 0.0) - assert not np.allclose(distances[1].grid_values, 0.0) + assert distances[0].grid_values.data.shape == (8, 40) + assert distances[1].grid_values.data.shape == (8, 40) + np.testing.assert_allclose(distances[0].grid_values.data, 0.0) + assert not np.allclose(distances[1].grid_values.data, 0.0) diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py index 28af0a0469..98ccd92bed 100644 --- a/n3fit/src/n3fit/vpinterface.py +++ b/n3fit/src/n3fit/vpinterface.py @@ -13,7 +13,7 @@ >>> pdf_model = pdfNN_layer_generator(nodes=[8], activations=['linear'], seed=0, flav_info=fake_fl) >>> n3pdf = N3PDF(pdf_model) >>> res = xplotting_grid(n3pdf, 1.6, fake_x) - >>> res.grid_values.shape + >>> res.grid_values.error_members().shape (1, 8, 3) diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py index 18b4ea2158..ae1eabd976 100644 --- a/validphys2/src/validphys/pdfgrids.py +++ b/validphys2/src/validphys/pdfgrids.py @@ -256,7 +256,7 @@ def distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(None))=None gr2_stats = xplotting_grids[normalize_to].grid_values cv2 = gr2_stats.central_value() sg2 = gr2_stats.std_error() - N2 = gr2_stats.get_members() + N2 = pdfs[normalize_to].get_members() newgrids = [] for grid, pdf in zip(xplotting_grids, pdfs): @@ -271,7 +271,7 @@ def distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(None))=None g_stats = grid.grid_values cv1 = g_stats.central_value() sg1 = g_stats.std_error() - N1 = g_stats.get_members() + N1 = pdf.get_members() # the distance definition distance = pdf.stats_class(np.sqrt((cv1-cv2)**2/(sg1**2/N1+sg2**2/N2))) @@ -297,7 +297,7 @@ def variance_distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(No gr2_stats = xplotting_grids[normalize_to].grid_values sg2 = gr2_stats.std_error() mo2 = gr2_stats.moment(4) - N2 = gr2_stats.get_members() + N2 = pdfs[normalize_to].get_members() s2 = (mo2-(N2-3)/(N2-1)*sg2**4)/N2 newgrids = [] @@ -313,7 +313,7 @@ def variance_distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(No g_stats = grid.grid_values sg1 = g_stats.std_error() mo1 = g_stats.moment(4) - N1 = g_stats.get_members() + N1 = pdf.get_members() s1 = (mo1-(N1-3)/(N1-1)*sg1**4)/N1 # the distance definition From 4a7d7ef4b9dcc4973abbd53f2421d560de3fcdbf Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 14 Feb 2022 12:07:05 +0100 Subject: [PATCH 7/7] removed NNPDFresult --- validphys2/src/validphys/results.py | 42 +++++++++++------------------ 1 file changed, 16 insertions(+), 26 deletions(-) diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index 052eb1f6bf..87f95450cd 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -26,7 +26,7 @@ check_two_dataspecs, ) -from validphys.core import DataSetSpec, PDF, DataGroupSpec +from validphys.core import DataSetSpec, PDF, DataGroupSpec, Stats from validphys.calcutils import ( all_chi2, central_chi2, @@ -36,7 +36,6 @@ ) from validphys.convolution import ( predictions, - central_predictions, PredictionsRequireCutsError, ) @@ -48,33 +47,13 @@ class Result: pass -# TODO: Eventually,only one of (NNPDFDataResult, StatsResult) should survive -class NNPDFDataResult(Result): - """A result fills its values from a pandas dataframe - For legacy (libNNPDF) compatibility, falls back to libNNPDF attributes""" - - def __init__(self, dataobj=None, central_value=None): - # This class is used by both validphys and libNNPDF objects - # when central_value is not explictly passed, fallback to - # libNNPDF object .get_cv() - if central_value is None: - central_value = dataobj.get_cv() - self._central_value = np.array(central_value).reshape(-1) - - @property - def central_value(self): - return self._central_value - - def __len__(self): - return len(self.central_value) - - class StatsResult(Result): def __init__(self, stats): self.stats = stats @property def rawdata(self): + """Returns the raw data with shape (Npoints, Npdf)""" return self.stats.data.T @property @@ -85,17 +64,28 @@ def central_value(self): def std_error(self): return self.stats.std_error() + def __len__(self): + """Returns the number of data points in the result""" + return self.rawdata.shape[0] + -class DataResult(NNPDFDataResult): - def __init__(self, dataobj, covmat, sqrtcovmat, central_value=None): - super().__init__(dataobj, central_value=central_value) +class DataResult(StatsResult): + """Holds the relevant information from a given dataset""" + def __init__(self, dataobj, covmat, sqrtcovmat): + self._central_value = dataobj.get_cv() + stats = Stats(self._central_value) self._covmat = covmat self._sqrtcovmat = sqrtcovmat + super().__init__(stats) @property def label(self): return "Data" + @property + def central_value(self): + return self._central_value + @property def std_error(self): return np.sqrt(np.diag(self.covmat))