Source code for validphys.coredata

Data containers backed by Python managed memory (Numpy arrays and Pandas

import dataclasses
import logging
from typing import Optional

import numpy as np
import pandas as pd

from reportengine.compat import yaml
from validphys.utils import generate_path_filtered_data

KIN_NAMES = ["kin1", "kin2", "kin3"]
log = logging.getLogger(__name__)

[docs]@dataclasses.dataclass(eq=False) class FKTableData: """ Data contained in an FKTable Parameters ---------- hadronic : bool Whether a hadronic (two PDFs) or a DIS (one PDF) convolution is needed. Q0 : float The scale at which the PDFs should be evaluated (in GeV). ndata : int The number of data points in the grid. xgrid : array, shape (nx) The points in x at which the PDFs should be evaluated. sigma : pd.DataFrame For hadronic data, the columns are the indexes in the ``NfxNf`` list of possible flavour combinations of two PDFs. The MultiIndex contains three keys, the data index, an index into ``xgrid`` for the first PDF and an idex into ``xgrid`` for the second PDF, indicating if the points in ``x`` where the PDF should be evaluated. For DIS data, the columns are indexes in the ``Nf`` list of flavours. The MultiIndex contains two keys, the data index and an index into ``xgrid`` indicating the points in ``x`` where the PDF should be evaluated. metadata : dict Other information contained in the FKTable. protected: bool When a fktable is protected cuts will not be applied. The most common use-case is when a total cross section is used as a normalization table for a differential cross section, in legacy code (<= NNPDF4.0) both fktables would be cut using the differential index. """ hadronic: bool Q0: float ndata: int xgrid: np.ndarray sigma: pd.DataFrame metadata: dict = dataclasses.field(default_factory=dict, repr=False) protected: bool = False
[docs] def with_cfactor(self, cfactor): """Returns a copy of the FKTableData object with cfactors applied to the fktable""" if all(c == 1.0 for c in cfactor): return self if len(cfactor) != self.ndata: if self.protected: cfactor = cfactor[0] else: name = self.metadata.get("target_dataset") raise ValueError( f"The length of cfactor for {name} differs from the number of datapoints in the grid" ) new_sigma = self.sigma.multiply(pd.Series(cfactor), axis=0, level=0) return dataclasses.replace(self, sigma=new_sigma)
[docs] def with_cuts(self, cuts): """Return a copy of the FKTabe with the cuts applied. The data index of the sigma operator (the outermost level), contains the data point that have been kept. The ndata property is updated to reflect the new number of datapoints. If cuts is None, return the object unmodified. Parameters ---------- cuts : array_like or validphys.core.Cuts or None. The cuts to be applied. Returns ------- res : FKTableData A copy of the FKtable with the cuts applies. Notes ----- The original number of points can be accessed with ``table.metadata['GridInfo'].ndata``. Examples -------- >>> from validphys.fkparser import load_fktable ... from validphys.loader import Loader ... l = Loader() ... ds = l.check_dataset('ATLASTTBARTOT', theoryid=53, cfac=('QCD',)) ... table = load_fktable(ds.fkspecs[0]) ... newtable = table.with_cuts([0,1]) >>> assert set(newtable.sigma.index.get_level_values(0)) == {0,1} >>> assert newtable.ndata == 2 >>> assert newtable.metadata['GridInfo'].ndata == 3 """ if hasattr(cuts, "load"): cuts = cuts.load() if cuts is None or self.protected: return self newndata = len(cuts) try: newsigma = self.sigma.loc[cuts] except KeyError as e: # This will be an ugly erorr msg, but it should be scary anyway log.error(f"Problem applying cuts to {self.metadata}") raise e return dataclasses.replace(self, ndata=newndata, sigma=newsigma)
@property def luminosity_mapping(self): """Return the flavour combinations that contribute to the fktable in the form of a single array The return shape is: (nbasis,) for DIS (nbasis*2,) for hadronic """ basis = self.sigma.columns.to_numpy() if self.hadronic: ret = np.zeros(14 * 14, dtype=bool) ret[basis] = True basis = np.array(np.where(ret.reshape(14, 14))).T.reshape(-1) return basis
[docs] def get_np_fktable(self): """Returns the fktable as a dense numpy array that can be directly manipulated with numpy The return shape is: (ndata, nx, nbasis) for DIS (ndata, nx, nx, nbasis) for hadronic where nx is the length of the xgrid and nbasis the number of flavour contributions that contribute """ # Read up the shape of the output table ndata = self.ndata nx = len(self.xgrid) nbasis = self.sigma.shape[1] if ndata == 0: if self.hadronic: return np.zeros((ndata, nbasis, nx, nx)) return np.zeros((ndata, nbasis, nx)) # Make the dataframe into a dense numpy array # First get the data index out of the way # this is necessary because cuts/shifts and for performance reasons # otherwise we will be putting things in a numpy array in very awkward orders ns = self.sigma.unstack(level=("data",), fill_value=0) x1 = ns.index.get_level_values(0) if self.hadronic: x2 = ns.index.get_level_values(1) fk_raw = np.zeros((nx, nx, ns.shape[1])) fk_raw[x2, x1, :] = ns.values # The output is (ndata, basis, x1, x2) fktable = fk_raw.reshape((nx, nx, nbasis, ndata)).T else: fk_raw = np.zeros((nx, ns.shape[1])) fk_raw[x1, :] = ns.values # The output is (ndata, basis, x1) fktable = fk_raw.reshape((nx, nbasis, ndata)).T return fktable
[docs]@dataclasses.dataclass(eq=False) class CFactorData: """ Data contained in a CFactor Parameters ---------- description : str Information on how the data was obtained. central_value : array, shape(ndata) The value of the cfactor for each data point. uncertainty : array, shape(ndata) The absolute uncertainty on the cfactor if available. """ description: str central_value: np.array uncertainty: np.array
[docs]@dataclasses.dataclass(eq=False) class CommonData: """ Data contained in Commondata files, relevant cuts applied. Parameters ---------- setname : str Name of the dataset ndata : int Number of data points commondataproc : str Process type, one of 21 options nkin : int Number of kinematics specified nsys : int Number of systematics commondata_table : pd.DataFrame Pandas dataframe containing the commondata systype_table : pd.DataFrame Pandas dataframe containing the systype index for each systematic alongside the uncertainty type (ADD/MULT/RAND) and name (CORR/UNCORR/THEORYCORR/SKIP) systematics_table: pd.DataFrame Panda dataframe containing the table of systematics """ setname: str ndata: int commondataproc: str nkin: int nsys: int commondata_table: pd.DataFrame = dataclasses.field(repr=False) systype_table: pd.DataFrame = dataclasses.field(repr=False) legacy: bool systematics_table: Optional[pd.DataFrame] = dataclasses.field(init=None, repr=False) legacy_name: Optional[str] = None kin_variables: Optional[list] = None def __post_init__(self): self.systematics_table = self.commondata_table.drop( columns=["process", "data", "stat"] + KIN_NAMES ) if self.legacy_name is None: self.legacy_name = self.setname # TODO: set for now commondataproc as a string as well self.commondataproc = str(self.commondataproc)
[docs] def with_cuts(self, cuts): """A method to return a CommonData object where an integer mask has been applied, keeping only data points which pass cuts. Note if the first data point passes cuts, the first entry of ``cuts`` should be ``0``. Paramters --------- cuts: list or validphys.core.Cuts or None """ # Ensure that the cuts we're applying applies to this dataset # only check, however, if the cuts is of type :py:class:`validphys.core.Cuts` if hasattr(cuts, "name") and self.setname != raise ValueError( f"The cuts provided are for {} which does not apply " f"to this commondata file: {self.setname}" ) if hasattr(cuts, "load"): cuts = cuts.load() if cuts is None: return self # We must shift the cuts up by 1 since a cut of 0 implies the first data point # while commondata indexing starts at 1. cuts = list(map(lambda x: x + 1, cuts)) newndata = len(cuts) new_commondata_table = self.commondata_table.loc[cuts] return dataclasses.replace(self, ndata=newndata, commondata_table=new_commondata_table)
@property def kinematics(self): return self.commondata_table[KIN_NAMES]
[docs] def get_kintable(self): return self.kinematics.values
@property def central_values(self): return self.commondata_table["data"]
[docs] def with_central_value(self, cv): tb = self.commondata_table.copy() tb["data"] = cv return dataclasses.replace(self, commondata_table=tb)
[docs] def get_cv(self): return self.central_values.values
@property def stat_errors(self): return self.commondata_table["stat"] @property def multiplicative_errors(self): """Returns the systematics which are multiplicative (systype is MULT) in a percentage format, with SKIP uncertainties removed. """ mult_systype = self.systype_table[self.systype_table["treatment"] == "MULT"] mult_table = self.systematics_table.filter(like="MULT") if self.legacy: # Needed in legacy because every uncertainty appears as both mult and add # so it is necessary to select the uncertainties that are to be consireded as MULT/ADD # Minus 1 because iloc starts from 0, while the systype counting starts from 1 mult_table = mult_table.iloc[:, mult_systype.index - 1] mult_table.columns = mult_systype["name"].to_numpy() return mult_table.loc[:, mult_table.columns != "SKIP"] @property def additive_errors(self): """Returns the systematics which are additive (systype is ADD) as absolute uncertainties (same units as data), with SKIP uncertainties removed. """ add_systype = self.systype_table[self.systype_table["treatment"] == "ADD"] add_table = self.systematics_table.filter(like="ADD") if self.legacy: # Minus 1 because iloc starts from 0, while the systype counting starts from 1 add_table = add_table.iloc[:, add_systype.index - 1] add_table.columns = add_systype["name"].to_numpy() return add_table.loc[:, add_table.columns != "SKIP"]
[docs] def systematic_errors(self, central_values=None): """Returns all systematic errors as absolute uncertainties, with a single column for each uncertainty. Converts :py:attr:`multiplicative_errors` to units of data and then appends onto :py:attr:`additive_errors`. By default uses the experimental central values to perform conversion, but the user can supply a 1-D array of central values, with length :py:attr:`self.ndata`, to use instead of the experimental central values to calculate the absolute contribution of the multiplicative systematics. Parameters ---------- central_values: None, np.array 1-D array containing alternative central values to combine with multiplicative uncertainties. This array must have length equal to :py:attr:`self.ndata`. By default ``central_values`` is None, and the central values of the commondata are used. Returns ------- systematic_errors: pd.DataFrame Dataframe containing systematic errors. """ if central_values is None: central_values = self.central_values.to_numpy() converted_mult_errors = self.multiplicative_errors * central_values[:, np.newaxis] / 100 return pd.concat((self.additive_errors, converted_mult_errors), axis=1)
[docs] def export_data(self, buffer): """Exports the central data defined by this commondata instance to the given buffer""" ret = {"data_central": self.central_values.tolist()} yaml.safe_dump(ret, buffer)
[docs] def export_uncertainties(self, buffer): """Exports the uncertainties defined by this commondata instance to the given buffer""" definitions = {} for idx, row in self.systype_table.iterrows(): if row["name"] != "SKIP": definitions[f"sys_{idx}"] = {"treatment": row["treatment"], "type": row["name"]} # Order the definitions by treatment as ADD, MULT # TODO: make it so that it corresponds to the original order exactly sorted_definitions = { k: v for k, v in sorted(definitions.items(), key=lambda item: item[1]["treatment"]) } bins = [] for idx, row in self.systematic_errors().iterrows(): tmp = {"stat": float(self.stat_errors[idx])} # Hope things come in the right order... for key_name, val in zip(sorted_definitions, row): tmp[key_name] = float(val) bins.append(tmp) sorted_definitions["stat"] = { "description": "Uncorrelated statistical uncertainties", "treatment": "ADD", "type": "UNCORR", } ret = {"definitions": sorted_definitions, "bins": bins} yaml.safe_dump(ret, buffer)
[docs] def export(self, folder_path): """Wrapper around export_data and export_uncertainties to write both uncertainties and data after filtering to a given folder """ folder_path.mkdir(exist_ok=True) # Get the same names as one would use for the filters data_path, unc_path = generate_path_filtered_data(folder_path, self.setname) # And attach it to the given folder data_path = folder_path / unc_path = folder_path / # Export data and uncertainties self.export_data("w", encoding="utf-8")) self.export_uncertainties("w", encoding="utf-8")) return data_path, unc_path