Source code for contur.data.build_covariance

import numpy as np
import contur
import contur.config.config as cfg
import contur.data.static_db as cdb
import yoda

[docs] class CovarianceBuilder(object): """ `ao` Yoda AO apply_min: apply the minimum number of systematic uncertainties criteria when determining whether or not to use the error breakdown for correlations. Class to handle retrieval of annotations/errors from YODA objects """ def __init__(self, orig_ao, ao, aos, isScaled, scaleFactor, apply_min=True,ignore_corrs=False): self.ao=ao self.hasBreakdown=self._getBreakdownAttr(orig_ao,apply_min) self.readMatrix =self._getMatrixAttr() self.nbins=self.ao.numPoints() self.cov=None self.uncov=None self.covM=None self.diagonal=None self.errorBreakdown = None self.scaleFactor = scaleFactor if self.hasBreakdown: if 'Estimate0D' in orig_ao.type(): self.errorBreakdown = dict([ (src, no.array([ orig_ao.errEnv(src) ])) for src in ao.sources() ]) else: self.errorBreakdown = dict([ (src, np.array([ b.errEnv(src) if b.hasSource(src) \ else 0.0 for b in orig_ao.bins() ])) \ for src in orig_ao.sources() ]) # always build the diagonal matrix. self.buildCovFromErrorBar() if self.readMatrix: # read the covariance matrix from the dictionary of analysis objects try: self.read_cov_matrix(isScaled,aos) except Exception as ex: cfg.contur_log.warning(ex) self.readMatrix = False #print("RM",orig_ao.path(),self.covM) else: if ('Estimate' in orig_ao.type() and orig_ao.dim() > 1): self.buildCovFromBreakdown(orig_ao,ignore_corrs) #print("BD",orig_ao.path(),self.covM) else: self.covM = self.diagonal.copy() #print("EB",orig_ao.path(),self.covM) def _getBreakdownAttr(self,ao,apply_min): """ return true if this AO has an error breakdown """ if not 'Estimate' in ao.type(): cfg.contur_log.debug("{} is of type {}. No breakdown".format(ao.path(),ao.type())) return False if apply_min and (ao.numErrs() if 'Estimate0D' in ao.type() \ else min(b.numErrs() for b in ao.bins())) < cfg.min_num_sys: return False return True def _getMatrixAttr(self): """ return true if this AO has a covariance matrix stored in another AO. """ if cfg.diag: return False self._covname = cdb.get_covariance_name(self.ao.path()) if self._covname: return True else: self._corrname = cdb.get_correlation_name(self.ao.path()) if self._corrname: return True return False
[docs] def read_cov_matrix(self,isScaled,aos): """ read the covariance matrix from another AO and return it. """ if not self.readMatrix: self.covM = None return # if it has alreader been read, don't do it again. if self.covM is not None: return if self._covname: is_cov = True name = self._covname cfg.contur_log.debug("reading covariance matrix {}".format(self._covname)) else: is_cov = False name = self._corrname cfg.contur_log.debug("reading correlation matrix {}".format(self._corrname)) # read the covariance matrix into an array. matrix_ao = yoda.plotting.utils.mkPlotFriendlyScatter(aos[name],includeOverflows=False,includeMaskedBins=False) try: # take the number of bins from the measurement nbins = self.ao.numPoints() nbins2 = int(np.sqrt(matrix_ao.numPoints())) # check this is consistent with the matrix if nbins != nbins2: raise cfg.ConturError("Inconsistent number of entries ({} vs {}) in cov matrix: {}. Will not use it.".format(nbins,nbins2,self._covname)) self.covM = np.zeros((nbins,nbins)) i = 0 j = 0 for z in matrix_ao.zVals(): self.covM[i][j] = z i=i+1 if i==nbins: i=0 j=j+1 except: cfg.contur_log.error("Failed to read {}".format(self._covname)) raise if not is_cov: # need to pre and post muliply by uncertainties cfg.contur_log.debug("Converting correlation to covariance. {}".format(self.covM)) yErrs = np.diag(self.ao.yErrAvgs()) self.covM = np.dot(yErrs, np.dot(self.covM, yErrs)) cfg.contur_log.debug("resulting matrix is: {}".format(self.covM)) elif isScaled: # need to scale the covariance. self.covM = self.covM * self.scaleFactor * self.scaleFactor
[docs] def buildCovFromBreakdown(self,orig_ao,ignore_corrs): """ Get the covariance, calculated by YODA from the error breakdown """ self.covM = np.array(orig_ao.covarianceMatrix(ignore_corrs,False,False,"^stat|^uncor|Data stat"))
# do this when yoda python name is fixed. # = np.array(orig_ao.covarianceMatrix(ignoreOffdiagonalTerms=ignore_corrs))
[docs] def buildCovFromErrorBar(self): """ Build the diagonal covariance from error bars """ dummyM = np.outer(range(self.nbins), range(self.nbins)) self.diagonal = np.zeros(dummyM.shape) systErrs = np.zeros(self.nbins) for ibin in range(self.nbins): #symmetrize the errors (be conservative - use the largest!) systErrs[ibin] = max(self.ao.point(ibin).errs(self.ao.dim()-1)) # Note that yoda will take the average (below) when it computes a covariance matrix from the error breakdown. # There can also be round errors depending on the precision of the error breakdown. #systErrs[ibin] = (abs(self.ao.points()[ibin].yErrs()[0])+abs(self.ao.points()[ibin].yErrs()[1]))*0.5 self.diagonal += np.diag(systErrs * systErrs)
[docs] def getErrorBreakdown(self): """ return the breakdown of (symmetrised) uncertainties """ return self.errorBreakdown