From e953892e406811ebb7a128727c06e607b69b50eb Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Wed, 10 May 2023 16:39:17 -0600 Subject: [PATCH 01/10] subdomain basic stats --- .../SubdomainBasicStatistics.py | 1565 +++++++++++++++++ 1 file changed, 1565 insertions(+) create mode 100644 ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py diff --git a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py new file mode 100644 index 0000000000..b3b5e1d40f --- /dev/null +++ b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py @@ -0,0 +1,1565 @@ +# Copyright 2017 Battelle Energy Alliance, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on July 10, 2013 + +@author: aalfonsi +""" +#External Modules--------------------------------------------------------------- +import numpy as np +import os +import copy +from collections import OrderedDict, defaultdict +import six +import xarray as xr +import scipy.stats as stats +#External Modules End----------------------------------------------------------- + +#Internal Modules--------------------------------------------------------------- +from .PostProcessorInterface import PostProcessorInterface +from .BasicStatistics import BasicStatistics +from ...utils import utils +from ...utils import InputData, InputTypes +from ...utils import mathUtils +from ... import Files +#Internal Modules End----------------------------------------------------------- + +class SubdomainBasicStatistics(PostProcessorInterface): + """ + BasicStatistics filter class. It computes all the most popular statistics + """ + + scalarVals = ['expectedValue', + 'minimum', + 'maximum', + 'median', + 'variance', + 'sigma', + 'percentile', + 'variationCoefficient', + 'skewness', + 'kurtosis', + 'samples', + 'higherPartialVariance', # Statistic metric not available yet + 'higherPartialSigma', # Statistic metric not available yet + 'lowerPartialSigma', # Statistic metric not available yet + 'lowerPartialVariance' # Statistic metric not available yet + ] + vectorVals = ['sensitivity', + 'covariance', + 'pearson', + 'spearman', + 'NormalizedSensitivity', + 'VarianceDependentSensitivity'] + # quantities that the standard error can be computed + steVals = ['expectedValue_ste', + 'median_ste', + 'variance_ste', + 'sigma_ste', + 'skewness_ste', + 'kurtosis_ste', + 'percentile_ste'] + + @classmethod + def getInputSpecification(cls): + """ + Method to get a reference to a class that specifies the input data for + class cls. + @ In, cls, the class for which we are retrieving the specification + @ Out, inputSpecification, InputData.ParameterInput, class to use for + specifying input of cls. + """ + ## This will replace the lines above + inputSpecification = super().getInputSpecification() + + for scalar in cls.scalarVals: + scalarSpecification = InputData.parameterInputFactory(scalar, contentType=InputTypes.StringListType) + if scalar == 'percentile': + #percent is a string type because otherwise we can't tell 95.0 from 95 + # which matters because the number is used in output. + scalarSpecification.addParam("percent", InputTypes.StringListType) + # percentile has additional "interpolation" parameter + scalarSpecification.addParam("interpolation", + param_type=InputTypes.makeEnumType("interpolation", + "interpolationType", + ["linear", "midpoint"]), + default="linear", + descr="""Interpolation method for percentile calculation. + 'linear' uses linear interpolation between nearest + data points while 'midpoint' uses the average of the + nearest data points.""") + if scalar == 'median': + # median has additional "interpolation" parameter + scalarSpecification.addParam("interpolation", + param_type=InputTypes.makeEnumType("interpolation", + "interpolationType", + ["linear", "midpoint"]), + default="linear", + descr="""Interpolation method for median calculation. 'linear' + uses linear interpolation between nearest data points + while 'midpoint' uses the average of the nearest data + points.""") + scalarSpecification.addParam("prefix", InputTypes.StringType) + inputSpecification.addSub(scalarSpecification) + + for vector in cls.vectorVals: + vectorSpecification = InputData.parameterInputFactory(vector) + vectorSpecification.addParam("prefix", InputTypes.StringType) + features = InputData.parameterInputFactory('features', + contentType=InputTypes.StringListType) + vectorSpecification.addSub(features) + targets = InputData.parameterInputFactory('targets', + contentType=InputTypes.StringListType) + vectorSpecification.addSub(targets) + inputSpecification.addSub(vectorSpecification) + + pivotParameterInput = InputData.parameterInputFactory('pivotParameter', contentType=InputTypes.StringType) + inputSpecification.addSub(pivotParameterInput) + + datasetInput = InputData.parameterInputFactory('dataset', contentType=InputTypes.BoolType) + inputSpecification.addSub(datasetInput) + + methodsToRunInput = InputData.parameterInputFactory("methodsToRun", contentType=InputTypes.StringType) + inputSpecification.addSub(methodsToRunInput) + + biasedInput = InputData.parameterInputFactory("biased", contentType=InputTypes.BoolType) + inputSpecification.addSub(biasedInput) + + multipleFeaturesInput = InputData.parameterInputFactory("multipleFeatures", contentType=InputTypes.BoolType) + inputSpecification.addSub(multipleFeaturesInput) + + return inputSpecification + + def __init__(self): + """ + Constructor + @ In, None + @ Out, None + """ + super().__init__() + self.parameters = {} # parameters dictionary (they are basically stored into a dictionary identified by tag "targets" + self.acceptedCalcParam = self.scalarVals + self.vectorVals + self.what = self.acceptedCalcParam # what needs to be computed... default...all + self.methodsToRun = [] # if a function is present, its outcome name is here stored... if it matches one of the known outcomes, the pp is going to use the function to compute it + self.printTag = 'PostProcessor BASIC STATISTIC' + self.biased = False # biased statistics? + self.pivotParameter = None # time-dependent statistics pivot parameter + self.pivotValue = None # time-dependent statistics pivot parameter values + self.dynamic = False # is it time-dependent? + self.sampleTag = None # Tag used to track samples + self.pbPresent = False # True if the ProbabilityWeight is available + self.realizationWeight = None # The joint probabilities + self.steMetaIndex = 'targets' # when Dataset is requested as output, the default index of ste metadata is ['targets', self.pivotParameter] + self.multipleFeatures = True # True if multiple features are employed in linear regression as feature inputs + self.sampleSize = None # number of sample size + self.calculations = {} + self.validDataType = ['PointSet', 'HistorySet', 'DataSet'] # The list of accepted types of DataObject + + def inputToInternal(self, currentInp): + """ + Method to convert an input object into the internal format that is + understandable by this pp. + @ In, currentInp, object, an object that needs to be converted + @ Out, (inputDataset, pbWeights), tuple, the dataset of inputs and the corresponding variable probability weight + """ + # The BasicStatistics postprocessor only accept DataObjects + self.dynamic = False + currentInput = currentInp [-1] if type(currentInp) == list else currentInp + if len(currentInput) == 0: + self.raiseAnError(IOError, "In post-processor " +self.name+" the input "+currentInput.name+" is empty.") + + pbWeights = None + if type(currentInput).__name__ == 'tuple': + return currentInput + # TODO: convert dict to dataset, I think this will be removed when DataSet is used by other entities that + # are currently using this Basic Statisitics PostProcessor. + if type(currentInput).__name__ == 'dict': + if 'targets' not in currentInput.keys(): + self.raiseAnError(IOError, 'Did not find targets in the input dictionary') + inputDataset = xr.Dataset() + for var, val in currentInput['targets'].items(): + inputDataset[var] = val + if 'metadata' in currentInput.keys(): + metadata = currentInput['metadata'] + self.pbPresent = True if 'ProbabilityWeight' in metadata else False + if self.pbPresent: + pbWeights = xr.Dataset() + self.realizationWeight = xr.Dataset() + self.realizationWeight['ProbabilityWeight'] = metadata['ProbabilityWeight']/metadata['ProbabilityWeight'].sum() + for target in self.parameters['targets']: + pbName = 'ProbabilityWeight-' + target + if pbName in metadata: + pbWeights[target] = metadata[pbName]/metadata[pbName].sum() + elif self.pbPresent: + pbWeights[target] = self.realizationWeight['ProbabilityWeight'] + else: + self.raiseAWarning('BasicStatistics postprocessor did not detect ProbabilityWeights! Assuming unit weights instead...') + else: + self.raiseAWarning('BasicStatistics postprocessor did not detect ProbabilityWeights! Assuming unit weights instead...') + if 'RAVEN_sample_ID' not in inputDataset.sizes.keys(): + self.raiseAWarning('BasicStatisitics postprocessor did not detect RAVEN_sample_ID! Assuming the first dimension of given data...') + self.sampleTag = utils.first(inputDataset.sizes.keys()) + return inputDataset, pbWeights + + if currentInput.type not in ['PointSet','HistorySet']: + self.raiseAnError(IOError, self, 'BasicStatistics postprocessor accepts PointSet and HistorySet only! Got ' + currentInput.type) + + # extract all required data from input DataObjects, an input dataset is constructed + dataSet = currentInput.asDataset() + try: + inputDataset = dataSet[self.parameters['targets']] + except KeyError: + missing = [var for var in self.parameters['targets'] if var not in dataSet] + self.raiseAnError(KeyError, "Variables: '{}' missing from dataset '{}'!".format(", ".join(missing),currentInput.name)) + self.sampleTag = currentInput.sampleTag + + if currentInput.type == 'HistorySet': + dims = inputDataset.sizes.keys() + if self.pivotParameter is None: + if len(dims) > 1: + self.raiseAnError(IOError, self, 'Time-dependent statistics is requested (HistorySet) but no pivotParameter \ + got inputted!') + elif self.pivotParameter not in dims: + self.raiseAnError(IOError, self, 'Pivot parameter', self.pivotParameter, 'is not the associated index for \ + requested variables', ','.join(self.parameters['targets'])) + else: + self.dynamic = True + if not currentInput.checkIndexAlignment(indexesToCheck=self.pivotParameter): + self.raiseAnError(IOError, "The data provided by the data objects", currentInput.name, "is not synchronized!") + self.pivotValue = inputDataset[self.pivotParameter].values + if self.pivotValue.size != len(inputDataset.groupby(self.pivotParameter)): + msg = "Duplicated values were identified in pivot parameter, please use the 'HistorySetSync'" + \ + " PostProcessor to syncronize your data before running 'BasicStatistics' PostProcessor." + self.raiseAnError(IOError, msg) + # extract all required meta data + metaVars = currentInput.getVars('meta') + self.pbPresent = True if 'ProbabilityWeight' in metaVars else False + if self.pbPresent: + pbWeights = xr.Dataset() + self.realizationWeight = dataSet[['ProbabilityWeight']]/dataSet[['ProbabilityWeight']].sum() + for target in self.parameters['targets']: + pbName = 'ProbabilityWeight-' + target + if pbName in metaVars: + pbWeights[target] = dataSet[pbName]/dataSet[pbName].sum() + elif self.pbPresent: + pbWeights[target] = self.realizationWeight['ProbabilityWeight'] + else: + self.raiseAWarning('BasicStatistics postprocessor did not detect ProbabilityWeights! Assuming unit weights instead...') + + return inputDataset, pbWeights + + def initialize(self, runInfo, inputs, initDict): + """ + Method to initialize the BasicStatistic pp. In here the working dir is + grepped. + @ In, runInfo, dict, dictionary of run info (e.g. working dir, etc) + @ In, inputs, list, list of inputs + @ In, initDict, dict, dictionary with initialization options + @ Out, None + """ + #construct a list of all the parameters that have requested values into self.allUsedParams + self.allUsedParams = set() + for metricName in self.scalarVals + self.vectorVals: + if metricName in self.toDo.keys(): + for entry in self.toDo[metricName]: + self.allUsedParams.update(entry['targets']) + try: + self.allUsedParams.update(entry['features']) + except KeyError: + pass + + #for backward compatibility, compile the full list of parameters used in Basic Statistics calculations + self.parameters['targets'] = list(self.allUsedParams) + super().initialize(runInfo, inputs, initDict) + inputObj = inputs[-1] if type(inputs) == list else inputs + if inputObj.type == 'HistorySet': + self.dynamic = True + inputMetaKeys = [] + outputMetaKeys = [] + for metric, infos in self.toDo.items(): + if metric in self.scalarVals + self.vectorVals: + steMetric = metric + '_ste' + if steMetric in self.steVals: + for info in infos: + prefix = info['prefix'] + for target in info['targets']: + if metric == 'percentile': + for strPercent in info['strPercent']: + metaVar = prefix + '_' + strPercent + '_ste_' + target if not self.outputDataset else metric + '_ste' + metaDim = inputObj.getDimensions(target) + if len(metaDim[target]) == 0: + inputMetaKeys.append(metaVar) + else: + outputMetaKeys.append(metaVar) + else: + metaVar = prefix + '_ste_' + target if not self.outputDataset else metric + '_ste' + metaDim = inputObj.getDimensions(target) + if len(metaDim[target]) == 0: + inputMetaKeys.append(metaVar) + else: + outputMetaKeys.append(metaVar) + metaParams = {} + if not self.outputDataset: + if len(outputMetaKeys) > 0: + metaParams = {key:[self.pivotParameter] for key in outputMetaKeys} + else: + if len(outputMetaKeys) > 0: + params = {} + for key in outputMetaKeys + inputMetaKeys: + # percentile standard error has additional index + if key == 'percentile_ste': + params[key] = [self.pivotParameter, self.steMetaIndex, 'percent'] + else: + params[key] = [self.pivotParameter, self.steMetaIndex] + metaParams.update(params) + elif len(inputMetaKeys) > 0: + params = {} + for key in inputMetaKeys: + # percentile standard error has additional index + if key == 'percentile_ste': + params[key] = [self.steMetaIndex, 'percent'] + else: + params[key] = [self.steMetaIndex] + metaParams.update(params) + metaKeys = inputMetaKeys + outputMetaKeys + self.addMetaKeys(metaKeys,metaParams) + + def _handleInput(self, paramInput, childVals=None): + """ + Function to handle the parsed paramInput for this class. + @ In, paramInput, ParameterInput, the already parsed input. + @ In, childVals, list, optional, quantities requested from child statistical object + @ Out, None + """ + if childVals is None: + childVals = [] + self.toDo = {} + for child in paramInput.subparts: + tag = child.getName() + #because percentile is strange (has attached parameters), we address it first + if tag in self.scalarVals + self.vectorVals: + if 'prefix' not in child.parameterValues: + self.raiseAnError(IOError, "No prefix is provided for node: ", tag) + #get the prefix + prefix = child.parameterValues['prefix'] + if tag == 'percentile': + #get targets + targets = set(child.value) + #what if user didn't give any targets? + if len(targets)<1: + self.raiseAWarning('No targets were specified in text of <'+tag+'>! Skipping metric...') + continue + #prepare storage dictionary, keys are percentiles, values are set(targets) + if tag not in self.toDo.keys(): + self.toDo[tag] = [] # list of {'targets':(), 'prefix':str, 'percent':str, 'interpolation':str} + if 'percent' not in child.parameterValues: + reqPercent = [0.05, 0.95] + strPercent = ['5','95'] + else: + reqPercent = set(utils.floatConversion(percent)/100. for percent in child.parameterValues['percent']) + strPercent = set(percent for percent in child.parameterValues['percent']) + if 'interpolation' not in child.parameterValues: + interpolation = 'linear' + else: + interpolation = child.parameterValues['interpolation'] + self.toDo[tag].append({'targets':set(targets), + 'prefix':prefix, + 'percent':reqPercent, + 'strPercent':strPercent, + 'interpolation':interpolation}) + # median also has an attached parameter + elif tag == 'median': + if tag not in self.toDo.keys(): + self.toDo[tag] = [] # list of {'targets':{}, 'prefix':str, 'interpolation':str} + if 'interpolation' not in child.parameterValues: + interpolation = 'linear' + else: + interpolation = child.parameterValues['interpolation'] + self.toDo[tag].append({'targets':set(child.value), + 'prefix':prefix, + 'interpolation':interpolation}) + elif tag in self.scalarVals: + if tag not in self.toDo.keys(): + self.toDo[tag] = [] # list of {'targets':(), 'prefix':str} + self.toDo[tag].append({'targets':set(child.value), + 'prefix':prefix}) + elif tag in self.vectorVals: + if tag not in self.toDo.keys(): + self.toDo[tag] = [] # list of {'targets':(),'features':(), 'prefix':str} + tnode = child.findFirst('targets') + if tnode is None: + self.raiseAnError('Request for vector value <'+tag+'> requires a "targets" node, and none was found!') + fnode = child.findFirst('features') + if fnode is None: + self.raiseAnError('Request for vector value <'+tag+'> requires a "features" node, and none was found!') + # we're storing toDo[tag] as a list of dictionaries. This is because the user might specify multiple + # nodes with the same metric (tag), but with different targets and features. For instance, the user might + # want the sensitivity of A and B to X and Y, and the sensitivity of C to W and Z, but not the sensitivity + # of A to W. If we didn't keep them separate, we could potentially waste a fair number of calculations. + self.toDo[tag].append({'targets':set(tnode.value), + 'features':set(fnode.value), + 'prefix':prefix}) + elif tag == "biased": + self.biased = child.value + elif tag == "pivotParameter": + self.pivotParameter = child.value + elif tag == "dataset": + self.outputDataset = child.value + elif tag == "multipleFeatures": + self.multipleFeatures = child.value + else: + if tag not in childVals: + self.raiseAWarning('Unrecognized node in BasicStatistics "',tag,'" has been ignored!') + + assert (len(self.toDo)>0), self.raiseAnError(IOError, 'BasicStatistics needs parameters to work on! Please check input for PP: ' + self.name) + + def _computePower(self, p, dataset): + """ + Compute the p-th power of weights + @ In, p, int, the power + @ In, dataset, xarray.Dataset, probability weights of all input variables + @ Out, pw, xarray.Dataset, the p-th power of weights + """ + pw = {} + coords = dataset.coords + for target, targValue in dataset.variables.items(): + ##remove index variable + if target in coords: + continue + pw[target] = np.power(targValue,p) + pw = xr.Dataset(data_vars=pw,coords=coords) + return pw + + def __computeVp(self,p,weights): + """ + Compute the sum of p-th power of weights + @ In, p, int, the power + @ In, weights, xarray.Dataset, probability weights of all input variables + @ Out, vp, xarray.Dataset, the sum of p-th power of weights + """ + vp = self._computePower(p,weights) + vp = vp.sum() + return vp + + def __computeEquivalentSampleSize(self,weights): + """ + Compute the equivalent sample size for given probability weights + @ In, weights, xarray.Dataset, probability weights of all input variables + @ Out, equivalentSize, xarray.Dataset, the equivalent sample size + """ + # The equivalent sample size for given samples, i.e. (sum of weights) squared / sum of the squared weights + # The definition of this quantity can be found: + # R. F. Potthoff, M. A. Woodbury and K. G. Manton, "'Equivalent SampleSize' and 'Equivalent Degrees of Freedom' + # Refinements for Inference Using Survey Weights Under Superpopulation Models", Journal of the American Statistical + # Association, Vol. 87, No. 418 (1992) + v1Square = self.__computeVp(1,weights)**2 + v2 = self.__computeVp(2,weights) + equivalentSize = v1Square/v2 + return equivalentSize + + def __computeUnbiasedCorrection(self,order,weightsOrN): + """ + Compute unbiased correction given weights and momement order + Reference paper: + Lorenzo Rimoldini, "Weighted skewness and kurtosis unbiased by sample size", http://arxiv.org/pdf/1304.6564.pdf + @ In, order, int, moment order + @ In, weightsOrN, xarray.Dataset or int, if xarray.Dataset -> weights else -> number of samples + @ Out, corrFactor, xarray.Dataset or int, xarray.Dataset (order <=3) or tuple of xarray.Dataset (order ==4), + the unbiased correction factor if weightsOrN is xarray.Dataset else integer + """ + if order > 4: + self.raiseAnError(RuntimeError,"computeUnbiasedCorrection is implemented for order <=4 only!") + if type(weightsOrN).__name__ not in ['int','int8','int16','int64','int32']: + if order == 2: + V1, v1Square, V2 = self.__computeVp(1, weightsOrN), self.__computeVp(1, weightsOrN)**2.0, self.__computeVp(2, weightsOrN) + corrFactor = v1Square/(v1Square-V2) + elif order == 3: + V1, v1Cubic, V2, V3 = self.__computeVp(1, weightsOrN), self.__computeVp(1, weightsOrN)**3.0, self.__computeVp(2, weightsOrN), self.__computeVp(3, weightsOrN) + corrFactor = v1Cubic/(v1Cubic-3.0*V2*V1+2.0*V3) + elif order == 4: + V1, v1Square, V2, V3, V4 = self.__computeVp(1, weightsOrN), self.__computeVp(1, weightsOrN)**2.0, self.__computeVp(2, weightsOrN), self.__computeVp(3, weightsOrN), self.__computeVp(4, weightsOrN) + numer1 = v1Square*(v1Square**2.0-3.0*v1Square*V2+2.0*V1*V3+3.0*V2**2.0-3.0*V4) + numer2 = 3.0*v1Square*(2.0*v1Square*V2-2.0*V1*V3-3.0*V2**2.0+3.0*V4) + denom = (v1Square-V2)*(v1Square**2.0-6.0*v1Square*V2+8.0*V1*V3+3.0*V2**2.0-6.0*V4) + corrFactor = numer1/denom ,numer2/denom + else: + if order == 2: + corrFactor = float(weightsOrN)/(float(weightsOrN)-1.0) + elif order == 3: + corrFactor = (float(weightsOrN)**2.0)/((float(weightsOrN)-1)*(float(weightsOrN)-2)) + elif order == 4: + corrFactor = (float(weightsOrN)*(float(weightsOrN)**2.0-2.0*float(weightsOrN)+3.0))/((float(weightsOrN)-1)*(float(weightsOrN)-2)*(float(weightsOrN)-3)),(3.0*float(weightsOrN)*(2.0*float(weightsOrN)-3.0))/((float(weightsOrN)-1)*(float(weightsOrN)-2)*(float(weightsOrN)-3)) + return corrFactor + + def _computeKurtosis(self, arrayIn, expValue, variance, pbWeight=None, dim=None): + """ + Method to compute the Kurtosis (fisher) of an array of observations + @ In, arrayIn, xarray.Dataset, the dataset from which the Kurtosis needs to be estimated + @ In, expValue, xarray.Dataset, expected value of arrayIn + @ In, variance, xarray.Dataset, variance of arrayIn + @ In, pbWeight, xarray.DataSet, optional, the reliability weights that correspond to the values in 'array'. + If not present, an unweighted approach is used + @ Out, result, xarray.Dataset, the Kurtosis of the dataset arrayIn. + """ + if dim is None: + dim = self.sampleTag + vr = self._computePower(2.0, variance) + if pbWeight is not None: + unbiasCorr = self.__computeUnbiasedCorrection(4,pbWeight) if not self.biased else 1.0 + vp = 1.0/self.__computeVp(1,pbWeight) + p4 = ((arrayIn - expValue)**4.0 * pbWeight).sum(dim=dim) + if not self.biased: + p2 = ((arrayIn - expValue)**2.0 * pbWeight).sum(dim=dim) + result = -3.0 + (p4*unbiasCorr[0]*vp - (p2*vp)**2.0 * unbiasCorr[1]) / vr + else: + result = -3.0 + (p4 * vp * unbiasCorr) / vr + else: + unbiasCorr = self.__computeUnbiasedCorrection(4,int(arrayIn.sizes[dim])) if not self.biased else 1.0 + vp = 1.0 / arrayIn.sizes[dim] + p4 = ((arrayIn - expValue)**4.0).sum(dim=dim) + if not self.biased: + p2 = (arrayIn - expValue).var(dim=dim) + result = -3.0 + (p4*unbiasCorr[0]*vp-p2**2.0*unbiasCorr[1]) / vr + else: + result = -3.0 + (p4*unbiasCorr*vp) / vr + return result + + def _computeSkewness(self, arrayIn, expValue, variance, pbWeight=None, dim=None): + """ + Method to compute the skewness of an array of observations + @ In, arrayIn, xarray.Dataset, the dataset from which the skewness needs to be estimated + @ In, expValue, xarray.Dataset, expected value of arrayIn + @ In, variance, xarray.Dataset, variance value of arrayIn + @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. + If not present, an unweighted approach is used + @ Out, result, xarray.Dataset, the skewness of the dataset arrayIn + """ + if dim is None: + dim = self.sampleTag + vr = self._computePower(1.5, variance) + if pbWeight is not None: + unbiasCorr = self.__computeUnbiasedCorrection(3,pbWeight) if not self.biased else 1.0 + vp = 1.0/self.__computeVp(1,pbWeight) + result = ((arrayIn - expValue)**3 * pbWeight).sum(dim=dim) * vp * unbiasCorr / vr + else: + unbiasCorr = self.__computeUnbiasedCorrection(3,int(arrayIn.sizes[dim])) if not self.biased else 1.0 + vp = 1.0 / arrayIn.sizes[dim] + result = ((arrayIn - expValue)**3).sum(dim=dim) * vp * unbiasCorr / vr + return result + + def _computeVariance(self, arrayIn, expValue, pbWeight=None, dim = None): + """ + Method to compute the Variance (fisher) of an array of observations + @ In, arrayIn, xarray.Dataset, the dataset from which the Variance needs to be estimated + @ In, expValue, xarray.Dataset, expected value of arrayIn + @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. + If not present, an unweighted approach is used + @ Out, result, xarray.Dataset, the Variance of the dataset arrayIn + """ + if dim is None: + dim = self.sampleTag + if pbWeight is not None: + unbiasCorr = self.__computeUnbiasedCorrection(2,pbWeight) if not self.biased else 1.0 + vp = 1.0/self.__computeVp(1,pbWeight) + result = ((arrayIn-expValue)**2 * pbWeight).sum(dim=dim) * vp * unbiasCorr + else: + unbiasCorr = self.__computeUnbiasedCorrection(2,int(arrayIn.sizes[dim])) if not self.biased else 1.0 + result = (arrayIn-expValue).var(dim=dim) * unbiasCorr + return result + + def _computeLowerPartialVariance(self, arrayIn, medValue, pbWeight=None, dim = None): + """ + Method to compute the lower partial variance of an array of observations + @ In, arrayIn, xarray.Dataset, the dataset from which the Variance needs to be estimated + @ In, medValue, xarray.Dataset, median value of arrayIn + @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. + If not present, an unweighted approach is used + @ Out, result, xarray.Dataset, the lower partial variance of the dataset arrayIn + """ + if dim is None: + dim = self.sampleTag + diff = (medValue-arrayIn).clip(min=0) + if pbWeight is not None: + vp = 1.0/self.__computeVp(1,pbWeight) + result = ((diff)**2 * pbWeight).sum(dim=dim) * vp + else: + result = diff.var(dim=dim) + return result + + def _computeHigherPartialVariance(self, arrayIn, medValue, pbWeight=None, dim = None): + """ + Method to compute the higher partial variance of an array of observations + @ In, arrayIn, xarray.Dataset, the dataset from which the Variance needs to be estimated + @ In, medValue, xarray.Dataset, median value of arrayIn + @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. + If not present, an unweighted approach is used + @ Out, result, xarray.Dataset, the higher partial variance of the dataset arrayIn + """ + if dim is None: + dim = self.sampleTag + diff = (arrayIn-medValue).clip(min=0) + if pbWeight is not None: + vp = 1.0/self.__computeVp(1,pbWeight) + result = ((diff)**2 * pbWeight).sum(dim=dim) * vp + else: + result = diff.var(dim=dim) + return result + + def _computeSigma(self,arrayIn,variance,pbWeight=None): + """ + Method to compute the sigma of an array of observations + @ In, arrayIn, xarray.Dataset, the dataset from which the standard deviation needs to be estimated + @ In, variance, xarray.Dataset, variance of arrayIn + @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. + If not present, an unweighted approach is used + @ Out, sigma, xarray.Dataset, the sigma of the dataset of arrayIn + """ + return np.sqrt(variance) + + def _computeWeightedPercentile(self,arrayIn,pbWeight,interpolation='linear',percent=[0.5]): + """ + Method to compute the weighted percentile in a array of data + @ In, arrayIn, list/numpy.array, the array of values from which the percentile needs to be estimated + @ In, pbWeight, list/numpy.array, the reliability weights that correspond to the values in 'array' + @ In, interpolation, str, 'linear' or 'midpoint' + @ In, percent, list/numpy.array, the percentile(s) that needs to be computed (between 0.01 and 1.0) + @ Out, result, list, the percentile(s) + """ + + # only do the argsort once for all requested percentiles + idxs = np.argsort(np.asarray(list(zip(pbWeight,arrayIn)))[:,1]) + # Inserting [0.0,arrayIn[idxs[0]]] is needed when few samples are generated and + # a percentile that is < that the first pb weight is requested. Otherwise the median + # is returned. + sortedWeightsAndPoints = np.insert(np.asarray(list(zip(pbWeight[idxs],arrayIn[idxs]))),0,[0.0,arrayIn[idxs[0]]],axis=0) + weightsCDF = np.cumsum(sortedWeightsAndPoints[:,0]) + weightsCDF /= weightsCDF[-1] + if interpolation == 'linear': + result = np.interp(percent, weightsCDF, sortedWeightsAndPoints[:, 1]).tolist() + elif interpolation == 'midpoint': + result = [self._computeSingleWeightedPercentile(pct, weightsCDF, sortedWeightsAndPoints) for pct in percent] + + return result + + def _computeSingleWeightedPercentile(self, pct, weightsCDF, sortedWeightsAndPoints): + """ + Method to compute a single percentile + @ In, pct, float, the percentile + @ In, weightsCDF, numpy.array, the cumulative sum of weights (CDF) + @ In, sortedWeightsAndPoints, numpy.array, array of weights and data points + @ Out, result, float, the percentile + """ + + # This step returns the index of the array which is < than the percentile, because + # the insertion create another entry, this index should shift to the bigger side + indexL = utils.first(np.asarray(weightsCDF >= pct).nonzero())[0] + # This step returns the indices (list of index) of the array which is > than the percentile + indexH = utils.first(np.asarray(weightsCDF > pct).nonzero()) + try: + # if the indices exists that means the desired percentile lies between two data points + # with index as indexL and indexH[0]. Calculate the midpoint of these two points + result = 0.5*(sortedWeightsAndPoints[indexL,1]+sortedWeightsAndPoints[indexH[0],1]) + except IndexError: + result = sortedWeightsAndPoints[indexL,1] + + return result + + def __runLocal(self, inputData): + """ + This method executes the postprocessor action. In this case, it computes all the requested statistical FOMs + @ In, inputData, tuple, (inputDataset, pbWeights), tuple, the dataset of inputs and the corresponding + variable probability weight + @ Out, outputSet or outputDict, xarray.Dataset or dict, dataset or dictionary containing the results + """ + inputDataset, pbWeights = inputData[0], inputData[1] + #storage dictionary for skipped metrics + self.skipped = {} + #construct a dict of required computations + needed = dict((metric,{'targets':set(),'percent':set(),'interpolation':''}) for metric in self.scalarVals) + needed.update(dict((metric,{'targets':set(),'features':set()}) for metric in self.vectorVals)) + for metric, params in self.toDo.items(): + if metric in self.scalarVals + self.vectorVals: + for entry in params: + needed[metric]['targets'].update(entry['targets']) + try: + needed[metric]['features'].update(entry['features']) + except KeyError: + pass + try: + needed[metric]['percent'].update(entry['percent']) + except KeyError: + pass + try: + needed[metric]['interpolation'] = entry['interpolation'] + except KeyError: + pass + # variable | needs | needed for + # -------------------------------------------------------------------- + # skewness needs | expectedValue,variance | + # kurtosis needs | expectedValue,variance | + # median needs | | lowerPartialVariance, higherPartialVariance + # percentile needs | expectedValue,sigma | + # maximum needs | | + # minimum needs | | + # covariance needs | | pearson,VarianceDependentSensitivity,NormalizedSensitivity + # NormalizedSensitivity | covariance,VarDepSens | + # VarianceDependentSensitivity | covariance | NormalizedSensitivity + # sensitivity needs | | + # pearson needs | covariance | + # sigma needs | variance | variationCoefficient + # variance | expectedValue | sigma, skewness, kurtosis + # expectedValue | | variance, variationCoefficient, skewness, kurtosis, + # | | lowerPartialVariance, higherPartialVariance + # lowerPartialVariance needs | expectedValue,median | lowerPartialSigma + # lowerPartialSigma needs | lowerPartialVariance | + # higherPartialVariance needs | expectedValue,median | higherPartialSigma + # higherPartialSigma needs | higherPartialVariance | + + # update needed dictionary when standard errors are requested + needed['expectedValue']['targets'].update(needed['sigma']['targets']) + needed['expectedValue']['targets'].update(needed['variationCoefficient']['targets']) + needed['expectedValue']['targets'].update(needed['variance']['targets']) + needed['expectedValue']['targets'].update(needed['median']['targets']) + needed['expectedValue']['targets'].update(needed['skewness']['targets']) + needed['expectedValue']['targets'].update(needed['kurtosis']['targets']) + needed['expectedValue']['targets'].update(needed['NormalizedSensitivity']['targets']) + needed['expectedValue']['targets'].update(needed['NormalizedSensitivity']['features']) + needed['expectedValue']['targets'].update(needed['percentile']['targets']) + needed['sigma']['targets'].update(needed['expectedValue']['targets']) + needed['sigma']['targets'].update(needed['percentile']['targets']) + needed['variance']['targets'].update(needed['sigma']['targets']) + needed['lowerPartialVariance']['targets'].update(needed['lowerPartialSigma']['targets']) + needed['higherPartialVariance']['targets'].update(needed['higherPartialSigma']['targets']) + needed['median']['targets'].update(needed['lowerPartialVariance']['targets']) + needed['median']['targets'].update(needed['higherPartialVariance']['targets']) + needed['covariance']['targets'].update(needed['NormalizedSensitivity']['targets']) + needed['covariance']['features'].update(needed['NormalizedSensitivity']['features']) + needed['VarianceDependentSensitivity']['targets'].update(needed['NormalizedSensitivity']['targets']) + needed['VarianceDependentSensitivity']['features'].update(needed['NormalizedSensitivity']['features']) + needed['covariance']['targets'].update(needed['pearson']['targets']) + needed['covariance']['features'].update(needed['pearson']['features']) + needed['covariance']['targets'].update(needed['VarianceDependentSensitivity']['targets']) + needed['covariance']['features'].update(needed['VarianceDependentSensitivity']['features']) + + for metric, params in needed.items(): + needed[metric]['targets'] = list(params['targets']) + try: + needed[metric]['features'] = list(params['features']) + except KeyError: + pass + + # + # BEGIN actual calculations + # + + calculations = {} + + ################# + # SCALAR VALUES # + ################# + # + # samples + # + self.sampleSize = inputDataset.sizes[self.sampleTag] + metric = 'samples' + if len(needed[metric]['targets']) > 0: + self.raiseADebug('Starting "'+metric+'"...') + if self.dynamic: + nt = inputDataset.sizes[self.pivotParameter] + sampleMat = np.zeros((len(self.parameters['targets']), len(self.pivotValue))) + sampleMat.fill(self.sampleSize) + samplesDA = xr.DataArray(sampleMat,dims=('targets', self.pivotParameter), coords={'targets':self.parameters['targets'], self.pivotParameter:self.pivotValue}) + else: + sampleMat = np.zeros(len(self.parameters['targets'])) + sampleMat.fill(self.sampleSize) + samplesDA = xr.DataArray(sampleMat,dims=('targets'), coords={'targets':self.parameters['targets']}) + self.calculations[metric] = samplesDA + calculations[metric] = samplesDA + # + # expected value + # + metric = 'expectedValue' + if len(needed[metric]['targets']) > 0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + if self.pbPresent: + relWeight = pbWeights[list(needed[metric]['targets'])] + equivalentSize = self.__computeEquivalentSampleSize(relWeight) + dataSet = dataSet * relWeight + expectedValueDS = dataSet.sum(dim = self.sampleTag) + calculations['equivalentSamples'] = equivalentSize + else: + expectedValueDS = dataSet.mean(dim = self.sampleTag) + self.calculations[metric] = expectedValueDS + calculations[metric] = expectedValueDS + # + # variance + # + metric = 'variance' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + meanSet = calculations['expectedValue'][list(needed[metric]['targets'])] + relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None + varianceDS = self._computeVariance(dataSet,meanSet,pbWeight=relWeight,dim=self.sampleTag) + calculations[metric] = varianceDS + # + # sigma + # + metric = 'sigma' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + sigmaDS = self._computePower(0.5,calculations['variance'][list(needed[metric]['targets'])]) + self.calculations[metric] = sigmaDS + calculations[metric] = sigmaDS + # + # coeff of variation (sigma/mu) + # + metric = 'variationCoefficient' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + calculations[metric] = calculations['sigma'][needed[metric]['targets']] / calculations['expectedValue'][needed[metric]['targets']] + # + # skewness + # + metric = 'skewness' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + meanSet = calculations['expectedValue'][list(needed[metric]['targets'])] + varianceSet = calculations['variance'][list(needed[metric]['targets'])] + relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None + calculations[metric] = self._computeSkewness(dataSet,meanSet,varianceSet,pbWeight=relWeight,dim=self.sampleTag) + # + # kurtosis + # + metric = 'kurtosis' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + meanSet = calculations['expectedValue'][list(needed[metric]['targets'])] + varianceSet = calculations['variance'][list(needed[metric]['targets'])] + relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None + calculations[metric] = self._computeKurtosis(dataSet,meanSet,varianceSet,pbWeight=relWeight,dim=self.sampleTag) + # + # median + # + metric = 'median' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + if self.pbPresent: + # if all weights are the same, calculate percentile with xarray, no need for _computeWeightedPercentile + allSameWeight = True + for target in needed[metric]['targets']: + targWeight = relWeight[target].values + if targWeight.min() != targWeight.max(): + allSameWeight = False + if allSameWeight: + medianSet = dataSet.median(dim=self.sampleTag) + else: + medianSet = xr.Dataset() + relWeight = pbWeights[list(needed[metric]['targets'])] + for target in needed[metric]['targets']: + targWeight = relWeight[target].values + targDa = dataSet[target] + if self.pivotParameter in targDa.sizes.keys(): + quantile = [self._computeWeightedPercentile(group.values,targWeight,needed[metric]['interpolation'],percent=[0.5])[0] for label,group in targDa.groupby(self.pivotParameter)] + else: + quantile = self._computeWeightedPercentile(targDa.values,targWeight,needed[metric]['interpolation'],percent=[0.5])[0] + if self.pivotParameter in targDa.sizes.keys(): + da = xr.DataArray(quantile,dims=(self.pivotParameter),coords={self.pivotParameter:self.pivotValue}) + else: + da = xr.DataArray(quantile) + medianSet[target] = da + else: + medianSet = dataSet.median(dim=self.sampleTag) + self.calculations[metric] = medianSet + calculations[metric] = medianSet + # + # lowerPartialVariance + # + metric = 'lowerPartialVariance' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + medianSet = calculations['median'][list(needed[metric]['targets'])] + relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None + lowerPartialVarianceDS = self._computeLowerPartialVariance(dataSet,medianSet,pbWeight=relWeight,dim=self.sampleTag) + + calculations[metric] = lowerPartialVarianceDS + # + # lowerPartialSigma + # + metric = 'lowerPartialSigma' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + lpsDS = self._computePower(0.5,calculations['lowerPartialVariance'][list(needed[metric]['targets'])]) + calculations[metric] = lpsDS + # + # higherPartialVariance + # + metric = 'higherPartialVariance' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + medianSet = calculations['median'][list(needed[metric]['targets'])] + relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None + higherPartialVarianceDS = self._computeHigherPartialVariance(dataSet,medianSet,pbWeight=relWeight,dim=self.sampleTag) + + calculations[metric] = lowerPartialVarianceDS + # + # higherPartialSigma + # + metric = 'higherPartialSigma' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + hpsDS = self._computePower(0.5,calculations['higherPartialVariance'][list(needed[metric]['targets'])]) + calculations[metric] = hpsDS + + ############################################################ + # Begin Standard Error Calculations + # + # Reference for standard error calculations (including percentile): + # B. Harding, C. Tremblay and D. Cousineau, "Standard errors: A review and evaluation of + # standard error estimators using Monte Carlo simulations", The Quantitative Methods of + # Psychology, Vol. 10, No. 2 (2014) + ############################################################ + metric = 'expectedValue' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting calculate standard error on"'+metric+'"...') + if self.pbPresent: + factor = self._computePower(0.5,calculations['equivalentSamples']) + else: + factor = np.sqrt(self.sampleSize) + calculations[metric+'_ste'] = calculations['sigma'][list(needed[metric]['targets'])]/factor + + metric = 'variance' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting calculate standard error on "'+metric+'"...') + varList = list(needed[metric]['targets']) + if self.pbPresent: + en = calculations['equivalentSamples'][varList] + factor = 2.0 /(en - 1.0) + factor = self._computePower(0.5,factor) + else: + factor = np.sqrt(2.0/(float(self.sampleSize) - 1.0)) + calculations[metric+'_ste'] = calculations['sigma'][varList]**2 * factor + + metric = 'sigma' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting calculate standard error on "'+metric+'"...') + varList = list(needed[metric]['targets']) + if self.pbPresent: + en = calculations['equivalentSamples'][varList] + factor = 2.0 * (en - 1.0) + factor = self._computePower(0.5,factor) + else: + factor = np.sqrt(2.0 * (float(self.sampleSize) - 1.0)) + calculations[metric+'_ste'] = calculations['sigma'][varList] / factor + + metric = 'median' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting calculate standard error on "'+metric+'"...') + varList = list(needed[metric]['targets']) + calculations[metric+'_ste'] = calculations['expectedValue_ste'][varList] * np.sqrt(np.pi/2.0) + + metric = 'skewness' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting calculate standard error on "'+metric+'"...') + varList = list(needed[metric]['targets']) + if self.pbPresent: + en = calculations['equivalentSamples'][varList] + factor = 6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.)) + factor = self._computePower(0.5,factor) + calculations[metric+'_ste'] = xr.full_like(calculations[metric],1.0) * factor + else: + en = float(self.sampleSize) + factor = np.sqrt(6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.))) + calculations[metric+'_ste'] = xr.full_like(calculations[metric],factor) + + metric = 'kurtosis' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting calculate standard error on "'+metric+'"...') + varList = list(needed[metric]['targets']) + if self.pbPresent: + en = calculations['equivalentSamples'][varList] + factor1 = self._computePower(0.5,6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.))) + factor2 = self._computePower(0.5,(en**2-1.)/((en-3.0)*(en+5.0))) + factor = 2.0 * factor1 * factor2 + calculations[metric+'_ste'] = xr.full_like(calculations[metric],1.0) * factor + else: + en = float(self.sampleSize) + factor = 2.0 * np.sqrt(6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.)))*np.sqrt((en**2-1.)/((en-3.0)*(en+5.0))) + calculations[metric+'_ste'] = xr.full_like(calculations[metric],factor) + ############################################################ + # End of Standard Error Calculations + ############################################################ + # + # maximum + # + metric = 'maximum' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + calculations[metric] = dataSet.max(dim=self.sampleTag) + # + # minimum + # + metric = 'minimum' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + calculations[metric] = dataSet.min(dim=self.sampleTag) + # + # percentile, this metric is handled differently + # + metric = 'percentile' + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + dataSet = inputDataset[list(needed[metric]['targets'])] + percent = list(needed[metric]['percent']) + # are there probability weights associated with the data? + if self.pbPresent: + relWeight = pbWeights[list(needed[metric]['targets'])] + # if all weights are the same, calculate percentile with xarray, no need for _computeWeightedPercentile + allSameWeight = True + for target in needed[metric]['targets']: + targWeight = relWeight[target].values + if targWeight.min() != targWeight.max(): + allSameWeight = False + if allSameWeight: + # all weights are the same, percentile can be calculated with xarray.DataSet + percentileSet = dataSet.quantile(percent,dim=self.sampleTag,interpolation=needed[metric]['interpolation']) + percentileSet = percentileSet.rename({'quantile': 'percent'}) + else: + # probability weights are not all the same + # xarray does not have capability to calculate weighted quantiles at present + # implement our own solution + percentileSet = xr.Dataset() + for target in needed[metric]['targets']: + targWeight = relWeight[target].values + targDa = dataSet[target] + if self.pivotParameter in targDa.sizes.keys(): + quantile = [] + for label, group in targDa.groupby(self.pivotParameter): + qtl = self._computeWeightedPercentile(group.values, targWeight, needed[metric]['interpolation'], percent=percent) + quantile.append(qtl) + da = xr.DataArray(quantile, dims=(self.pivotParameter, 'percent'), coords={'percent': percent, self.pivotParameter: self.pivotValue}) + else: + quantile = self._computeWeightedPercentile(targDa.values, targWeight, needed[metric]['interpolation'], percent=percent) + da = xr.DataArray(quantile, dims=('percent'), coords={'percent': percent}) + + percentileSet[target] = da + + # TODO: remove when complete + # interpolation: {'linear', 'lower', 'higher','midpoint','nearest'}, do not try to use 'linear' or 'midpoint' + # The xarray.Dataset.where() will not return the corrrect solution + # 'lower' is used for consistent + # using xarray.Dataset.sel(**{'quantile':reqPercent}) to retrieve the quantile values + #dataSetWeighted = dataSet * relWeight + #percentileSet = dataSet.where(dataSetWeighted==dataSetWeighted.quantile(percent,dim=self.sampleTag,interpolation='lower')).mean(self.sampleTag) + else: + percentileSet = dataSet.quantile(percent,dim=self.sampleTag,interpolation=needed[metric]['interpolation']) + percentileSet = percentileSet.rename({'quantile':'percent'}) + calculations[metric] = percentileSet + + # because percentile is different, calculate standard error here + # standard error calculation uses the standard normal formulation for speed + self.raiseADebug('Starting calculate standard error on "'+metric+'"...') + norm = stats.norm + factor = np.sqrt(np.asarray(percent)*(1.0 - np.asarray(percent)))/norm.pdf(norm.ppf(percent)) + sigmaAdjusted = calculations['sigma'][list(needed[metric]['targets'])]/np.sqrt(calculations['equivalentSamples'][list(needed[metric]['targets'])]) + sigmaAdjusted = sigmaAdjusted.expand_dims(dim={'percent': percent}) + factor = xr.DataArray(data=factor, dims='percent', coords={'percent': percent}) + calculations[metric + '_ste'] = sigmaAdjusted*factor + + # # TODO: this is the KDE method, it is a more accurate method of calculating standard error + # # for percentile, but the computation time is too long. IF this computation can be sped up, + # # implement it here: + # percentileSteSet = xr.Dataset() + # calculatedPercentiles = calculations[metric] + # relWeight = pbWeights[list(needed[metric]['targets'])] + # for target in needed[metric]['targets']: + # targWeight = relWeight[target].values + # en = calculations['equivalentSamples'][target].values + # targDa = dataSet[target] + # if self.pivotParameter in targDa.sizes.keys(): + # percentileSte = np.zeros((len(self.pivotValue), len(percent))) # array + # for i, (label, group) in enumerate(targDa.groupby(self.pivotParameter)): # array + # if group.values.min() == group.values.max(): + # subPercentileSte = np.array([0.0]*len(percent)) + # else: + # # get KDE + # kde = stats.gaussian_kde(group.values, weights=targWeight) + # vals = calculatedPercentiles[target].sel(**{'percent': percent, self.pivotParameter: label}).values + # factor = np.sqrt(np.asarray(percent)*(1.0 - np.asarray(percent))/en) + # subPercentileSte = factor/kde(vals) + # percentileSte[i, :] = subPercentileSte + # da = xr.DataArray(percentileSte, dims=(self.pivotParameter, 'percent'), coords={self.pivotParameter: self.pivotValue, 'percent': percent}) + # percentileSteSet[target] = da + # else: + # calcPercentiles = calculatedPercentiles[target] + # if targDa.values.min() == targDa.values.max(): + # # distribution is a delta function, so no KDE construction + # percentileSte = list(np.zeros(calcPercentiles.shape)) + # else: + # # get KDE + # kde = stats.gaussian_kde(targDa.values, weights=targWeight) + # factor = np.sqrt(np.array(percent)*(1.0 - np.array(percent))/en) + # percentileSte = list(factor/kde(calcPercentiles.values)) + # da = xr.DataArray(percentileSte, dims=('percent'), coords={'percent': percent}) + # percentileSteSet[target] = da + # calculations[metric+'_ste'] = percentileSteSet + + def startVector(metric): + """ + Common method among all metrics for establishing parameters + @ In, metric, string, the name of the statistics metric to calculate + @ Out, targets, list(str), list of target parameter names (evaluate metrics for these) + @ Out, features, list(str), list of feature parameter names (evaluate with respect to these) + @ Out, skip, bool, if True it means either features or parameters were missing, so don't calculate anything + """ + # default to skipping, change that if we find criteria + targets = [] + features = [] + skip = True + if len(needed[metric]['targets'])>0: + self.raiseADebug('Starting "'+metric+'"...') + targets = list(needed[metric]['targets']) + features = list(needed[metric]['features']) + skip = False #True only if we don't have targets and features + if skip: + if metric not in self.skipped.keys(): + self.skipped[metric] = True + return targets,features,skip + + ################# + # VECTOR VALUES # + ################# + # + # sensitivity matrix + # + metric = 'sensitivity' + targets,features,skip = startVector(metric) + #NOTE sklearn expects the transpose of what we usually do in RAVEN, so #samples by #features + if not skip: + #for sensitivity matrix, we don't use numpy/scipy methods to calculate matrix operations, + #so we loop over targets and features + params = list(set(targets).union(set(features))) + dataSet = inputDataset[params] + relWeight = pbWeights[params] if self.pbPresent else None + intersectionSet = set(targets) & set(features) + if self.pivotParameter in dataSet.sizes.keys(): + dataSet = dataSet.to_array().transpose(self.pivotParameter,self.sampleTag,'variable') + featSet = dataSet.sel(**{'variable':features}).values + targSet = dataSet.sel(**{'variable':targets}).values + pivotVals = dataSet.coords[self.pivotParameter].values + da = None + for i in range(len(pivotVals)): + ds = self.sensitivityCalculation(features,targets,featSet[i,:,:],targSet[i,:,:],intersectionSet) + da = ds if da is None else xr.concat([da,ds], dim=self.pivotParameter) + da.coords[self.pivotParameter] = pivotVals + else: + # construct target and feature matrices + dataSet = dataSet.to_array().transpose(self.sampleTag,'variable') + featSet = dataSet.sel(**{'variable':features}).values + targSet = dataSet.sel(**{'variable':targets}).values + da = self.sensitivityCalculation(features,targets,featSet,targSet,intersectionSet) + calculations[metric] = da + # + # covariance matrix + # + metric = 'covariance' + targets,features,skip = startVector(metric) + if not skip: + # because the C implementation is much faster than picking out individual values, + # we do the full covariance matrix with all the targets and features. + # FIXME adding an alternative for users to choose pick OR do all, defaulting to something smart + # dependent on the percentage of the full matrix desired, would be better. + # IF this is fixed, make sure all the features and targets are requested for all the metrics + # dependent on this metric + params = list(set(targets).union(set(features))) + dataSet = inputDataset[params] + relWeight = pbWeights[params] if self.pbPresent else None + if self.pbPresent: + fact = (self.__computeUnbiasedCorrection(2, self.realizationWeight)).to_array().values if not self.biased else 1.0 + meanSet = (dataSet * relWeight).sum(dim = self.sampleTag) + else: + meanSet = dataSet.mean(dim = self.sampleTag) + fact = 1.0 / (float(dataSet.sizes[self.sampleTag]) - 1.0) if not self.biased else 1.0 / float(dataSet.sizes[self.sampleTag]) + targVars = list(dataSet.data_vars) + varianceSet = self._computeVariance(dataSet,meanSet,pbWeight=relWeight,dim=self.sampleTag) + dataSet = dataSet - meanSet + if self.pivotParameter in dataSet.sizes.keys(): + ds = None + paramDA = dataSet.to_array().transpose(self.pivotParameter,'variable',self.sampleTag).values + varianceDA = varianceSet[targVars].to_array().transpose(self.pivotParameter,'variable').values + pivotVals = dataSet.coords[self.pivotParameter].values + for i in range(len(pivotVals)): + # construct target and feature matrices + paramSamples = paramDA[i,...] + da = self.covarianceCalculation(paramDA[i,...],fact,varianceDA[i,:],targVars) + ds = da if ds is None else xr.concat([ds,da], dim=self.pivotParameter) + ds.coords[self.pivotParameter] = pivotVals + calculations[metric] = ds + else: + # construct target and feature matrices + paramSamples = dataSet.to_array().transpose('variable',self.sampleTag).values + varianceDA = varianceSet[targVars].to_array().values + da = self.covarianceCalculation(paramSamples,fact,varianceDA,targVars) + calculations[metric] = da + + def getCovarianceSubset(desired): + """ + @ In, desired, list(str), list of parameters to extract from covariance matrix + @ Out, reducedCov, xarray.DataArray, reduced covariance matrix + """ + if self.pivotParameter in desired: + self.raiseAnError(RuntimeError, 'The pivotParameter "{}" is among the parameters requested for performing statistics. Please remove!'.format(self.pivotParameter)) + reducedCov = calculations['covariance'].sel(**{'targets':desired,'features':desired}) + return reducedCov + # + # pearson matrix + # + # see comments in covariance for notes on C implementation + metric = 'pearson' + targets,features,skip = startVector(metric) + if not skip: + params = list(set(targets).union(set(features))) + reducedCovar = getCovarianceSubset(params) + targCoords = reducedCovar.coords['targets'].values + if self.pivotParameter in reducedCovar.sizes.keys(): + pivotCoords = reducedCovar.coords[self.pivotParameter].values + ds = None + for label, group in reducedCovar.groupby(self.pivotParameter): + corrMatrix = self.corrCoeff(group.values) + da = xr.DataArray(corrMatrix, dims=('targets','features'), coords={'targets':targCoords,'features':targCoords}) + ds = da if ds is None else xr.concat([ds,da], dim=self.pivotParameter) + ds.coords[self.pivotParameter] = pivotCoords + calculations[metric] = ds + else: + corrMatrix = self.corrCoeff(reducedCovar.values) + da = xr.DataArray(corrMatrix, dims=('targets','features'), coords={'targets':targCoords,'features':targCoords}) + calculations[metric] = da + # + # spearman matrix + # + # see RAVEN theory manual for a detailed explaination + # of the formulation used here + # + metric = 'spearman' + targets,features,skip = startVector(metric) + #NOTE sklearn expects the transpose of what we usually do in RAVEN, so #samples by #features + if not skip: + #for spearman matrix, we don't use numpy/scipy methods to calculate matrix operations, + #so we loop over targets and features + params = list(set(targets).union(set(features))) + dataSet = inputDataset[params] + relWeight = pbWeights[params] if self.pbPresent else None + if self.pivotParameter in dataSet.sizes.keys(): + dataSet = dataSet.to_array().transpose(self.pivotParameter,self.sampleTag,'variable') + featSet = dataSet.sel(**{'variable':features}).values + targSet = dataSet.sel(**{'variable':targets}).values + pivotVals = dataSet.coords[self.pivotParameter].values + da = None + for i in range(len(pivotVals)): + ds = self.spearmanCorrelation(features,targets,featSet[i,:,:],targSet[i,:,:],relWeight) + da = ds if da is None else xr.concat([da,ds], dim=self.pivotParameter) + da.coords[self.pivotParameter] = pivotVals + else: + # construct target and feature matrices + dataSet = dataSet.to_array().transpose(self.sampleTag,'variable') + featSet = dataSet.sel(**{'variable':features}).values + targSet = dataSet.sel(**{'variable':targets}).values + da = self.spearmanCorrelation(features,targets,featSet,targSet,relWeight) + calculations[metric] = da + # + # VarianceDependentSensitivity matrix + # The formula for this calculation is coming from: http://www.math.uah.edu/stat/expect/Matrices.html + # The best linear predictor: L(Y|X) = expectedValue(Y) + cov(Y,X) * [vc(X)]^(-1) * [X-expectedValue(X)] + # where Y is a vector of outputs, and X is a vector of inputs, cov(Y,X) is the covariance matrix of Y and X, + # vc(X) is the covariance matrix of X with itself. + # The variance dependent sensitivity matrix is defined as: cov(Y,X) * [vc(X)]^(-1) + metric = 'VarianceDependentSensitivity' + targets,features,skip = startVector(metric) + if not skip: + params = list(set(targets).union(set(features))) + reducedCovar = getCovarianceSubset(params) + targCoords = reducedCovar.coords['targets'].values + if self.pivotParameter in reducedCovar.sizes.keys(): + pivotCoords = reducedCovar.coords[self.pivotParameter].values + ds = None + for label, group in reducedCovar.groupby(self.pivotParameter): + da = self.varianceDepSenCalculation(targCoords,group.values) + ds = da if ds is None else xr.concat([ds,da], dim=self.pivotParameter) + ds.coords[self.pivotParameter] = pivotCoords + calculations[metric] = ds + else: + da = self.varianceDepSenCalculation(targCoords,reducedCovar.values) + calculations[metric] = da + + # + # Normalized variance dependent sensitivity matrix + # variance dependent sensitivity normalized by the mean (% change of output)/(% change of input) + # + metric = 'NormalizedSensitivity' + targets,features,skip = startVector(metric) + if not skip: + params = list(set(targets).union(set(features))) + reducedSen = calculations['VarianceDependentSensitivity'].sel(**{'targets':params,'features':params}) + meanDA = calculations['expectedValue'][params].to_array() + meanDA = meanDA.rename({'variable':'targets'}) + reducedSen /= meanDA + meanDA = meanDA.rename({'targets':'features'}) + reducedSen *= meanDA + calculations[metric] = reducedSen + + for metric, ds in calculations.items(): + if metric in self.scalarVals + self.steVals +['equivalentSamples'] and metric !='samples': + calculations[metric] = ds.to_array().rename({'variable':'targets'}) + # in here we fill the NaN with "nan". In this way, we are sure that even if + # there might be NaN in any raw for a certain timestep we do not drop the variable + # In the past, in a condition such as: + # time, A, B, C + # 0, 1, NaN, 1 + # 1, 1, 0.5, 1 + # 2, 1, 2.0, 2 + # the variable B would have been dropped (in the printing stage) + # with this modification, this should not happen anymore + outputSet = xr.Dataset(data_vars=calculations).fillna("nan") + + if self.outputDataset: + # Add 'RAVEN_sample_ID' to output dataset for consistence + if 'RAVEN_sample_ID' not in outputSet.sizes.keys(): + outputSet = outputSet.expand_dims('RAVEN_sample_ID') + outputSet['RAVEN_sample_ID'] = [0] + return outputSet + else: + outputDict = {} + for metric, requestList in self.toDo.items(): + for targetDict in requestList: + prefix = targetDict['prefix'].strip() + for target in targetDict['targets']: + if metric in self.scalarVals and metric != 'percentile': + varName = prefix + '_' + target + outputDict[varName] = np.atleast_1d(outputSet[metric].sel(**{'targets':target})) + steMetric = metric + '_ste' + if steMetric in self.steVals: + metaVar = prefix + '_ste_' + target + outputDict[metaVar] = np.atleast_1d(outputSet[steMetric].sel(**{'targets':target})) + elif metric == 'percentile': + for percent in targetDict['strPercent']: + varName = '_'.join([prefix,percent,target]) + percentVal = float(percent)/100. + outputDict[varName] = np.atleast_1d(outputSet[metric].sel(**{'targets':target,'percent':percentVal})) + steMetric = metric + '_ste' + if steMetric in self.steVals: + metaVar = '_'.join([prefix,percent,'ste',target]) + outputDict[metaVar] = np.atleast_1d(outputSet[steMetric].sel(**{'targets':target,'percent':percentVal})) + else: + #check if it was skipped for some reason + skip = self.skipped.get(metric, None) + if skip is not None: + self.raiseADebug('Metric',metric,'was skipped for parameters',targetDict,'! See warnings for details. Ignoring...') + continue + if metric in self.vectorVals: + for feature in targetDict['features']: + varName = '_'.join([prefix,target,feature]) + outputDict[varName] = np.atleast_1d(outputSet[metric].sel(**{'targets':target,'features':feature})) + if self.pivotParameter in outputSet.sizes.keys(): + outputDict[self.pivotParameter] = np.atleast_1d(self.pivotValue) + + return outputDict + + def corrCoeff(self, covM): + """ + This method calculates the correlation coefficient Matrix (pearson) for the given data. + Unbiased unweighted covariance matrix, weights is None, bias is 0 (default) + Biased unweighted covariance matrix, weights is None, bias is 1 + Unbiased weighted covariance matrix, weights is not None, bias is 0 + Biased weighted covariance matrix, weights is not None, bias is 1 + can be calcuated depending on the selection of the inputs. + @ In, covM, numpy.array, [#targets,#targets] covariance matrix + @ Out, covM, numpy.array, [#targets,#targets] correlation matrix + """ + try: + d = np.diag(covM) + except ValueError: + # scalar covariance + # nan if incorrect value (nan, inf, 0), 1 otherwise + return covM / covM + stdDev = np.sqrt(d) + covM /= stdDev[:,None] + covM /= stdDev[None,:] + return covM + + def sensitivityCalculation(self,featVars, targVars, featSamples, targSamples, intersectionSet): + """ + This method computes the sensitivity coefficients based on the SciKitLearn LinearRegression method + @ In, featVars, list, list of feature variables + @ In, targVars, list, list of target variables + @ In, featSamples, numpy.ndarray, [#samples, #features] array of features + @ In, targSamples, numpy.ndarray, [#samples, #targets] array of targets + @ In, intersectionSet, boolean, True if some target variables are in the list of features + @ Out, da, xarray.DataArray, contains the calculations of sensitivity coefficients + """ + from sklearn.linear_model import LinearRegression + if self.multipleFeatures: + # intersectionSet is flag that used to check the relationship between the features and targets. + # If True, part of the target variables are listed in teh feature set, then multivariate linear + # regression should not be used, and a loop over the target set is required. + # If False, which means there is no overlap between the target set and feature set. + # mutivariate linear regression can be used. However, for both cases, co-linearity check should be + # added for the feature set. ~ wangc + + if not intersectionSet: + condNumber = np.linalg.cond(featSamples) + if condNumber > 30.: + self.raiseAWarning("Condition Number: {:10.4f} > 30.0. Detected SEVERE multicollinearity problem. Sensitivity might be incorrect!".format(condNumber)) + senMatrix = LinearRegression().fit(featSamples,targSamples).coef_ + else: + # Target variables are in feature variables list, multi-target linear regression can not be used + # Since the 'multi-colinearity' exists, we need to loop over target variables + # TODO: Some general methods need to be implemented in order to handle the 'multi-colinearity' -- wangc + senMatrix = np.zeros((len(targVars), len(featVars))) + for p, targ in enumerate(targVars): + ind = list(featVars).index(targ) if targ in featVars else None + if ind is not None: + featMat = np.delete(featSamples,ind,axis=1) + else: + featMat = featSamples + regCoeff = LinearRegression().fit(featMat, targSamples[:,p]).coef_ + condNumber = np.linalg.cond(featMat) + if condNumber > 30.: + self.raiseAWarning("Condition Number: {:10.4f} > 30.0. Detected SEVERE multicollinearity problem. Sensitivity might be incorrect!".format(condNumber)) + if ind is not None: + regCoeff = np.insert(regCoeff,ind,1.0) + senMatrix[p,:] = regCoeff + else: + senMatrix = np.zeros((len(targVars), len(featVars))) + for p, feat in enumerate(featVars): + regCoeff = LinearRegression().fit(featSamples[:,p].reshape(-1,1),targSamples).coef_ + senMatrix[:,p] = regCoeff[:,0] + da = xr.DataArray(senMatrix, dims=('targets','features'), coords={'targets':targVars,'features':featVars}) + + return da + + def covarianceCalculation(self,paramSamples,fact,variance,targVars): + """ + This method computes the covariance of given sample matrix + @ In, paramSamples, numpy.ndarray, [#parameters, #samples], array of parameters + @ In, fact, float, the unbiase correction factor + @ In, variance, numpy.ndarray, [#parameters], variance of parameters + @ In, targVars, list, the list of parameters + @ Out, da, xarray.DataArray, contains the calculations of covariance + """ + if self.pbPresent: + paramSamplesT = (paramSamples*self.realizationWeight['ProbabilityWeight'].values).T + else: + paramSamplesT = paramSamples.T + cov = np.dot(paramSamples, paramSamplesT.conj()) + cov *= fact + np.fill_diagonal(cov,variance) + da = xr.DataArray(cov, dims=('targets','features'), coords={'targets':targVars,'features':targVars}) + return da + + def varianceDepSenCalculation(self,targCoords, cov): + """ + This method computes the covariance of given sample matrix + @ In, targCoords, list, the list of parameters + @ In, cov, numpy.ndarray, the covariance of parameters + @ Out, da, xarray.DataArray, contains the calculations of variance dependent sensitivities + """ + senMatrix = np.zeros((len(targCoords), len(targCoords))) + if self.multipleFeatures: + for p, param in enumerate(targCoords): + covX = np.delete(cov,p,axis=0) + covX = np.delete(covX,p,axis=1) + covYX = np.delete(cov[p,:],p) + sensCoef = np.dot(covYX,np.linalg.pinv(covX)) + sensCoef = np.insert(sensCoef,p,1.0) + senMatrix[p,:] = sensCoef + else: + for p, param in enumerate(targCoords): + covX = cov[p,p] + covYX = cov[:,p] + sensCoef = covYX / covX + senMatrix[:,p] = sensCoef + da = xr.DataArray(senMatrix, dims=('targets','features'), coords={'targets':targCoords,'features':targCoords}) + return da + + def spearmanCorrelation(self, featVars, targVars, featSamples, targSamples, pbWeights): + """ + This method computes the spearman correlation coefficients + @ In, featVars, list, list of feature variables + @ In, targVars, list, list of target variables + @ In, featSamples, numpy.ndarray, [#samples, #features] array of features + @ In, targSamples, numpy.ndarray, [#samples, #targets] array of targets + @ In, pbWeights, dataset, probability weights + @ Out, da, xarray.DataArray, contains the calculations of spearman coefficients + """ + spearmanMat = np.zeros((len(targVars), len(featVars))) + wf, wt = None, None + # compute unbiased factor + if self.pbPresent: + fact = (self.__computeUnbiasedCorrection(2, self.realizationWeight)).to_array().values if not self.biased else 1.0 + vp = self.__computeVp(1,self.realizationWeight)['ProbabilityWeight'].values + varianceFactor = fact*(1.0/vp) + else: + fact = 1.0 / (float(featSamples.shape[0]) - 1.0) if not self.biased else 1.0 / float(featSamples.shape[0]) + varianceFactor = fact + + for tidx, target in enumerate(targVars): + for fidx, feat in enumerate(featVars): + if self.pbPresent: + wf, wt = np.asarray(pbWeights[feat]), np.asarray(pbWeights[target]) + rankFeature, rankTarget = mathUtils.rankData(featSamples[:,fidx],wf), mathUtils.rankData(targSamples[:,tidx],wt) + # compute covariance of the ranked features + cov = np.cov(rankFeature, y=rankTarget, aweights=wt) + covF = np.cov(rankFeature,y=rankFeature, aweights=wf) + covT = np.cov(rankTarget,y=rankTarget, aweights=wf) + # apply correction factor (for biased or unbiased) (off diagonal) + cov[~np.eye(2,dtype=bool)] *= fact + covF[~np.eye(2,dtype=bool)] *= fact + covT[~np.eye(2,dtype=bool)] *= fact + # apply correction factor (for biased or unbiased) (diagonal) + cov[np.eye(2,dtype=bool)] *= varianceFactor + covF[~np.eye(2,dtype=bool)] *= varianceFactor + covT[~np.eye(2,dtype=bool)] *= varianceFactor + # now we can compute the pearson of such pairs + spearman = (cov / np.sqrt(covF * covT))[-1,0] + spearmanMat[tidx,fidx] = spearman + + da = xr.DataArray(spearmanMat, dims=('targets','features'), coords={'targets':targVars,'features':featVars}) + return da + + def run(self, inputIn): + """ + This method executes the postprocessor action. In this case, it computes all the requested statistical FOMs + @ In, inputIn, object, object contained the data to process. (inputToInternal output) + @ Out, outputSet, xarray.Dataset or dictionary, dataset or dictionary containing the results + """ + inputData = self.inputToInternal(inputIn) + outputSet = self.__runLocal(inputData) + return outputSet + + def collectOutput(self, finishedJob, output): + """ + Function to place all of the computed data into the output object + @ In, finishedJob, JobHandler External or Internal instance, A JobHandler object that is in charge of running this post-processor + @ In, output, dataObjects, The object where we want to place our computed results + @ Out, None + """ + super().collectOutput(finishedJob, output) From e225c0ca18ba0e59c84962735918f14332d99b36 Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Fri, 12 May 2023 15:05:12 -0600 Subject: [PATCH 02/10] updated grid entity --- ravenframework/GridEntities.py | 80 +- .../Models/PostProcessors/BasicStatistics.py | 1 - .../Models/PostProcessors/Factory.py | 1 + .../SubdomainBasicStatistics.py | 1432 +---------------- ravenframework/Samplers/Grid.py | 3 +- ravenframework/utils/mathUtils.py | 1 - 6 files changed, 125 insertions(+), 1393 deletions(-) diff --git a/ravenframework/GridEntities.py b/ravenframework/GridEntities.py index c4968e52d3..6e188b3cce 100644 --- a/ravenframework/GridEntities.py +++ b/ravenframework/GridEntities.py @@ -58,15 +58,16 @@ def __init__(self): self.gridInitDict = {} @classmethod - def _readMoreXml(cls, xmlNode, dimensionTags=None, dimTagsPrefix=None): + def _handleInput(self, paramInput, dimensionTags=None, dimTagsPrefix=None): """ - XML reader for the grid statement. - @ In, xmlNode, xml.etree.ElementTree.Element, XML element node that represents the portion of the input that belongs to this class + Function to handle the parsed paramInput for this class. + @ In, paramInput, ParameterInput, the already parsed input. @ In, dimensionTag, list, optional, names of the tag that represents the grid dimensions @ In, dimTagsPrefix, dict, optional, eventual prefix to use for defining the dimName @ Out, None """ + @classmethod def initialize(cls, initDictionary=None): """ @@ -252,16 +253,16 @@ def __init__(self): self.constructTensor = False # True if we need to construct the tensor product of the the ND grid (full grid) or just the iterator (False) self.uniqueCellNumber = 0 # number of unique cells self.gridIterator = None # the grid iterator - self.gridInitDict = {} # dictionary with initialization grid info from _readMoreXML. If None, the "initialize" method will look for all the information in the in Dictionary + self.gridInitDict = {} # dictionary with initialization grid info from _handleInput. If None, the "initialize" method will look for all the information in the in Dictionary self.volumetricRatio = None # volumetric ratio (optional if steplenght is read or passed in initDict) self.nVar = None self.dimName = None self.gridInitDict = {} - def _readMoreXml(self, xmlNode, dimensionTags=None, dimTagsPrefix=None): + def _handleInput(self, paramInput, dimensionTags=None, dimTagsPrefix=None): """ - XML reader for the grid statement. - @ In, xmlNode, xml.etree.ElementTree.Element, XML element node that represents the portion of the input that belongs to this class + Function to handle the parsed paramInput for the grid (single) instance. + @ In, paramInput, ParameterInput, the already parsed input. @ In, dimensionTag, list, optional, names of the tag that represents the grid dimensions @ In, dimTagsPrefix, dict, optional, eventual prefix to use for defining the dimName @ Out, None @@ -269,22 +270,22 @@ def _readMoreXml(self, xmlNode, dimensionTags=None, dimTagsPrefix=None): self.gridInitDict = {'dimensionNames':[],'lowerBounds':{},'upperBounds':{},'stepLength':{}} gridInfo = {} dimInfo = {} - for child in xmlNode: + for child in paramInput.subparts: self.dimName = None if dimensionTags is not None: - if child.tag in dimensionTags: - self.dimName = child.attrib['name'] + if child.getName() in dimensionTags: + self.dimName = child.parameterValues['name'] if dimTagsPrefix is not None: - self.dimName = dimTagsPrefix[child.tag] + self.dimName if child.tag in dimTagsPrefix else self.dimName - if child.tag == "grid": - gridInfo[self.dimName] = self._readGridStructure(child, xmlNode) - for childChild in child: - if childChild.tag == "grid": + self.dimName = dimTagsPrefix[child.getName()] + self.dimName if child.getName() in dimTagsPrefix else self.dimName + if child.getName() == "grid": + gridInfo[self.dimName] = self._readGridStructure(child, paramInput) + for childChild in child.subparts: + if childChild.getName() == "grid": gridInfo[self.dimName] = self._readGridStructure(childChild,child) - if 'dim' in childChild.attrib: + if 'dim' in childChild.parameterValues: dimID = str(len(self.gridInitDict['dimensionNames']) + 1) if self.dimName is None else self.dimName try: - dimInfo[dimID] = [int(childChild.attrib['dim']),None] + dimInfo[dimID] = [int(childChild.parameterValues['dim']),None] except ValueError: self.raiseAnError(ValueError, "cannot convert 'dim' attribute in integer!") # check for globalGrid type of structure @@ -309,51 +310,51 @@ def _readMoreXml(self, xmlNode, dimensionTags=None, dimTagsPrefix=None): def _readGridStructure(self, child, parent): """ This method is aimed to read the grid structure in the xml node - @ In, child, xml.etree.ElementTree.Element, the xml node containing the grid info - @ In, parent, xml.etree.ElementTree.Element, the xml node that contains the node in which the grid info are defined + @ In, child, ParameterInput, the already parsed input containing the grid info + @ In, parent, ParameterInput, the already parsed input containing the node in which the grid info are defined @ Out, gridStruct, tuple, the grid structure read ((type, construction type, upper and lower bounds), gridName) """ - if child.tag =='grid': + if child.getName() =='grid': gridStruct, gridName = self._fillGrid(child) if self.dimName is None: self.dimName = str(len(self.gridInitDict['dimensionNames'])+1) - if parent.tag != 'globalGrid': + if parent.getName() != 'globalGrid': self.gridInitDict['dimensionNames'].append(self.dimName) else: if gridName is None: self.raiseAnError(IOError, 'grid defined in globalGrid block must have the attribute "name"!') - self.dimName = parent.tag + ':' + gridName + self.dimName = parent.getName() + ':' + gridName return gridStruct def _fillGrid(self, child): """ This method is aimed to fill the grid structure from an XML node - @ In, child, xml.etree.ElementTree.Element, the xml node containing the grid info + @ In, child, ParameterInput, the already parsed input containing the grid info @ Out, gridStruct, tuple, the grid structure read ((type, construction type, upper and lower bounds), gridName) """ constrType = None - if 'construction' in child.attrib.keys(): - constrType = child.attrib['construction'] - if 'type' not in child.attrib.keys(): + if 'construction' in child.parameterValues.keys(): + constrType = child.parameterValues['construction'] + if 'type' not in child.parameterValues.keys(): self.raiseAnError(IOError, "Each XML node needs to have the attribute type!!!!") nameGrid = None if constrType in ['custom','equal']: - bounds = [floatConversion(element) for element in child.text.split()] + bounds = [floatConversion(element) for element in child.value.split()] bounds.sort() lower, upper = min(bounds), max(bounds) - if 'name' in child.attrib.keys(): - nameGrid = child.attrib['name'] + if 'name' in child.parameterValues.keys(): + nameGrid = child.parameterValues['name'] if constrType == 'custom': - gridStruct = (child.attrib['type'],constrType,bounds),nameGrid + gridStruct = (child.parameterValues['type'],constrType,bounds),nameGrid elif constrType == 'equal': if len(bounds) != 2: - self.raiseAnError(IOError, f'body of grid XML node needs to contain 2 values (lower and upper bounds).Tag = {child.tag}') - if 'steps' not in child.attrib.keys(): + self.raiseAnError(IOError, f'body of grid XML node needs to contain 2 values (lower and upper bounds).Tag = {child.getName()}') + if 'steps' not in child.parameterValues.keys(): self.raiseAnError(IOError, 'the attribute step needs to be inputted when "construction" attribute == equal!') - gridStruct = (child.attrib['type'],constrType,np.linspace(lower,upper,partialEval(child.attrib['steps'])+1)),nameGrid - elif child.attrib['type'] == 'globalGrid': - gridStruct = (child.attrib['type'],constrType,child.text),nameGrid + gridStruct = (child.parameterValues['type'],constrType,np.linspace(lower,upper,partialEval(child.parameterValues['steps'])+1)),nameGrid + elif child.parameterValues['type'] == 'globalGrid': + gridStruct = (child.parameterValues['type'],constrType,child.value),nameGrid else: self.raiseAnError(IOError, f'construction type unknown! Got: {constrType}') @@ -759,15 +760,16 @@ def __len__(self): return totalLength - def _readMoreXml(self, xmlNode, dimensionTags=None, dimTagsPrefix=None): + + def _handleInput(self, paramInput, dimensionTags=None, dimTagsPrefix=None): """ - XML reader for the Multi-grid statement. - @ In, xmlNode, xml.etree.ElementTree, XML element node that represents the portion of the input that belongs to this class + Function to handle the parsed paramInput for the grid (multi-grid) instance. + @ In, paramInput, ParameterInput, the already parsed input. @ In, dimensionTag, list, optional, names of the tag that represents the grid dimensions @ In, dimTagsPrefix, dict, optional, eventual prefix to use for defining the dimName @ Out, None """ - self.grid.getrootnode().get("grid")._readMoreXml(xmlNode, dimensionTags, dimTagsPrefix) + self.grid.getrootnode().get("grid")._handleInput(paramInput, dimensionTags, dimTagsPrefix) def initialize(self, initDictionary=None): """ diff --git a/ravenframework/Models/PostProcessors/BasicStatistics.py b/ravenframework/Models/PostProcessors/BasicStatistics.py index 3a54d20c48..34450cca98 100644 --- a/ravenframework/Models/PostProcessors/BasicStatistics.py +++ b/ravenframework/Models/PostProcessors/BasicStatistics.py @@ -31,7 +31,6 @@ from ...utils import utils from ...utils import InputData, InputTypes from ...utils import mathUtils -from ... import Files #Internal Modules End----------------------------------------------------------- class BasicStatistics(PostProcessorInterface): diff --git a/ravenframework/Models/PostProcessors/Factory.py b/ravenframework/Models/PostProcessors/Factory.py index b74904937e..683cfe0d2f 100644 --- a/ravenframework/Models/PostProcessors/Factory.py +++ b/ravenframework/Models/PostProcessors/Factory.py @@ -19,6 +19,7 @@ from .PostProcessorInterface import PostProcessorInterface from .PostProcessorReadyInterface import PostProcessorReadyInterface from .BasicStatistics import BasicStatistics +from .SubdomainBasicStatistics import SubdomainBasicStatistics from .LimitSurface import LimitSurface from .Metric import Metric from .DataMining import DataMining diff --git a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py index b3b5e1d40f..47dd663676 100644 --- a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py +++ b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -Created on July 10, 2013 +Created on May 10, 2023 @author: aalfonsi """ @@ -32,45 +32,13 @@ from ...utils import utils from ...utils import InputData, InputTypes from ...utils import mathUtils -from ... import Files #Internal Modules End----------------------------------------------------------- class SubdomainBasicStatistics(PostProcessorInterface): """ - BasicStatistics filter class. It computes all the most popular statistics + Subdomain basic statitistics class. It computes all statistics on subdomains """ - scalarVals = ['expectedValue', - 'minimum', - 'maximum', - 'median', - 'variance', - 'sigma', - 'percentile', - 'variationCoefficient', - 'skewness', - 'kurtosis', - 'samples', - 'higherPartialVariance', # Statistic metric not available yet - 'higherPartialSigma', # Statistic metric not available yet - 'lowerPartialSigma', # Statistic metric not available yet - 'lowerPartialVariance' # Statistic metric not available yet - ] - vectorVals = ['sensitivity', - 'covariance', - 'pearson', - 'spearman', - 'NormalizedSensitivity', - 'VarianceDependentSensitivity'] - # quantities that the standard error can be computed - steVals = ['expectedValue_ste', - 'median_ste', - 'variance_ste', - 'sigma_ste', - 'skewness_ste', - 'kurtosis_ste', - 'percentile_ste'] - @classmethod def getInputSpecification(cls): """ @@ -80,65 +48,24 @@ class cls. @ Out, inputSpecification, InputData.ParameterInput, class to use for specifying input of cls. """ - ## This will replace the lines above - inputSpecification = super().getInputSpecification() - - for scalar in cls.scalarVals: - scalarSpecification = InputData.parameterInputFactory(scalar, contentType=InputTypes.StringListType) - if scalar == 'percentile': - #percent is a string type because otherwise we can't tell 95.0 from 95 - # which matters because the number is used in output. - scalarSpecification.addParam("percent", InputTypes.StringListType) - # percentile has additional "interpolation" parameter - scalarSpecification.addParam("interpolation", - param_type=InputTypes.makeEnumType("interpolation", - "interpolationType", - ["linear", "midpoint"]), - default="linear", - descr="""Interpolation method for percentile calculation. - 'linear' uses linear interpolation between nearest - data points while 'midpoint' uses the average of the - nearest data points.""") - if scalar == 'median': - # median has additional "interpolation" parameter - scalarSpecification.addParam("interpolation", - param_type=InputTypes.makeEnumType("interpolation", - "interpolationType", - ["linear", "midpoint"]), - default="linear", - descr="""Interpolation method for median calculation. 'linear' - uses linear interpolation between nearest data points - while 'midpoint' uses the average of the nearest data - points.""") - scalarSpecification.addParam("prefix", InputTypes.StringType) - inputSpecification.addSub(scalarSpecification) - - for vector in cls.vectorVals: - vectorSpecification = InputData.parameterInputFactory(vector) - vectorSpecification.addParam("prefix", InputTypes.StringType) - features = InputData.parameterInputFactory('features', - contentType=InputTypes.StringListType) - vectorSpecification.addSub(features) - targets = InputData.parameterInputFactory('targets', - contentType=InputTypes.StringListType) - vectorSpecification.addSub(targets) - inputSpecification.addSub(vectorSpecification) - - pivotParameterInput = InputData.parameterInputFactory('pivotParameter', contentType=InputTypes.StringType) - inputSpecification.addSub(pivotParameterInput) - - datasetInput = InputData.parameterInputFactory('dataset', contentType=InputTypes.BoolType) - inputSpecification.addSub(datasetInput) - - methodsToRunInput = InputData.parameterInputFactory("methodsToRun", contentType=InputTypes.StringType) - inputSpecification.addSub(methodsToRunInput) - - biasedInput = InputData.parameterInputFactory("biased", contentType=InputTypes.BoolType) - inputSpecification.addSub(biasedInput) - - multipleFeaturesInput = InputData.parameterInputFactory("multipleFeatures", contentType=InputTypes.BoolType) - inputSpecification.addSub(multipleFeaturesInput) - + # We get the input specs from the basic statistics and we just add + # the subdomain info + inputSpecification = BasicStatistics.getInputSpecification() + + subdomainInputs = InputData.parameterInputFactory("subdomain", printPriority=100, + descr='defines the subdomain specs to be used for the subdomain statistics') + variableInput = InputData.parameterInputFactory("variable", printPriority=80, + descr="defines the variables to be used for the subdomain statistics.") + variableInput.addParam("name", InputTypes.StringNoLeadingSpacesType, + descr=r"""Name of the variable for this grid/subdomain. \nb As for the other objects, + this is the name that can be used to refer to this specific entity from other input blocks""") + gridInput = InputData.parameterInputFactory("grid", contentType=InputTypes.StringType) + gridInput.addParam("type", InputTypes.StringType) + gridInput.addParam("construction", InputTypes.StringType) + gridInput.addParam("steps", InputTypes.IntegerType) + variableInput.addSub(gridInput) + subdomainInputs.addSub(variableInput) + inputSpecification.addSub(subdomainInputs) return inputSpecification def __init__(self): @@ -148,23 +75,10 @@ def __init__(self): @ Out, None """ super().__init__() - self.parameters = {} # parameters dictionary (they are basically stored into a dictionary identified by tag "targets" - self.acceptedCalcParam = self.scalarVals + self.vectorVals - self.what = self.acceptedCalcParam # what needs to be computed... default...all - self.methodsToRun = [] # if a function is present, its outcome name is here stored... if it matches one of the known outcomes, the pp is going to use the function to compute it - self.printTag = 'PostProcessor BASIC STATISTIC' - self.biased = False # biased statistics? - self.pivotParameter = None # time-dependent statistics pivot parameter - self.pivotValue = None # time-dependent statistics pivot parameter values - self.dynamic = False # is it time-dependent? - self.sampleTag = None # Tag used to track samples - self.pbPresent = False # True if the ProbabilityWeight is available - self.realizationWeight = None # The joint probabilities - self.steMetaIndex = 'targets' # when Dataset is requested as output, the default index of ste metadata is ['targets', self.pivotParameter] - self.multipleFeatures = True # True if multiple features are employed in linear regression as feature inputs - self.sampleSize = None # number of sample size - self.calculations = {} - self.validDataType = ['PointSet', 'HistorySet', 'DataSet'] # The list of accepted types of DataObject + from ...Models.PostProcessors import factory as ppFactory # delay import to allow definition + self.stat = ppFactory.returnInstance('BasicStatistics') + self.validDataType = ['PointSet', 'HistorySet', 'DataSet'] + self.printTag = 'PostProcessor SUBDOMAIN STATISTICS' def inputToInternal(self, currentInp): """ @@ -268,72 +182,8 @@ def initialize(self, runInfo, inputs, initDict): @ In, initDict, dict, dictionary with initialization options @ Out, None """ - #construct a list of all the parameters that have requested values into self.allUsedParams - self.allUsedParams = set() - for metricName in self.scalarVals + self.vectorVals: - if metricName in self.toDo.keys(): - for entry in self.toDo[metricName]: - self.allUsedParams.update(entry['targets']) - try: - self.allUsedParams.update(entry['features']) - except KeyError: - pass + self.stat.intialize(runInfo, inputs, initDict) - #for backward compatibility, compile the full list of parameters used in Basic Statistics calculations - self.parameters['targets'] = list(self.allUsedParams) - super().initialize(runInfo, inputs, initDict) - inputObj = inputs[-1] if type(inputs) == list else inputs - if inputObj.type == 'HistorySet': - self.dynamic = True - inputMetaKeys = [] - outputMetaKeys = [] - for metric, infos in self.toDo.items(): - if metric in self.scalarVals + self.vectorVals: - steMetric = metric + '_ste' - if steMetric in self.steVals: - for info in infos: - prefix = info['prefix'] - for target in info['targets']: - if metric == 'percentile': - for strPercent in info['strPercent']: - metaVar = prefix + '_' + strPercent + '_ste_' + target if not self.outputDataset else metric + '_ste' - metaDim = inputObj.getDimensions(target) - if len(metaDim[target]) == 0: - inputMetaKeys.append(metaVar) - else: - outputMetaKeys.append(metaVar) - else: - metaVar = prefix + '_ste_' + target if not self.outputDataset else metric + '_ste' - metaDim = inputObj.getDimensions(target) - if len(metaDim[target]) == 0: - inputMetaKeys.append(metaVar) - else: - outputMetaKeys.append(metaVar) - metaParams = {} - if not self.outputDataset: - if len(outputMetaKeys) > 0: - metaParams = {key:[self.pivotParameter] for key in outputMetaKeys} - else: - if len(outputMetaKeys) > 0: - params = {} - for key in outputMetaKeys + inputMetaKeys: - # percentile standard error has additional index - if key == 'percentile_ste': - params[key] = [self.pivotParameter, self.steMetaIndex, 'percent'] - else: - params[key] = [self.pivotParameter, self.steMetaIndex] - metaParams.update(params) - elif len(inputMetaKeys) > 0: - params = {} - for key in inputMetaKeys: - # percentile standard error has additional index - if key == 'percentile_ste': - params[key] = [self.steMetaIndex, 'percent'] - else: - params[key] = [self.steMetaIndex] - metaParams.update(params) - metaKeys = inputMetaKeys + outputMetaKeys - self.addMetaKeys(metaKeys,metaParams) def _handleInput(self, paramInput, childVals=None): """ @@ -342,338 +192,72 @@ def _handleInput(self, paramInput, childVals=None): @ In, childVals, list, optional, quantities requested from child statistical object @ Out, None """ - if childVals is None: - childVals = [] - self.toDo = {} - for child in paramInput.subparts: - tag = child.getName() - #because percentile is strange (has attached parameters), we address it first - if tag in self.scalarVals + self.vectorVals: - if 'prefix' not in child.parameterValues: - self.raiseAnError(IOError, "No prefix is provided for node: ", tag) - #get the prefix - prefix = child.parameterValues['prefix'] - if tag == 'percentile': - #get targets - targets = set(child.value) - #what if user didn't give any targets? - if len(targets)<1: - self.raiseAWarning('No targets were specified in text of <'+tag+'>! Skipping metric...') - continue - #prepare storage dictionary, keys are percentiles, values are set(targets) - if tag not in self.toDo.keys(): - self.toDo[tag] = [] # list of {'targets':(), 'prefix':str, 'percent':str, 'interpolation':str} - if 'percent' not in child.parameterValues: - reqPercent = [0.05, 0.95] - strPercent = ['5','95'] - else: - reqPercent = set(utils.floatConversion(percent)/100. for percent in child.parameterValues['percent']) - strPercent = set(percent for percent in child.parameterValues['percent']) - if 'interpolation' not in child.parameterValues: - interpolation = 'linear' - else: - interpolation = child.parameterValues['interpolation'] - self.toDo[tag].append({'targets':set(targets), - 'prefix':prefix, - 'percent':reqPercent, - 'strPercent':strPercent, - 'interpolation':interpolation}) - # median also has an attached parameter - elif tag == 'median': - if tag not in self.toDo.keys(): - self.toDo[tag] = [] # list of {'targets':{}, 'prefix':str, 'interpolation':str} - if 'interpolation' not in child.parameterValues: - interpolation = 'linear' - else: - interpolation = child.parameterValues['interpolation'] - self.toDo[tag].append({'targets':set(child.value), - 'prefix':prefix, - 'interpolation':interpolation}) - elif tag in self.scalarVals: - if tag not in self.toDo.keys(): - self.toDo[tag] = [] # list of {'targets':(), 'prefix':str} - self.toDo[tag].append({'targets':set(child.value), - 'prefix':prefix}) - elif tag in self.vectorVals: - if tag not in self.toDo.keys(): - self.toDo[tag] = [] # list of {'targets':(),'features':(), 'prefix':str} - tnode = child.findFirst('targets') - if tnode is None: - self.raiseAnError('Request for vector value <'+tag+'> requires a "targets" node, and none was found!') - fnode = child.findFirst('features') - if fnode is None: - self.raiseAnError('Request for vector value <'+tag+'> requires a "features" node, and none was found!') - # we're storing toDo[tag] as a list of dictionaries. This is because the user might specify multiple - # nodes with the same metric (tag), but with different targets and features. For instance, the user might - # want the sensitivity of A and B to X and Y, and the sensitivity of C to W and Z, but not the sensitivity - # of A to W. If we didn't keep them separate, we could potentially waste a fair number of calculations. - self.toDo[tag].append({'targets':set(tnode.value), - 'features':set(fnode.value), - 'prefix':prefix}) - elif tag == "biased": - self.biased = child.value - elif tag == "pivotParameter": - self.pivotParameter = child.value - elif tag == "dataset": - self.outputDataset = child.value - elif tag == "multipleFeatures": - self.multipleFeatures = child.value - else: - if tag not in childVals: - self.raiseAWarning('Unrecognized node in BasicStatistics "',tag,'" has been ignored!') - - assert (len(self.toDo)>0), self.raiseAnError(IOError, 'BasicStatistics needs parameters to work on! Please check input for PP: ' + self.name) - - def _computePower(self, p, dataset): - """ - Compute the p-th power of weights - @ In, p, int, the power - @ In, dataset, xarray.Dataset, probability weights of all input variables - @ Out, pw, xarray.Dataset, the p-th power of weights - """ - pw = {} - coords = dataset.coords - for target, targValue in dataset.variables.items(): - ##remove index variable - if target in coords: - continue - pw[target] = np.power(targValue,p) - pw = xr.Dataset(data_vars=pw,coords=coords) - return pw - - def __computeVp(self,p,weights): - """ - Compute the sum of p-th power of weights - @ In, p, int, the power - @ In, weights, xarray.Dataset, probability weights of all input variables - @ Out, vp, xarray.Dataset, the sum of p-th power of weights - """ - vp = self._computePower(p,weights) - vp = vp.sum() - return vp - - def __computeEquivalentSampleSize(self,weights): - """ - Compute the equivalent sample size for given probability weights - @ In, weights, xarray.Dataset, probability weights of all input variables - @ Out, equivalentSize, xarray.Dataset, the equivalent sample size - """ - # The equivalent sample size for given samples, i.e. (sum of weights) squared / sum of the squared weights - # The definition of this quantity can be found: - # R. F. Potthoff, M. A. Woodbury and K. G. Manton, "'Equivalent SampleSize' and 'Equivalent Degrees of Freedom' - # Refinements for Inference Using Survey Weights Under Superpopulation Models", Journal of the American Statistical - # Association, Vol. 87, No. 418 (1992) - v1Square = self.__computeVp(1,weights)**2 - v2 = self.__computeVp(2,weights) - equivalentSize = v1Square/v2 - return equivalentSize + # initialize basic stats + subdomain = paramInput.popSub('subdomain') + if subdomain is None: + self.raiseAnError(IOError,' tag not found!') + self.stat._handleInput(paramInput, childVals) + + for child in subdomain.subparts: + if child.tag == 'variable': + varName = child.parameterValues['name'] + for childChild in child.subparts: + + + + + # variable for tracking if distributions or functions have been declared + foundDistOrFunc = False + # store variable name for re-use + varName = child.parameterValues['name'] + # set shape if present + if 'shape' in child.parameterValues: + self.variableShapes[varName] = child.parameterValues['shape'] + # read subnodes + for childChild in child.subparts: + if childChild.getName() == 'distribution': + # can only have a distribution if doesn't already have a distribution or function + if foundDistOrFunc: + self.raiseAnError(IOError, 'A sampled variable cannot have both a distribution and a function, or more than one of either!') + else: + foundDistOrFunc = True + # name of the distribution to sample + toBeSampled = childChild.value + varData = {} + varData['name'] = childChild.value + # variable dimensionality + if 'dim' not in childChild.parameterValues: + dim = 1 + else: + dim = childChild.parameterValues['dim'] + varData['dim'] = dim + # set up mapping for variable to distribution + self.variables2distributionsMapping[varName] = varData + # flag distribution as needing to be sampled + self.toBeSampled[prefix + varName] = toBeSampled + elif childChild.getName() == 'function': + # can only have a function if doesn't already have a distribution or function + if not foundDistOrFunc: + foundDistOrFunc = True + else: + self.raiseAnError(IOError, 'A sampled variable cannot have both a distribution and a function!') + # function name + toBeSampled = childChild.value + # track variable as a functional sample + self.dependentSample[prefix + varName] = toBeSampled - def __computeUnbiasedCorrection(self,order,weightsOrN): - """ - Compute unbiased correction given weights and momement order - Reference paper: - Lorenzo Rimoldini, "Weighted skewness and kurtosis unbiased by sample size", http://arxiv.org/pdf/1304.6564.pdf - @ In, order, int, moment order - @ In, weightsOrN, xarray.Dataset or int, if xarray.Dataset -> weights else -> number of samples - @ Out, corrFactor, xarray.Dataset or int, xarray.Dataset (order <=3) or tuple of xarray.Dataset (order ==4), - the unbiased correction factor if weightsOrN is xarray.Dataset else integer - """ - if order > 4: - self.raiseAnError(RuntimeError,"computeUnbiasedCorrection is implemented for order <=4 only!") - if type(weightsOrN).__name__ not in ['int','int8','int16','int64','int32']: - if order == 2: - V1, v1Square, V2 = self.__computeVp(1, weightsOrN), self.__computeVp(1, weightsOrN)**2.0, self.__computeVp(2, weightsOrN) - corrFactor = v1Square/(v1Square-V2) - elif order == 3: - V1, v1Cubic, V2, V3 = self.__computeVp(1, weightsOrN), self.__computeVp(1, weightsOrN)**3.0, self.__computeVp(2, weightsOrN), self.__computeVp(3, weightsOrN) - corrFactor = v1Cubic/(v1Cubic-3.0*V2*V1+2.0*V3) - elif order == 4: - V1, v1Square, V2, V3, V4 = self.__computeVp(1, weightsOrN), self.__computeVp(1, weightsOrN)**2.0, self.__computeVp(2, weightsOrN), self.__computeVp(3, weightsOrN), self.__computeVp(4, weightsOrN) - numer1 = v1Square*(v1Square**2.0-3.0*v1Square*V2+2.0*V1*V3+3.0*V2**2.0-3.0*V4) - numer2 = 3.0*v1Square*(2.0*v1Square*V2-2.0*V1*V3-3.0*V2**2.0+3.0*V4) - denom = (v1Square-V2)*(v1Square**2.0-6.0*v1Square*V2+8.0*V1*V3+3.0*V2**2.0-6.0*V4) - corrFactor = numer1/denom ,numer2/denom - else: - if order == 2: - corrFactor = float(weightsOrN)/(float(weightsOrN)-1.0) - elif order == 3: - corrFactor = (float(weightsOrN)**2.0)/((float(weightsOrN)-1)*(float(weightsOrN)-2)) - elif order == 4: - corrFactor = (float(weightsOrN)*(float(weightsOrN)**2.0-2.0*float(weightsOrN)+3.0))/((float(weightsOrN)-1)*(float(weightsOrN)-2)*(float(weightsOrN)-3)),(3.0*float(weightsOrN)*(2.0*float(weightsOrN)-3.0))/((float(weightsOrN)-1)*(float(weightsOrN)-2)*(float(weightsOrN)-3)) - return corrFactor - def _computeKurtosis(self, arrayIn, expValue, variance, pbWeight=None, dim=None): - """ - Method to compute the Kurtosis (fisher) of an array of observations - @ In, arrayIn, xarray.Dataset, the dataset from which the Kurtosis needs to be estimated - @ In, expValue, xarray.Dataset, expected value of arrayIn - @ In, variance, xarray.Dataset, variance of arrayIn - @ In, pbWeight, xarray.DataSet, optional, the reliability weights that correspond to the values in 'array'. - If not present, an unweighted approach is used - @ Out, result, xarray.Dataset, the Kurtosis of the dataset arrayIn. - """ - if dim is None: - dim = self.sampleTag - vr = self._computePower(2.0, variance) - if pbWeight is not None: - unbiasCorr = self.__computeUnbiasedCorrection(4,pbWeight) if not self.biased else 1.0 - vp = 1.0/self.__computeVp(1,pbWeight) - p4 = ((arrayIn - expValue)**4.0 * pbWeight).sum(dim=dim) - if not self.biased: - p2 = ((arrayIn - expValue)**2.0 * pbWeight).sum(dim=dim) - result = -3.0 + (p4*unbiasCorr[0]*vp - (p2*vp)**2.0 * unbiasCorr[1]) / vr - else: - result = -3.0 + (p4 * vp * unbiasCorr) / vr - else: - unbiasCorr = self.__computeUnbiasedCorrection(4,int(arrayIn.sizes[dim])) if not self.biased else 1.0 - vp = 1.0 / arrayIn.sizes[dim] - p4 = ((arrayIn - expValue)**4.0).sum(dim=dim) - if not self.biased: - p2 = (arrayIn - expValue).var(dim=dim) - result = -3.0 + (p4*unbiasCorr[0]*vp-p2**2.0*unbiasCorr[1]) / vr else: - result = -3.0 + (p4*unbiasCorr*vp) / vr - return result - - def _computeSkewness(self, arrayIn, expValue, variance, pbWeight=None, dim=None): - """ - Method to compute the skewness of an array of observations - @ In, arrayIn, xarray.Dataset, the dataset from which the skewness needs to be estimated - @ In, expValue, xarray.Dataset, expected value of arrayIn - @ In, variance, xarray.Dataset, variance value of arrayIn - @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. - If not present, an unweighted approach is used - @ Out, result, xarray.Dataset, the skewness of the dataset arrayIn - """ - if dim is None: - dim = self.sampleTag - vr = self._computePower(1.5, variance) - if pbWeight is not None: - unbiasCorr = self.__computeUnbiasedCorrection(3,pbWeight) if not self.biased else 1.0 - vp = 1.0/self.__computeVp(1,pbWeight) - result = ((arrayIn - expValue)**3 * pbWeight).sum(dim=dim) * vp * unbiasCorr / vr - else: - unbiasCorr = self.__computeUnbiasedCorrection(3,int(arrayIn.sizes[dim])) if not self.biased else 1.0 - vp = 1.0 / arrayIn.sizes[dim] - result = ((arrayIn - expValue)**3).sum(dim=dim) * vp * unbiasCorr / vr - return result - - def _computeVariance(self, arrayIn, expValue, pbWeight=None, dim = None): - """ - Method to compute the Variance (fisher) of an array of observations - @ In, arrayIn, xarray.Dataset, the dataset from which the Variance needs to be estimated - @ In, expValue, xarray.Dataset, expected value of arrayIn - @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. - If not present, an unweighted approach is used - @ Out, result, xarray.Dataset, the Variance of the dataset arrayIn - """ - if dim is None: - dim = self.sampleTag - if pbWeight is not None: - unbiasCorr = self.__computeUnbiasedCorrection(2,pbWeight) if not self.biased else 1.0 - vp = 1.0/self.__computeVp(1,pbWeight) - result = ((arrayIn-expValue)**2 * pbWeight).sum(dim=dim) * vp * unbiasCorr - else: - unbiasCorr = self.__computeUnbiasedCorrection(2,int(arrayIn.sizes[dim])) if not self.biased else 1.0 - result = (arrayIn-expValue).var(dim=dim) * unbiasCorr - return result - - def _computeLowerPartialVariance(self, arrayIn, medValue, pbWeight=None, dim = None): - """ - Method to compute the lower partial variance of an array of observations - @ In, arrayIn, xarray.Dataset, the dataset from which the Variance needs to be estimated - @ In, medValue, xarray.Dataset, median value of arrayIn - @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. - If not present, an unweighted approach is used - @ Out, result, xarray.Dataset, the lower partial variance of the dataset arrayIn - """ - if dim is None: - dim = self.sampleTag - diff = (medValue-arrayIn).clip(min=0) - if pbWeight is not None: - vp = 1.0/self.__computeVp(1,pbWeight) - result = ((diff)**2 * pbWeight).sum(dim=dim) * vp - else: - result = diff.var(dim=dim) - return result + if tag not in childVals: + self.raiseAWarning('Unrecognized node in BasicStatistics "',tag,'" has been ignored!') - def _computeHigherPartialVariance(self, arrayIn, medValue, pbWeight=None, dim = None): - """ - Method to compute the higher partial variance of an array of observations - @ In, arrayIn, xarray.Dataset, the dataset from which the Variance needs to be estimated - @ In, medValue, xarray.Dataset, median value of arrayIn - @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. - If not present, an unweighted approach is used - @ Out, result, xarray.Dataset, the higher partial variance of the dataset arrayIn - """ - if dim is None: - dim = self.sampleTag - diff = (arrayIn-medValue).clip(min=0) - if pbWeight is not None: - vp = 1.0/self.__computeVp(1,pbWeight) - result = ((diff)**2 * pbWeight).sum(dim=dim) * vp - else: - result = diff.var(dim=dim) - return result - def _computeSigma(self,arrayIn,variance,pbWeight=None): - """ - Method to compute the sigma of an array of observations - @ In, arrayIn, xarray.Dataset, the dataset from which the standard deviation needs to be estimated - @ In, variance, xarray.Dataset, variance of arrayIn - @ In, pbWeight, xarray.Dataset, optional, the reliability weights that correspond to dataset arrayIn. - If not present, an unweighted approach is used - @ Out, sigma, xarray.Dataset, the sigma of the dataset of arrayIn - """ - return np.sqrt(variance) - def _computeWeightedPercentile(self,arrayIn,pbWeight,interpolation='linear',percent=[0.5]): - """ - Method to compute the weighted percentile in a array of data - @ In, arrayIn, list/numpy.array, the array of values from which the percentile needs to be estimated - @ In, pbWeight, list/numpy.array, the reliability weights that correspond to the values in 'array' - @ In, interpolation, str, 'linear' or 'midpoint' - @ In, percent, list/numpy.array, the percentile(s) that needs to be computed (between 0.01 and 1.0) - @ Out, result, list, the percentile(s) - """ - # only do the argsort once for all requested percentiles - idxs = np.argsort(np.asarray(list(zip(pbWeight,arrayIn)))[:,1]) - # Inserting [0.0,arrayIn[idxs[0]]] is needed when few samples are generated and - # a percentile that is < that the first pb weight is requested. Otherwise the median - # is returned. - sortedWeightsAndPoints = np.insert(np.asarray(list(zip(pbWeight[idxs],arrayIn[idxs]))),0,[0.0,arrayIn[idxs[0]]],axis=0) - weightsCDF = np.cumsum(sortedWeightsAndPoints[:,0]) - weightsCDF /= weightsCDF[-1] - if interpolation == 'linear': - result = np.interp(percent, weightsCDF, sortedWeightsAndPoints[:, 1]).tolist() - elif interpolation == 'midpoint': - result = [self._computeSingleWeightedPercentile(pct, weightsCDF, sortedWeightsAndPoints) for pct in percent] - return result - def _computeSingleWeightedPercentile(self, pct, weightsCDF, sortedWeightsAndPoints): - """ - Method to compute a single percentile - @ In, pct, float, the percentile - @ In, weightsCDF, numpy.array, the cumulative sum of weights (CDF) - @ In, sortedWeightsAndPoints, numpy.array, array of weights and data points - @ Out, result, float, the percentile - """ - # This step returns the index of the array which is < than the percentile, because - # the insertion create another entry, this index should shift to the bigger side - indexL = utils.first(np.asarray(weightsCDF >= pct).nonzero())[0] - # This step returns the indices (list of index) of the array which is > than the percentile - indexH = utils.first(np.asarray(weightsCDF > pct).nonzero()) - try: - # if the indices exists that means the desired percentile lies between two data points - # with index as indexL and indexH[0]. Calculate the midpoint of these two points - result = 0.5*(sortedWeightsAndPoints[indexL,1]+sortedWeightsAndPoints[indexH[0],1]) - except IndexError: - result = sortedWeightsAndPoints[indexL,1] - return result def __runLocal(self, inputData): """ @@ -683,867 +267,13 @@ def __runLocal(self, inputData): @ Out, outputSet or outputDict, xarray.Dataset or dict, dataset or dictionary containing the results """ inputDataset, pbWeights = inputData[0], inputData[1] - #storage dictionary for skipped metrics - self.skipped = {} - #construct a dict of required computations - needed = dict((metric,{'targets':set(),'percent':set(),'interpolation':''}) for metric in self.scalarVals) - needed.update(dict((metric,{'targets':set(),'features':set()}) for metric in self.vectorVals)) - for metric, params in self.toDo.items(): - if metric in self.scalarVals + self.vectorVals: - for entry in params: - needed[metric]['targets'].update(entry['targets']) - try: - needed[metric]['features'].update(entry['features']) - except KeyError: - pass - try: - needed[metric]['percent'].update(entry['percent']) - except KeyError: - pass - try: - needed[metric]['interpolation'] = entry['interpolation'] - except KeyError: - pass - # variable | needs | needed for - # -------------------------------------------------------------------- - # skewness needs | expectedValue,variance | - # kurtosis needs | expectedValue,variance | - # median needs | | lowerPartialVariance, higherPartialVariance - # percentile needs | expectedValue,sigma | - # maximum needs | | - # minimum needs | | - # covariance needs | | pearson,VarianceDependentSensitivity,NormalizedSensitivity - # NormalizedSensitivity | covariance,VarDepSens | - # VarianceDependentSensitivity | covariance | NormalizedSensitivity - # sensitivity needs | | - # pearson needs | covariance | - # sigma needs | variance | variationCoefficient - # variance | expectedValue | sigma, skewness, kurtosis - # expectedValue | | variance, variationCoefficient, skewness, kurtosis, - # | | lowerPartialVariance, higherPartialVariance - # lowerPartialVariance needs | expectedValue,median | lowerPartialSigma - # lowerPartialSigma needs | lowerPartialVariance | - # higherPartialVariance needs | expectedValue,median | higherPartialSigma - # higherPartialSigma needs | higherPartialVariance | - - # update needed dictionary when standard errors are requested - needed['expectedValue']['targets'].update(needed['sigma']['targets']) - needed['expectedValue']['targets'].update(needed['variationCoefficient']['targets']) - needed['expectedValue']['targets'].update(needed['variance']['targets']) - needed['expectedValue']['targets'].update(needed['median']['targets']) - needed['expectedValue']['targets'].update(needed['skewness']['targets']) - needed['expectedValue']['targets'].update(needed['kurtosis']['targets']) - needed['expectedValue']['targets'].update(needed['NormalizedSensitivity']['targets']) - needed['expectedValue']['targets'].update(needed['NormalizedSensitivity']['features']) - needed['expectedValue']['targets'].update(needed['percentile']['targets']) - needed['sigma']['targets'].update(needed['expectedValue']['targets']) - needed['sigma']['targets'].update(needed['percentile']['targets']) - needed['variance']['targets'].update(needed['sigma']['targets']) - needed['lowerPartialVariance']['targets'].update(needed['lowerPartialSigma']['targets']) - needed['higherPartialVariance']['targets'].update(needed['higherPartialSigma']['targets']) - needed['median']['targets'].update(needed['lowerPartialVariance']['targets']) - needed['median']['targets'].update(needed['higherPartialVariance']['targets']) - needed['covariance']['targets'].update(needed['NormalizedSensitivity']['targets']) - needed['covariance']['features'].update(needed['NormalizedSensitivity']['features']) - needed['VarianceDependentSensitivity']['targets'].update(needed['NormalizedSensitivity']['targets']) - needed['VarianceDependentSensitivity']['features'].update(needed['NormalizedSensitivity']['features']) - needed['covariance']['targets'].update(needed['pearson']['targets']) - needed['covariance']['features'].update(needed['pearson']['features']) - needed['covariance']['targets'].update(needed['VarianceDependentSensitivity']['targets']) - needed['covariance']['features'].update(needed['VarianceDependentSensitivity']['features']) - - for metric, params in needed.items(): - needed[metric]['targets'] = list(params['targets']) - try: - needed[metric]['features'] = list(params['features']) - except KeyError: - pass - - # - # BEGIN actual calculations - # - - calculations = {} - ################# - # SCALAR VALUES # - ################# - # - # samples - # - self.sampleSize = inputDataset.sizes[self.sampleTag] - metric = 'samples' - if len(needed[metric]['targets']) > 0: - self.raiseADebug('Starting "'+metric+'"...') - if self.dynamic: - nt = inputDataset.sizes[self.pivotParameter] - sampleMat = np.zeros((len(self.parameters['targets']), len(self.pivotValue))) - sampleMat.fill(self.sampleSize) - samplesDA = xr.DataArray(sampleMat,dims=('targets', self.pivotParameter), coords={'targets':self.parameters['targets'], self.pivotParameter:self.pivotValue}) - else: - sampleMat = np.zeros(len(self.parameters['targets'])) - sampleMat.fill(self.sampleSize) - samplesDA = xr.DataArray(sampleMat,dims=('targets'), coords={'targets':self.parameters['targets']}) - self.calculations[metric] = samplesDA - calculations[metric] = samplesDA - # - # expected value - # - metric = 'expectedValue' - if len(needed[metric]['targets']) > 0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - if self.pbPresent: - relWeight = pbWeights[list(needed[metric]['targets'])] - equivalentSize = self.__computeEquivalentSampleSize(relWeight) - dataSet = dataSet * relWeight - expectedValueDS = dataSet.sum(dim = self.sampleTag) - calculations['equivalentSamples'] = equivalentSize - else: - expectedValueDS = dataSet.mean(dim = self.sampleTag) - self.calculations[metric] = expectedValueDS - calculations[metric] = expectedValueDS - # - # variance - # - metric = 'variance' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - meanSet = calculations['expectedValue'][list(needed[metric]['targets'])] - relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None - varianceDS = self._computeVariance(dataSet,meanSet,pbWeight=relWeight,dim=self.sampleTag) - calculations[metric] = varianceDS - # - # sigma - # - metric = 'sigma' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - sigmaDS = self._computePower(0.5,calculations['variance'][list(needed[metric]['targets'])]) - self.calculations[metric] = sigmaDS - calculations[metric] = sigmaDS - # - # coeff of variation (sigma/mu) - # - metric = 'variationCoefficient' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - calculations[metric] = calculations['sigma'][needed[metric]['targets']] / calculations['expectedValue'][needed[metric]['targets']] - # - # skewness - # - metric = 'skewness' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - meanSet = calculations['expectedValue'][list(needed[metric]['targets'])] - varianceSet = calculations['variance'][list(needed[metric]['targets'])] - relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None - calculations[metric] = self._computeSkewness(dataSet,meanSet,varianceSet,pbWeight=relWeight,dim=self.sampleTag) - # - # kurtosis - # - metric = 'kurtosis' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - meanSet = calculations['expectedValue'][list(needed[metric]['targets'])] - varianceSet = calculations['variance'][list(needed[metric]['targets'])] - relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None - calculations[metric] = self._computeKurtosis(dataSet,meanSet,varianceSet,pbWeight=relWeight,dim=self.sampleTag) - # - # median - # - metric = 'median' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - if self.pbPresent: - # if all weights are the same, calculate percentile with xarray, no need for _computeWeightedPercentile - allSameWeight = True - for target in needed[metric]['targets']: - targWeight = relWeight[target].values - if targWeight.min() != targWeight.max(): - allSameWeight = False - if allSameWeight: - medianSet = dataSet.median(dim=self.sampleTag) - else: - medianSet = xr.Dataset() - relWeight = pbWeights[list(needed[metric]['targets'])] - for target in needed[metric]['targets']: - targWeight = relWeight[target].values - targDa = dataSet[target] - if self.pivotParameter in targDa.sizes.keys(): - quantile = [self._computeWeightedPercentile(group.values,targWeight,needed[metric]['interpolation'],percent=[0.5])[0] for label,group in targDa.groupby(self.pivotParameter)] - else: - quantile = self._computeWeightedPercentile(targDa.values,targWeight,needed[metric]['interpolation'],percent=[0.5])[0] - if self.pivotParameter in targDa.sizes.keys(): - da = xr.DataArray(quantile,dims=(self.pivotParameter),coords={self.pivotParameter:self.pivotValue}) - else: - da = xr.DataArray(quantile) - medianSet[target] = da - else: - medianSet = dataSet.median(dim=self.sampleTag) - self.calculations[metric] = medianSet - calculations[metric] = medianSet - # - # lowerPartialVariance - # - metric = 'lowerPartialVariance' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - medianSet = calculations['median'][list(needed[metric]['targets'])] - relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None - lowerPartialVarianceDS = self._computeLowerPartialVariance(dataSet,medianSet,pbWeight=relWeight,dim=self.sampleTag) - - calculations[metric] = lowerPartialVarianceDS - # - # lowerPartialSigma - # - metric = 'lowerPartialSigma' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - lpsDS = self._computePower(0.5,calculations['lowerPartialVariance'][list(needed[metric]['targets'])]) - calculations[metric] = lpsDS - # - # higherPartialVariance - # - metric = 'higherPartialVariance' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - medianSet = calculations['median'][list(needed[metric]['targets'])] - relWeight = pbWeights[list(needed[metric]['targets'])] if self.pbPresent else None - higherPartialVarianceDS = self._computeHigherPartialVariance(dataSet,medianSet,pbWeight=relWeight,dim=self.sampleTag) - - calculations[metric] = lowerPartialVarianceDS - # - # higherPartialSigma - # - metric = 'higherPartialSigma' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - hpsDS = self._computePower(0.5,calculations['higherPartialVariance'][list(needed[metric]['targets'])]) - calculations[metric] = hpsDS - - ############################################################ - # Begin Standard Error Calculations - # - # Reference for standard error calculations (including percentile): - # B. Harding, C. Tremblay and D. Cousineau, "Standard errors: A review and evaluation of - # standard error estimators using Monte Carlo simulations", The Quantitative Methods of - # Psychology, Vol. 10, No. 2 (2014) - ############################################################ - metric = 'expectedValue' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting calculate standard error on"'+metric+'"...') - if self.pbPresent: - factor = self._computePower(0.5,calculations['equivalentSamples']) - else: - factor = np.sqrt(self.sampleSize) - calculations[metric+'_ste'] = calculations['sigma'][list(needed[metric]['targets'])]/factor - - metric = 'variance' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting calculate standard error on "'+metric+'"...') - varList = list(needed[metric]['targets']) - if self.pbPresent: - en = calculations['equivalentSamples'][varList] - factor = 2.0 /(en - 1.0) - factor = self._computePower(0.5,factor) - else: - factor = np.sqrt(2.0/(float(self.sampleSize) - 1.0)) - calculations[metric+'_ste'] = calculations['sigma'][varList]**2 * factor - - metric = 'sigma' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting calculate standard error on "'+metric+'"...') - varList = list(needed[metric]['targets']) - if self.pbPresent: - en = calculations['equivalentSamples'][varList] - factor = 2.0 * (en - 1.0) - factor = self._computePower(0.5,factor) - else: - factor = np.sqrt(2.0 * (float(self.sampleSize) - 1.0)) - calculations[metric+'_ste'] = calculations['sigma'][varList] / factor - - metric = 'median' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting calculate standard error on "'+metric+'"...') - varList = list(needed[metric]['targets']) - calculations[metric+'_ste'] = calculations['expectedValue_ste'][varList] * np.sqrt(np.pi/2.0) - - metric = 'skewness' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting calculate standard error on "'+metric+'"...') - varList = list(needed[metric]['targets']) - if self.pbPresent: - en = calculations['equivalentSamples'][varList] - factor = 6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.)) - factor = self._computePower(0.5,factor) - calculations[metric+'_ste'] = xr.full_like(calculations[metric],1.0) * factor - else: - en = float(self.sampleSize) - factor = np.sqrt(6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.))) - calculations[metric+'_ste'] = xr.full_like(calculations[metric],factor) - - metric = 'kurtosis' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting calculate standard error on "'+metric+'"...') - varList = list(needed[metric]['targets']) - if self.pbPresent: - en = calculations['equivalentSamples'][varList] - factor1 = self._computePower(0.5,6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.))) - factor2 = self._computePower(0.5,(en**2-1.)/((en-3.0)*(en+5.0))) - factor = 2.0 * factor1 * factor2 - calculations[metric+'_ste'] = xr.full_like(calculations[metric],1.0) * factor - else: - en = float(self.sampleSize) - factor = 2.0 * np.sqrt(6.*en*(en-1.)/((en-2.)*(en+1.)*(en+3.)))*np.sqrt((en**2-1.)/((en-3.0)*(en+5.0))) - calculations[metric+'_ste'] = xr.full_like(calculations[metric],factor) - ############################################################ - # End of Standard Error Calculations - ############################################################ - # - # maximum - # - metric = 'maximum' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - calculations[metric] = dataSet.max(dim=self.sampleTag) - # - # minimum - # - metric = 'minimum' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - calculations[metric] = dataSet.min(dim=self.sampleTag) - # - # percentile, this metric is handled differently - # - metric = 'percentile' - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - dataSet = inputDataset[list(needed[metric]['targets'])] - percent = list(needed[metric]['percent']) - # are there probability weights associated with the data? - if self.pbPresent: - relWeight = pbWeights[list(needed[metric]['targets'])] - # if all weights are the same, calculate percentile with xarray, no need for _computeWeightedPercentile - allSameWeight = True - for target in needed[metric]['targets']: - targWeight = relWeight[target].values - if targWeight.min() != targWeight.max(): - allSameWeight = False - if allSameWeight: - # all weights are the same, percentile can be calculated with xarray.DataSet - percentileSet = dataSet.quantile(percent,dim=self.sampleTag,interpolation=needed[metric]['interpolation']) - percentileSet = percentileSet.rename({'quantile': 'percent'}) - else: - # probability weights are not all the same - # xarray does not have capability to calculate weighted quantiles at present - # implement our own solution - percentileSet = xr.Dataset() - for target in needed[metric]['targets']: - targWeight = relWeight[target].values - targDa = dataSet[target] - if self.pivotParameter in targDa.sizes.keys(): - quantile = [] - for label, group in targDa.groupby(self.pivotParameter): - qtl = self._computeWeightedPercentile(group.values, targWeight, needed[metric]['interpolation'], percent=percent) - quantile.append(qtl) - da = xr.DataArray(quantile, dims=(self.pivotParameter, 'percent'), coords={'percent': percent, self.pivotParameter: self.pivotValue}) - else: - quantile = self._computeWeightedPercentile(targDa.values, targWeight, needed[metric]['interpolation'], percent=percent) - da = xr.DataArray(quantile, dims=('percent'), coords={'percent': percent}) - - percentileSet[target] = da - - # TODO: remove when complete - # interpolation: {'linear', 'lower', 'higher','midpoint','nearest'}, do not try to use 'linear' or 'midpoint' - # The xarray.Dataset.where() will not return the corrrect solution - # 'lower' is used for consistent - # using xarray.Dataset.sel(**{'quantile':reqPercent}) to retrieve the quantile values - #dataSetWeighted = dataSet * relWeight - #percentileSet = dataSet.where(dataSetWeighted==dataSetWeighted.quantile(percent,dim=self.sampleTag,interpolation='lower')).mean(self.sampleTag) - else: - percentileSet = dataSet.quantile(percent,dim=self.sampleTag,interpolation=needed[metric]['interpolation']) - percentileSet = percentileSet.rename({'quantile':'percent'}) - calculations[metric] = percentileSet - - # because percentile is different, calculate standard error here - # standard error calculation uses the standard normal formulation for speed - self.raiseADebug('Starting calculate standard error on "'+metric+'"...') - norm = stats.norm - factor = np.sqrt(np.asarray(percent)*(1.0 - np.asarray(percent)))/norm.pdf(norm.ppf(percent)) - sigmaAdjusted = calculations['sigma'][list(needed[metric]['targets'])]/np.sqrt(calculations['equivalentSamples'][list(needed[metric]['targets'])]) - sigmaAdjusted = sigmaAdjusted.expand_dims(dim={'percent': percent}) - factor = xr.DataArray(data=factor, dims='percent', coords={'percent': percent}) - calculations[metric + '_ste'] = sigmaAdjusted*factor - - # # TODO: this is the KDE method, it is a more accurate method of calculating standard error - # # for percentile, but the computation time is too long. IF this computation can be sped up, - # # implement it here: - # percentileSteSet = xr.Dataset() - # calculatedPercentiles = calculations[metric] - # relWeight = pbWeights[list(needed[metric]['targets'])] - # for target in needed[metric]['targets']: - # targWeight = relWeight[target].values - # en = calculations['equivalentSamples'][target].values - # targDa = dataSet[target] - # if self.pivotParameter in targDa.sizes.keys(): - # percentileSte = np.zeros((len(self.pivotValue), len(percent))) # array - # for i, (label, group) in enumerate(targDa.groupby(self.pivotParameter)): # array - # if group.values.min() == group.values.max(): - # subPercentileSte = np.array([0.0]*len(percent)) - # else: - # # get KDE - # kde = stats.gaussian_kde(group.values, weights=targWeight) - # vals = calculatedPercentiles[target].sel(**{'percent': percent, self.pivotParameter: label}).values - # factor = np.sqrt(np.asarray(percent)*(1.0 - np.asarray(percent))/en) - # subPercentileSte = factor/kde(vals) - # percentileSte[i, :] = subPercentileSte - # da = xr.DataArray(percentileSte, dims=(self.pivotParameter, 'percent'), coords={self.pivotParameter: self.pivotValue, 'percent': percent}) - # percentileSteSet[target] = da - # else: - # calcPercentiles = calculatedPercentiles[target] - # if targDa.values.min() == targDa.values.max(): - # # distribution is a delta function, so no KDE construction - # percentileSte = list(np.zeros(calcPercentiles.shape)) - # else: - # # get KDE - # kde = stats.gaussian_kde(targDa.values, weights=targWeight) - # factor = np.sqrt(np.array(percent)*(1.0 - np.array(percent))/en) - # percentileSte = list(factor/kde(calcPercentiles.values)) - # da = xr.DataArray(percentileSte, dims=('percent'), coords={'percent': percent}) - # percentileSteSet[target] = da - # calculations[metric+'_ste'] = percentileSteSet - - def startVector(metric): - """ - Common method among all metrics for establishing parameters - @ In, metric, string, the name of the statistics metric to calculate - @ Out, targets, list(str), list of target parameter names (evaluate metrics for these) - @ Out, features, list(str), list of feature parameter names (evaluate with respect to these) - @ Out, skip, bool, if True it means either features or parameters were missing, so don't calculate anything - """ - # default to skipping, change that if we find criteria - targets = [] - features = [] - skip = True - if len(needed[metric]['targets'])>0: - self.raiseADebug('Starting "'+metric+'"...') - targets = list(needed[metric]['targets']) - features = list(needed[metric]['features']) - skip = False #True only if we don't have targets and features - if skip: - if metric not in self.skipped.keys(): - self.skipped[metric] = True - return targets,features,skip - - ################# - # VECTOR VALUES # - ################# - # - # sensitivity matrix - # - metric = 'sensitivity' - targets,features,skip = startVector(metric) - #NOTE sklearn expects the transpose of what we usually do in RAVEN, so #samples by #features - if not skip: - #for sensitivity matrix, we don't use numpy/scipy methods to calculate matrix operations, - #so we loop over targets and features - params = list(set(targets).union(set(features))) - dataSet = inputDataset[params] - relWeight = pbWeights[params] if self.pbPresent else None - intersectionSet = set(targets) & set(features) - if self.pivotParameter in dataSet.sizes.keys(): - dataSet = dataSet.to_array().transpose(self.pivotParameter,self.sampleTag,'variable') - featSet = dataSet.sel(**{'variable':features}).values - targSet = dataSet.sel(**{'variable':targets}).values - pivotVals = dataSet.coords[self.pivotParameter].values - da = None - for i in range(len(pivotVals)): - ds = self.sensitivityCalculation(features,targets,featSet[i,:,:],targSet[i,:,:],intersectionSet) - da = ds if da is None else xr.concat([da,ds], dim=self.pivotParameter) - da.coords[self.pivotParameter] = pivotVals - else: - # construct target and feature matrices - dataSet = dataSet.to_array().transpose(self.sampleTag,'variable') - featSet = dataSet.sel(**{'variable':features}).values - targSet = dataSet.sel(**{'variable':targets}).values - da = self.sensitivityCalculation(features,targets,featSet,targSet,intersectionSet) - calculations[metric] = da - # - # covariance matrix - # - metric = 'covariance' - targets,features,skip = startVector(metric) - if not skip: - # because the C implementation is much faster than picking out individual values, - # we do the full covariance matrix with all the targets and features. - # FIXME adding an alternative for users to choose pick OR do all, defaulting to something smart - # dependent on the percentage of the full matrix desired, would be better. - # IF this is fixed, make sure all the features and targets are requested for all the metrics - # dependent on this metric - params = list(set(targets).union(set(features))) - dataSet = inputDataset[params] - relWeight = pbWeights[params] if self.pbPresent else None - if self.pbPresent: - fact = (self.__computeUnbiasedCorrection(2, self.realizationWeight)).to_array().values if not self.biased else 1.0 - meanSet = (dataSet * relWeight).sum(dim = self.sampleTag) - else: - meanSet = dataSet.mean(dim = self.sampleTag) - fact = 1.0 / (float(dataSet.sizes[self.sampleTag]) - 1.0) if not self.biased else 1.0 / float(dataSet.sizes[self.sampleTag]) - targVars = list(dataSet.data_vars) - varianceSet = self._computeVariance(dataSet,meanSet,pbWeight=relWeight,dim=self.sampleTag) - dataSet = dataSet - meanSet - if self.pivotParameter in dataSet.sizes.keys(): - ds = None - paramDA = dataSet.to_array().transpose(self.pivotParameter,'variable',self.sampleTag).values - varianceDA = varianceSet[targVars].to_array().transpose(self.pivotParameter,'variable').values - pivotVals = dataSet.coords[self.pivotParameter].values - for i in range(len(pivotVals)): - # construct target and feature matrices - paramSamples = paramDA[i,...] - da = self.covarianceCalculation(paramDA[i,...],fact,varianceDA[i,:],targVars) - ds = da if ds is None else xr.concat([ds,da], dim=self.pivotParameter) - ds.coords[self.pivotParameter] = pivotVals - calculations[metric] = ds - else: - # construct target and feature matrices - paramSamples = dataSet.to_array().transpose('variable',self.sampleTag).values - varianceDA = varianceSet[targVars].to_array().values - da = self.covarianceCalculation(paramSamples,fact,varianceDA,targVars) - calculations[metric] = da - - def getCovarianceSubset(desired): - """ - @ In, desired, list(str), list of parameters to extract from covariance matrix - @ Out, reducedCov, xarray.DataArray, reduced covariance matrix - """ - if self.pivotParameter in desired: - self.raiseAnError(RuntimeError, 'The pivotParameter "{}" is among the parameters requested for performing statistics. Please remove!'.format(self.pivotParameter)) - reducedCov = calculations['covariance'].sel(**{'targets':desired,'features':desired}) - return reducedCov - # - # pearson matrix - # - # see comments in covariance for notes on C implementation - metric = 'pearson' - targets,features,skip = startVector(metric) - if not skip: - params = list(set(targets).union(set(features))) - reducedCovar = getCovarianceSubset(params) - targCoords = reducedCovar.coords['targets'].values - if self.pivotParameter in reducedCovar.sizes.keys(): - pivotCoords = reducedCovar.coords[self.pivotParameter].values - ds = None - for label, group in reducedCovar.groupby(self.pivotParameter): - corrMatrix = self.corrCoeff(group.values) - da = xr.DataArray(corrMatrix, dims=('targets','features'), coords={'targets':targCoords,'features':targCoords}) - ds = da if ds is None else xr.concat([ds,da], dim=self.pivotParameter) - ds.coords[self.pivotParameter] = pivotCoords - calculations[metric] = ds - else: - corrMatrix = self.corrCoeff(reducedCovar.values) - da = xr.DataArray(corrMatrix, dims=('targets','features'), coords={'targets':targCoords,'features':targCoords}) - calculations[metric] = da - # - # spearman matrix - # - # see RAVEN theory manual for a detailed explaination - # of the formulation used here - # - metric = 'spearman' - targets,features,skip = startVector(metric) - #NOTE sklearn expects the transpose of what we usually do in RAVEN, so #samples by #features - if not skip: - #for spearman matrix, we don't use numpy/scipy methods to calculate matrix operations, - #so we loop over targets and features - params = list(set(targets).union(set(features))) - dataSet = inputDataset[params] - relWeight = pbWeights[params] if self.pbPresent else None - if self.pivotParameter in dataSet.sizes.keys(): - dataSet = dataSet.to_array().transpose(self.pivotParameter,self.sampleTag,'variable') - featSet = dataSet.sel(**{'variable':features}).values - targSet = dataSet.sel(**{'variable':targets}).values - pivotVals = dataSet.coords[self.pivotParameter].values - da = None - for i in range(len(pivotVals)): - ds = self.spearmanCorrelation(features,targets,featSet[i,:,:],targSet[i,:,:],relWeight) - da = ds if da is None else xr.concat([da,ds], dim=self.pivotParameter) - da.coords[self.pivotParameter] = pivotVals - else: - # construct target and feature matrices - dataSet = dataSet.to_array().transpose(self.sampleTag,'variable') - featSet = dataSet.sel(**{'variable':features}).values - targSet = dataSet.sel(**{'variable':targets}).values - da = self.spearmanCorrelation(features,targets,featSet,targSet,relWeight) - calculations[metric] = da - # - # VarianceDependentSensitivity matrix - # The formula for this calculation is coming from: http://www.math.uah.edu/stat/expect/Matrices.html - # The best linear predictor: L(Y|X) = expectedValue(Y) + cov(Y,X) * [vc(X)]^(-1) * [X-expectedValue(X)] - # where Y is a vector of outputs, and X is a vector of inputs, cov(Y,X) is the covariance matrix of Y and X, - # vc(X) is the covariance matrix of X with itself. - # The variance dependent sensitivity matrix is defined as: cov(Y,X) * [vc(X)]^(-1) - metric = 'VarianceDependentSensitivity' - targets,features,skip = startVector(metric) - if not skip: - params = list(set(targets).union(set(features))) - reducedCovar = getCovarianceSubset(params) - targCoords = reducedCovar.coords['targets'].values - if self.pivotParameter in reducedCovar.sizes.keys(): - pivotCoords = reducedCovar.coords[self.pivotParameter].values - ds = None - for label, group in reducedCovar.groupby(self.pivotParameter): - da = self.varianceDepSenCalculation(targCoords,group.values) - ds = da if ds is None else xr.concat([ds,da], dim=self.pivotParameter) - ds.coords[self.pivotParameter] = pivotCoords - calculations[metric] = ds - else: - da = self.varianceDepSenCalculation(targCoords,reducedCovar.values) - calculations[metric] = da - - # - # Normalized variance dependent sensitivity matrix - # variance dependent sensitivity normalized by the mean (% change of output)/(% change of input) - # - metric = 'NormalizedSensitivity' - targets,features,skip = startVector(metric) - if not skip: - params = list(set(targets).union(set(features))) - reducedSen = calculations['VarianceDependentSensitivity'].sel(**{'targets':params,'features':params}) - meanDA = calculations['expectedValue'][params].to_array() - meanDA = meanDA.rename({'variable':'targets'}) - reducedSen /= meanDA - meanDA = meanDA.rename({'targets':'features'}) - reducedSen *= meanDA - calculations[metric] = reducedSen - - for metric, ds in calculations.items(): - if metric in self.scalarVals + self.steVals +['equivalentSamples'] and metric !='samples': - calculations[metric] = ds.to_array().rename({'variable':'targets'}) - # in here we fill the NaN with "nan". In this way, we are sure that even if - # there might be NaN in any raw for a certain timestep we do not drop the variable - # In the past, in a condition such as: - # time, A, B, C - # 0, 1, NaN, 1 - # 1, 1, 0.5, 1 - # 2, 1, 2.0, 2 - # the variable B would have been dropped (in the printing stage) - # with this modification, this should not happen anymore - outputSet = xr.Dataset(data_vars=calculations).fillna("nan") - - if self.outputDataset: - # Add 'RAVEN_sample_ID' to output dataset for consistence - if 'RAVEN_sample_ID' not in outputSet.sizes.keys(): - outputSet = outputSet.expand_dims('RAVEN_sample_ID') - outputSet['RAVEN_sample_ID'] = [0] - return outputSet - else: - outputDict = {} - for metric, requestList in self.toDo.items(): - for targetDict in requestList: - prefix = targetDict['prefix'].strip() - for target in targetDict['targets']: - if metric in self.scalarVals and metric != 'percentile': - varName = prefix + '_' + target - outputDict[varName] = np.atleast_1d(outputSet[metric].sel(**{'targets':target})) - steMetric = metric + '_ste' - if steMetric in self.steVals: - metaVar = prefix + '_ste_' + target - outputDict[metaVar] = np.atleast_1d(outputSet[steMetric].sel(**{'targets':target})) - elif metric == 'percentile': - for percent in targetDict['strPercent']: - varName = '_'.join([prefix,percent,target]) - percentVal = float(percent)/100. - outputDict[varName] = np.atleast_1d(outputSet[metric].sel(**{'targets':target,'percent':percentVal})) - steMetric = metric + '_ste' - if steMetric in self.steVals: - metaVar = '_'.join([prefix,percent,'ste',target]) - outputDict[metaVar] = np.atleast_1d(outputSet[steMetric].sel(**{'targets':target,'percent':percentVal})) - else: - #check if it was skipped for some reason - skip = self.skipped.get(metric, None) - if skip is not None: - self.raiseADebug('Metric',metric,'was skipped for parameters',targetDict,'! See warnings for details. Ignoring...') - continue - if metric in self.vectorVals: - for feature in targetDict['features']: - varName = '_'.join([prefix,target,feature]) - outputDict[varName] = np.atleast_1d(outputSet[metric].sel(**{'targets':target,'features':feature})) - if self.pivotParameter in outputSet.sizes.keys(): - outputDict[self.pivotParameter] = np.atleast_1d(self.pivotValue) - - return outputDict - - def corrCoeff(self, covM): - """ - This method calculates the correlation coefficient Matrix (pearson) for the given data. - Unbiased unweighted covariance matrix, weights is None, bias is 0 (default) - Biased unweighted covariance matrix, weights is None, bias is 1 - Unbiased weighted covariance matrix, weights is not None, bias is 0 - Biased weighted covariance matrix, weights is not None, bias is 1 - can be calcuated depending on the selection of the inputs. - @ In, covM, numpy.array, [#targets,#targets] covariance matrix - @ Out, covM, numpy.array, [#targets,#targets] correlation matrix - """ - try: - d = np.diag(covM) - except ValueError: - # scalar covariance - # nan if incorrect value (nan, inf, 0), 1 otherwise - return covM / covM - stdDev = np.sqrt(d) - covM /= stdDev[:,None] - covM /= stdDev[None,:] - return covM - - def sensitivityCalculation(self,featVars, targVars, featSamples, targSamples, intersectionSet): - """ - This method computes the sensitivity coefficients based on the SciKitLearn LinearRegression method - @ In, featVars, list, list of feature variables - @ In, targVars, list, list of target variables - @ In, featSamples, numpy.ndarray, [#samples, #features] array of features - @ In, targSamples, numpy.ndarray, [#samples, #targets] array of targets - @ In, intersectionSet, boolean, True if some target variables are in the list of features - @ Out, da, xarray.DataArray, contains the calculations of sensitivity coefficients - """ - from sklearn.linear_model import LinearRegression - if self.multipleFeatures: - # intersectionSet is flag that used to check the relationship between the features and targets. - # If True, part of the target variables are listed in teh feature set, then multivariate linear - # regression should not be used, and a loop over the target set is required. - # If False, which means there is no overlap between the target set and feature set. - # mutivariate linear regression can be used. However, for both cases, co-linearity check should be - # added for the feature set. ~ wangc - if not intersectionSet: - condNumber = np.linalg.cond(featSamples) - if condNumber > 30.: - self.raiseAWarning("Condition Number: {:10.4f} > 30.0. Detected SEVERE multicollinearity problem. Sensitivity might be incorrect!".format(condNumber)) - senMatrix = LinearRegression().fit(featSamples,targSamples).coef_ - else: - # Target variables are in feature variables list, multi-target linear regression can not be used - # Since the 'multi-colinearity' exists, we need to loop over target variables - # TODO: Some general methods need to be implemented in order to handle the 'multi-colinearity' -- wangc - senMatrix = np.zeros((len(targVars), len(featVars))) - for p, targ in enumerate(targVars): - ind = list(featVars).index(targ) if targ in featVars else None - if ind is not None: - featMat = np.delete(featSamples,ind,axis=1) - else: - featMat = featSamples - regCoeff = LinearRegression().fit(featMat, targSamples[:,p]).coef_ - condNumber = np.linalg.cond(featMat) - if condNumber > 30.: - self.raiseAWarning("Condition Number: {:10.4f} > 30.0. Detected SEVERE multicollinearity problem. Sensitivity might be incorrect!".format(condNumber)) - if ind is not None: - regCoeff = np.insert(regCoeff,ind,1.0) - senMatrix[p,:] = regCoeff - else: - senMatrix = np.zeros((len(targVars), len(featVars))) - for p, feat in enumerate(featVars): - regCoeff = LinearRegression().fit(featSamples[:,p].reshape(-1,1),targSamples).coef_ - senMatrix[:,p] = regCoeff[:,0] - da = xr.DataArray(senMatrix, dims=('targets','features'), coords={'targets':targVars,'features':featVars}) - - return da - - def covarianceCalculation(self,paramSamples,fact,variance,targVars): - """ - This method computes the covariance of given sample matrix - @ In, paramSamples, numpy.ndarray, [#parameters, #samples], array of parameters - @ In, fact, float, the unbiase correction factor - @ In, variance, numpy.ndarray, [#parameters], variance of parameters - @ In, targVars, list, the list of parameters - @ Out, da, xarray.DataArray, contains the calculations of covariance - """ - if self.pbPresent: - paramSamplesT = (paramSamples*self.realizationWeight['ProbabilityWeight'].values).T - else: - paramSamplesT = paramSamples.T - cov = np.dot(paramSamples, paramSamplesT.conj()) - cov *= fact - np.fill_diagonal(cov,variance) - da = xr.DataArray(cov, dims=('targets','features'), coords={'targets':targVars,'features':targVars}) - return da + return outputDict - def varianceDepSenCalculation(self,targCoords, cov): - """ - This method computes the covariance of given sample matrix - @ In, targCoords, list, the list of parameters - @ In, cov, numpy.ndarray, the covariance of parameters - @ Out, da, xarray.DataArray, contains the calculations of variance dependent sensitivities - """ - senMatrix = np.zeros((len(targCoords), len(targCoords))) - if self.multipleFeatures: - for p, param in enumerate(targCoords): - covX = np.delete(cov,p,axis=0) - covX = np.delete(covX,p,axis=1) - covYX = np.delete(cov[p,:],p) - sensCoef = np.dot(covYX,np.linalg.pinv(covX)) - sensCoef = np.insert(sensCoef,p,1.0) - senMatrix[p,:] = sensCoef - else: - for p, param in enumerate(targCoords): - covX = cov[p,p] - covYX = cov[:,p] - sensCoef = covYX / covX - senMatrix[:,p] = sensCoef - da = xr.DataArray(senMatrix, dims=('targets','features'), coords={'targets':targCoords,'features':targCoords}) - return da - def spearmanCorrelation(self, featVars, targVars, featSamples, targSamples, pbWeights): - """ - This method computes the spearman correlation coefficients - @ In, featVars, list, list of feature variables - @ In, targVars, list, list of target variables - @ In, featSamples, numpy.ndarray, [#samples, #features] array of features - @ In, targSamples, numpy.ndarray, [#samples, #targets] array of targets - @ In, pbWeights, dataset, probability weights - @ Out, da, xarray.DataArray, contains the calculations of spearman coefficients - """ - spearmanMat = np.zeros((len(targVars), len(featVars))) - wf, wt = None, None - # compute unbiased factor - if self.pbPresent: - fact = (self.__computeUnbiasedCorrection(2, self.realizationWeight)).to_array().values if not self.biased else 1.0 - vp = self.__computeVp(1,self.realizationWeight)['ProbabilityWeight'].values - varianceFactor = fact*(1.0/vp) - else: - fact = 1.0 / (float(featSamples.shape[0]) - 1.0) if not self.biased else 1.0 / float(featSamples.shape[0]) - varianceFactor = fact - for tidx, target in enumerate(targVars): - for fidx, feat in enumerate(featVars): - if self.pbPresent: - wf, wt = np.asarray(pbWeights[feat]), np.asarray(pbWeights[target]) - rankFeature, rankTarget = mathUtils.rankData(featSamples[:,fidx],wf), mathUtils.rankData(targSamples[:,tidx],wt) - # compute covariance of the ranked features - cov = np.cov(rankFeature, y=rankTarget, aweights=wt) - covF = np.cov(rankFeature,y=rankFeature, aweights=wf) - covT = np.cov(rankTarget,y=rankTarget, aweights=wf) - # apply correction factor (for biased or unbiased) (off diagonal) - cov[~np.eye(2,dtype=bool)] *= fact - covF[~np.eye(2,dtype=bool)] *= fact - covT[~np.eye(2,dtype=bool)] *= fact - # apply correction factor (for biased or unbiased) (diagonal) - cov[np.eye(2,dtype=bool)] *= varianceFactor - covF[~np.eye(2,dtype=bool)] *= varianceFactor - covT[~np.eye(2,dtype=bool)] *= varianceFactor - # now we can compute the pearson of such pairs - spearman = (cov / np.sqrt(covF * covT))[-1,0] - spearmanMat[tidx,fidx] = spearman - da = xr.DataArray(spearmanMat, dims=('targets','features'), coords={'targets':targVars,'features':featVars}) - return da def run(self, inputIn): """ diff --git a/ravenframework/Samplers/Grid.py b/ravenframework/Samplers/Grid.py index 89de6451a7..88327bbca2 100644 --- a/ravenframework/Samplers/Grid.py +++ b/ravenframework/Samplers/Grid.py @@ -100,7 +100,8 @@ def localInputAndChecks(self, xmlNode, paramInput): self.raiseAnError(IOError,'limit is not used in Grid sampler') self.limit = 1 # FIXME: THIS READ MORE XML MUST BE CONVERTED IN THE INPUTPARAMETER COLLECTOR!!!!!!! - self.gridEntity._readMoreXml(xmlNode,dimensionTags=["variable", "Distribution"], dimTagsPrefix={"Distribution": ""}) + self.gridEntity._handleInput(paramInput, dimensionTags=["variable", "Distribution"], dimTagsPrefix={"Distribution": ""}) + grdInfo = self.gridEntity.returnParameter("gridInfo") for axis, value in grdInfo.items(): self.gridInfo[axis] = value[0] diff --git a/ravenframework/utils/mathUtils.py b/ravenframework/utils/mathUtils.py index fefedabecb..6520e6b08f 100644 --- a/ravenframework/utils/mathUtils.py +++ b/ravenframework/utils/mathUtils.py @@ -17,7 +17,6 @@ created on 03/26/2015 @author: senrs """ - import math import functools import copy From a44ea7fab797bcb024547d150e23b89679e49091 Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Mon, 22 May 2023 10:23:58 -0600 Subject: [PATCH 03/10] added tests --- ravenframework/GridEntities.py | 66 +++ .../Models/PostProcessors/BasicStatistics.py | 25 +- .../Models/PostProcessors/LimitSurface.py | 13 - .../PostProcessors/PostProcessorInterface.py | 2 +- .../SubdomainBasicStatistics.py | 242 ++++------- ravenframework/Samplers/Grid.py | 1 - ravenframework/Samplers/LimitSurfaceSearch.py | 2 +- .../subdomainStats/subdomainSensitivity.csv | 5 + .../timedepSubdomainStats_dump.csv | 3 + .../timedepSubdomainStats_dump_0.csv | 4 + .../sensitivity_subdomain.xml | 381 ++++++++++++++++++ .../SubdomainBasicStatistics/tests | 18 + .../time_dep_subdomain.xml | 134 ++++++ 13 files changed, 709 insertions(+), 187 deletions(-) create mode 100644 tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainStats/subdomainSensitivity.csv create mode 100644 tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump.csv create mode 100644 tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump_0.csv create mode 100644 tests/framework/PostProcessors/SubdomainBasicStatistics/sensitivity_subdomain.xml create mode 100644 tests/framework/PostProcessors/SubdomainBasicStatistics/tests create mode 100644 tests/framework/PostProcessors/SubdomainBasicStatistics/time_dep_subdomain.xml diff --git a/ravenframework/GridEntities.py b/ravenframework/GridEntities.py index 6e188b3cce..91b500577f 100644 --- a/ravenframework/GridEntities.py +++ b/ravenframework/GridEntities.py @@ -18,6 +18,7 @@ """ import abc import sys +from numbers import Number # External Modules---------------------------------------------------------------------------------- import numpy as np @@ -389,6 +390,9 @@ def initialize(self, initDictionary=None): initDict = initDictionary if initDictionary is not None else {} computeCells = bool(initDict.get('computeCells', False)) self.constructTensor = bool(initDict['constructTensor']) if 'constructTensor' in initDict else False + # assert that compute cells are requested if and only if constructTensor is True + assert not(not self.constructTensor and computeCells), 'computeCells can be only computed if the tensor grid is activated!!' + if len(self.gridInitDict.keys()) != 0: readKeys = list(self.gridInitDict.keys()) if initDict is not None: @@ -702,6 +706,47 @@ def returnCoordinateFromIndex(self, multiDimIndex, returnDict=False, recastMetho return coordinates + def returnCellIdsWithCoordinates(self, returnDict=False, recastMethods={}): + """ + Method to return the container of the cell ids with the coordinates of the vertices in the real space (real coordinates) + @ In, returnDict, bool, optional, flag to request the output of the coordinates in dictionary format or not. + if True a dict ({cellID: {dimName1: + coordinate1,dimName2:coordinate2,etc} } is returned + if False a tuple is riturned {cellID:(coordinate1,coordinate2,etc)} + @ In, recastMethods, dict, optional, dictionary containing the methods that need to be used for trasforming the coordinates + ex. {'dimName1':[methodToTransformCoordinate,*args]} + @ Out, cellIdsAndVerteces, dict, dict containing the cell ids and the coordinates of the verteces + """ + # if no cell ids have been computed, it returns an empty dict + cellIdsAndVerteces = {} + for cid in self.gridContainer['cellIDs']: + cellIdsAndVerteces[cid] = [] + for vertex in self.gridContainer['cellIDs'][cid]: + cellIdsAndVerteces[cid].append(self.returnCoordinateFromIndex(vertex, returnDict, recastMethods)) + return cellIdsAndVerteces + + def returnCellsMidPoints(self, returnDict=False, recastMethods={}): + """ + Method to return the container of the cell ids with the coordinate of the mid point in the real space (real coordinate) + @ In, returnDict, bool, optional, flag to request the output of the coordinates in dictionary format or not. + if True a dict ({cellID: {dimName1: + coordinate1,dimName2:coordinate2,etc} } is returned + if False a tuple is riturned {cellID:(coordinate1,coordinate2,etc)} + @ In, recastMethods, dict, optional, dictionary containing the methods that need to be used for trasforming the coordinates + ex. {'dimName1':[methodToTransformCoordinate,*args]} + @ Out, cellIdsAndMidPoints, dict, dict containing the cell ids and the coordinates of the verteces + """ + cellsAndVerteces = self.returnCellIdsWithCoordinates(False, recastMethods) + # if no cell ids have been computed, it returns an empty dict + cellIdsAndMidPoints = {} + for cid in cellsAndVerteces: + midPoint = np.atleast_2d(cellsAndVerteces[cid]).mean(axis=0) + if returnDict: + cellIdsAndMidPoints[cid] = {d: midPoint[cnt] for cnt, d in enumerate(self.gridContainer['dimensionNames'])} + else: + cellIdsAndMidPoints[cid] = midPoint + return cellIdsAndMidPoints + def flush(self): """ Reset GridEntity attributes to allow rerunning a workflow @@ -1049,6 +1094,27 @@ def returnCoordinateFromIndex(self, multiDimNDIndex, returnDict=False, recastMet return node.get('grid').returnCoordinateFromIndex(multiDimIndex, returnDict, recastMethods) + def returnCoordinateFromIndex(self, multiDimNDIndex, returnDict=False, recastMethods={}): + """ + Method to return a point in the grid. This method will return the coordinates of the point is requested by multiDimIndex + In addition, it advances the iterator in order to point to the following coordinate + @ In, multiDimNDIndex, tuple, tuple containing the Id of the point needs to be returned (e.g. 3 dim grid, (xID,yID,zID)) + @ In, returnDict, bool, optional, flag to request the output in dictionary format or not. + if True a dict ( {dimName1: + coordinate1,dimName2:coordinate2,etc} is returned + if False a tuple is riturned (coordinate1,coordinate2,etc) + @ In, recastMethods, dict, optional, dictionary containing the methods that need to be used for trasforming the coordinates + ex. {'dimName1':[methodToTransformCoordinate,*args]} + @ Out, coordinate, tuple or dict, tuple (if returnDict=False) or dict (if returnDict=True) containing the coordinates + """ + if isinstance(multiDimNDIndex[0], Number): + level, multiDimIndex = self.multiGridIterator[0], multiDimNDIndex + else: + level, multiDimIndex = multiDimNDIndex[0], multiDimNDIndex[1] + node = self.grid.find(self.mappingLevelName[level]) + + return node.get('grid').returnCoordinateFromIndex(multiDimIndex, returnDict, recastMethods) + def returnParameter(self, parameterName, nodeName=None): """ Method to return one of the initialization parameters diff --git a/ravenframework/Models/PostProcessors/BasicStatistics.py b/ravenframework/Models/PostProcessors/BasicStatistics.py index 34450cca98..25f29dab6e 100644 --- a/ravenframework/Models/PostProcessors/BasicStatistics.py +++ b/ravenframework/Models/PostProcessors/BasicStatistics.py @@ -154,7 +154,7 @@ def __init__(self): self.biased = False # biased statistics? self.pivotParameter = None # time-dependent statistics pivot parameter self.pivotValue = None # time-dependent statistics pivot parameter values - self.dynamic = False # is it time-dependent? + self.dynamic = None # is it time-dependent? self.sampleTag = None # Tag used to track samples self.pbPresent = False # True if the ProbabilityWeight is available self.realizationWeight = None # The joint probabilities @@ -172,13 +172,23 @@ def inputToInternal(self, currentInp): @ Out, (inputDataset, pbWeights), tuple, the dataset of inputs and the corresponding variable probability weight """ # The BasicStatistics postprocessor only accept DataObjects - self.dynamic = False + if self.dynamic is None: + self.dynamic = False currentInput = currentInp [-1] if type(currentInp) == list else currentInp if len(currentInput) == 0: self.raiseAnError(IOError, "In post-processor " +self.name+" the input "+currentInput.name+" is empty.") pbWeights = None if type(currentInput).__name__ == 'tuple': + # if tuple we check that we already have a dataset + # and store the probability weights + if len(currentInput) != 2: + self.raiseAnError(RuntimeError, "If tuple is sent in, the dataset and the pb weights must be sent in!") + if type(currentInput[0]).__name__ != 'Dataset' or type(currentInput[1]).__name__ != 'Dataset': + self.raiseAnError(RuntimeError, "If tuple is sent in, the elements must be Dataset!") + if 'ProbabilityWeight' in currentInput[1]: + self.realizationWeight = xr.Dataset() + self.realizationWeight['ProbabilityWeight'] = currentInput[1]['ProbabilityWeight'] return currentInput # TODO: convert dict to dataset, I think this will be removed when DataSet is used by other entities that # are currently using this Basic Statisitics PostProcessor. @@ -193,8 +203,9 @@ def inputToInternal(self, currentInp): self.pbPresent = True if 'ProbabilityWeight' in metadata else False if self.pbPresent: pbWeights = xr.Dataset() + pbWeights['ProbabilityWeight'] = metadata['ProbabilityWeight']/metadata['ProbabilityWeight'].sum() self.realizationWeight = xr.Dataset() - self.realizationWeight['ProbabilityWeight'] = metadata['ProbabilityWeight']/metadata['ProbabilityWeight'].sum() + self.realizationWeight['ProbabilityWeight'] = pbWeights['ProbabilityWeight'] for target in self.parameters['targets']: pbName = 'ProbabilityWeight-' + target if pbName in metadata: @@ -246,6 +257,7 @@ def inputToInternal(self, currentInp): if self.pbPresent: pbWeights = xr.Dataset() self.realizationWeight = dataSet[['ProbabilityWeight']]/dataSet[['ProbabilityWeight']].sum() + pbWeights['ProbabilityWeight'] = self.realizationWeight['ProbabilityWeight'] for target in self.parameters['targets']: pbName = 'ProbabilityWeight-' + target if pbName in metaVars: @@ -492,7 +504,7 @@ def __computeUnbiasedCorrection(self,order,weightsOrN): denom = (v1Square-V2)*(v1Square**2.0-6.0*v1Square*V2+8.0*V1*V3+3.0*V2**2.0-6.0*V4) corrFactor = numer1/denom ,numer2/denom else: - if order == 2: + if order == 2: corrFactor = float(weightsOrN)/(float(weightsOrN)-1.0) elif order == 3: corrFactor = (float(weightsOrN)**2.0)/((float(weightsOrN)-1)*(float(weightsOrN)-2)) @@ -673,7 +685,7 @@ def _computeSingleWeightedPercentile(self, pct, weightsCDF, sortedWeightsAndPoin return result - def __runLocal(self, inputData): + def _runLocal(self, inputData): """ This method executes the postprocessor action. In this case, it computes all the requested statistical FOMs @ In, inputData, tuple, (inputDataset, pbWeights), tuple, the dataset of inputs and the corresponding @@ -1550,7 +1562,8 @@ def run(self, inputIn): @ Out, outputSet, xarray.Dataset or dictionary, dataset or dictionary containing the results """ inputData = self.inputToInternal(inputIn) - outputSet = self.__runLocal(inputData) + outputSet = self._runLocal(inputData) + print(outputSet) return outputSet def collectOutput(self, finishedJob, output): diff --git a/ravenframework/Models/PostProcessors/LimitSurface.py b/ravenframework/Models/PostProcessors/LimitSurface.py index 751aa1e619..c5232fe96b 100644 --- a/ravenframework/Models/PostProcessors/LimitSurface.py +++ b/ravenframework/Models/PostProcessors/LimitSurface.py @@ -27,7 +27,6 @@ from ...utils import InputData, InputTypes, utils, mathUtils from ...SupervisedLearning import factory as romFactory from ... import GridEntities -from ... import Files #Internal Modules End-------------------------------------------------------------------------------- class LimitSurface(PostProcessorInterface): @@ -139,18 +138,6 @@ def _initializeLSpp(self, runInfo, inputs, initDict): self.indexes = index if self.indexes == -1: self.raiseAnError(IOError, 'LimitSurface PostProcessor needs a PointSet as INPUT!!!!!!') - #else: - # # check if parameters are contained in the data - # inpKeys = self.inputs[self.indexes].getParaKeys("inputs") - # outKeys = self.inputs[self.indexes].getParaKeys("outputs") - # self.paramType = {} - # for param in self.parameters['targets']: - # if param not in inpKeys + outKeys: - # self.raiseAnError(IOError, 'LimitSurface PostProcessor: The param ' + param + ' not contained in Data ' + self.inputs[self.indexes].name + ' !') - # if param in inpKeys: - # self.paramType[param] = 'inputs' - # else: - # self.paramType[param] = 'outputs' if self.bounds == None: dataSet = self.inputs[self.indexes].asDataset() self.bounds = {"lowerBounds":{},"upperBounds":{}} diff --git a/ravenframework/Models/PostProcessors/PostProcessorInterface.py b/ravenframework/Models/PostProcessors/PostProcessorInterface.py index ed21f8cd3f..f540d2177a 100644 --- a/ravenframework/Models/PostProcessors/PostProcessorInterface.py +++ b/ravenframework/Models/PostProcessors/PostProcessorInterface.py @@ -66,7 +66,7 @@ def __init__(self): ## However, the DataObject.load can not be directly used to collect single realization ## One possible solution is all postpocessors return a list of realizations, and we only ## use addRealization method to add the collections into the DataObjects - self.outputMultipleRealizations = False + def _handleInput(self, paramInput): """ diff --git a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py index 47dd663676..d5ffc7fd85 100644 --- a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py +++ b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py @@ -60,7 +60,7 @@ class cls. descr=r"""Name of the variable for this grid/subdomain. \nb As for the other objects, this is the name that can be used to refer to this specific entity from other input blocks""") gridInput = InputData.parameterInputFactory("grid", contentType=InputTypes.StringType) - gridInput.addParam("type", InputTypes.StringType) + gridInput.addParam("type", InputTypes.makeEnumType("type", "selectionType",['value'])) gridInput.addParam("construction", InputTypes.StringType) gridInput.addParam("steps", InputTypes.IntegerType) variableInput.addSub(gridInput) @@ -75,9 +75,13 @@ def __init__(self): @ Out, None """ super().__init__() - from ...Models.PostProcessors import factory as ppFactory # delay import to allow definition + # delay import to allow definition + from ...Models.PostProcessors import factory as ppFactory + from ... import GridEntities self.stat = ppFactory.returnInstance('BasicStatistics') + self.gridEntity = GridEntities.factory.returnInstance("GridEntity") self.validDataType = ['PointSet', 'HistorySet', 'DataSet'] + self.outputMultipleRealizations = True self.printTag = 'PostProcessor SUBDOMAIN STATISTICS' def inputToInternal(self, currentInp): @@ -85,93 +89,47 @@ def inputToInternal(self, currentInp): Method to convert an input object into the internal format that is understandable by this pp. @ In, currentInp, object, an object that needs to be converted - @ Out, (inputDataset, pbWeights), tuple, the dataset of inputs and the corresponding variable probability weight + @ Out, cellBasedData, dict, the dict of datasets of inputs and the corresponding variable probability weight (cellId:(processedDataSet, pbWeights)) """ - # The BasicStatistics postprocessor only accept DataObjects + cellBasedData = {} + cellIDs = self.gridEntity.returnCellIdsWithCoordinates() + dimensionNames = self.gridEntity.returnParameter('dimensionNames') self.dynamic = False currentInput = currentInp [-1] if type(currentInp) == list else currentInp if len(currentInput) == 0: self.raiseAnError(IOError, "In post-processor " +self.name+" the input "+currentInput.name+" is empty.") - - pbWeights = None - if type(currentInput).__name__ == 'tuple': - return currentInput - # TODO: convert dict to dataset, I think this will be removed when DataSet is used by other entities that - # are currently using this Basic Statisitics PostProcessor. - if type(currentInput).__name__ == 'dict': - if 'targets' not in currentInput.keys(): - self.raiseAnError(IOError, 'Did not find targets in the input dictionary') - inputDataset = xr.Dataset() - for var, val in currentInput['targets'].items(): - inputDataset[var] = val - if 'metadata' in currentInput.keys(): - metadata = currentInput['metadata'] - self.pbPresent = True if 'ProbabilityWeight' in metadata else False - if self.pbPresent: - pbWeights = xr.Dataset() - self.realizationWeight = xr.Dataset() - self.realizationWeight['ProbabilityWeight'] = metadata['ProbabilityWeight']/metadata['ProbabilityWeight'].sum() - for target in self.parameters['targets']: - pbName = 'ProbabilityWeight-' + target - if pbName in metadata: - pbWeights[target] = metadata[pbName]/metadata[pbName].sum() - elif self.pbPresent: - pbWeights[target] = self.realizationWeight['ProbabilityWeight'] - else: - self.raiseAWarning('BasicStatistics postprocessor did not detect ProbabilityWeights! Assuming unit weights instead...') - else: - self.raiseAWarning('BasicStatistics postprocessor did not detect ProbabilityWeights! Assuming unit weights instead...') - if 'RAVEN_sample_ID' not in inputDataset.sizes.keys(): - self.raiseAWarning('BasicStatisitics postprocessor did not detect RAVEN_sample_ID! Assuming the first dimension of given data...') - self.sampleTag = utils.first(inputDataset.sizes.keys()) - return inputDataset, pbWeights - if currentInput.type not in ['PointSet','HistorySet']: - self.raiseAnError(IOError, self, 'BasicStatistics postprocessor accepts PointSet and HistorySet only! Got ' + currentInput.type) + self.raiseAnError(IOError, self, 'This Postprocessor accepts PointSet and HistorySet only! Got ' + currentInput.type) # extract all required data from input DataObjects, an input dataset is constructed dataSet = currentInput.asDataset() - try: - inputDataset = dataSet[self.parameters['targets']] - except KeyError: - missing = [var for var in self.parameters['targets'] if var not in dataSet] - self.raiseAnError(KeyError, "Variables: '{}' missing from dataset '{}'!".format(", ".join(missing),currentInput.name)) - self.sampleTag = currentInput.sampleTag - - if currentInput.type == 'HistorySet': - dims = inputDataset.sizes.keys() - if self.pivotParameter is None: - if len(dims) > 1: - self.raiseAnError(IOError, self, 'Time-dependent statistics is requested (HistorySet) but no pivotParameter \ - got inputted!') - elif self.pivotParameter not in dims: - self.raiseAnError(IOError, self, 'Pivot parameter', self.pivotParameter, 'is not the associated index for \ - requested variables', ','.join(self.parameters['targets'])) - else: - self.dynamic = True - if not currentInput.checkIndexAlignment(indexesToCheck=self.pivotParameter): - self.raiseAnError(IOError, "The data provided by the data objects", currentInput.name, "is not synchronized!") - self.pivotValue = inputDataset[self.pivotParameter].values - if self.pivotValue.size != len(inputDataset.groupby(self.pivotParameter)): - msg = "Duplicated values were identified in pivot parameter, please use the 'HistorySetSync'" + \ - " PostProcessor to syncronize your data before running 'BasicStatistics' PostProcessor." - self.raiseAnError(IOError, msg) - # extract all required meta data - metaVars = currentInput.getVars('meta') - self.pbPresent = True if 'ProbabilityWeight' in metaVars else False - if self.pbPresent: - pbWeights = xr.Dataset() - self.realizationWeight = dataSet[['ProbabilityWeight']]/dataSet[['ProbabilityWeight']].sum() - for target in self.parameters['targets']: - pbName = 'ProbabilityWeight-' + target - if pbName in metaVars: - pbWeights[target] = dataSet[pbName]/dataSet[pbName].sum() - elif self.pbPresent: - pbWeights[target] = self.realizationWeight['ProbabilityWeight'] - else: - self.raiseAWarning('BasicStatistics postprocessor did not detect ProbabilityWeights! Assuming unit weights instead...') - - return inputDataset, pbWeights + processedDataSet, pbWeights = self.stat.inputToInternal(currentInput) + for cellId, verteces in cellIDs.items(): + # create masks + maskDataset = None + # hyperBox shape(number dimensions, number of verteces) + hyperBox = np.atleast_2d(verteces).T + for dim, dimName in enumerate(dimensionNames): + if maskDataset is None: + maskDataset = (dataSet[dimName] >= min(hyperBox[dim])) & (dataSet[dimName] <= max(hyperBox[dim])) + else: + maskDataset = maskDataset & (dataSet[dimName] >= min(hyperBox[dim])) & (dataSet[dimName] <= max(hyperBox[dim])) + + # the following is the cropped dataset that we need to use for the subdomain + #cellDataset = dataSet.where(maskDataset, drop=True) + cellDataset = processedDataSet.where(maskDataset, drop=True) + cellPbWeights = pbWeights.where(maskDataset, drop=True) + # check if at least sample is available (for scalar quantities) and at least 2 samples for derivative quantities + setWhat = set(self.stat.what) + minimumNumberOfSamples = 2 if len(setWhat.intersection(set(self.stat.vectorVals))) > 0 else 1 + if len(cellDataset[currentInput.sampleTag]) < minimumNumberOfSamples: + self.raiseAnError(RuntimeError,"Number of samples in cell " + f"{cellId} < {minimumNumberOfSamples}. Found {len(cellDataset[currentInput.sampleTag])}" + " samples within the cell. Please make the evaluation grid coarser or increase number of samples!") + + # store datasets + cellBasedData[cellId] = cellDataset, cellPbWeights + return cellBasedData def initialize(self, runInfo, inputs, initDict): """ @@ -182,108 +140,62 @@ def initialize(self, runInfo, inputs, initDict): @ In, initDict, dict, dictionary with initialization options @ Out, None """ - self.stat.intialize(runInfo, inputs, initDict) + self.stat.initialize(runInfo, inputs, initDict) + self.gridEntity.initialize({'computeCells': True,'constructTensor': True,}) + # check that the grid is defined only in the input space. + dims = self.gridEntity.returnParameter("dimensionNames") + inputs = inputs[-1].getVars("input") + if not all(item in inputs for item in dims): + unset = ', '.join(list(set(dims) - set(inputs))) + self.raiseAnError(RuntimeError, "Subdomain grid must be defined on the input space only (inputs)." + f"The following variables '{unset}' are not part of the input space of DataObject {inputs[-1].name}!") - def _handleInput(self, paramInput, childVals=None): + def _handleInput(self, paramInput): """ Function to handle the parsed paramInput for this class. @ In, paramInput, ParameterInput, the already parsed input. - @ In, childVals, list, optional, quantities requested from child statistical object @ Out, None """ # initialize basic stats - subdomain = paramInput.popSub('subdomain') + subdomain = paramInput.findFirst('subdomain') if subdomain is None: self.raiseAnError(IOError,' tag not found!') - self.stat._handleInput(paramInput, childVals) - - for child in subdomain.subparts: - if child.tag == 'variable': - varName = child.parameterValues['name'] - for childChild in child.subparts: - - - - - # variable for tracking if distributions or functions have been declared - foundDistOrFunc = False - # store variable name for re-use - varName = child.parameterValues['name'] - # set shape if present - if 'shape' in child.parameterValues: - self.variableShapes[varName] = child.parameterValues['shape'] - # read subnodes - for childChild in child.subparts: - if childChild.getName() == 'distribution': - # can only have a distribution if doesn't already have a distribution or function - if foundDistOrFunc: - self.raiseAnError(IOError, 'A sampled variable cannot have both a distribution and a function, or more than one of either!') - else: - foundDistOrFunc = True - # name of the distribution to sample - toBeSampled = childChild.value - varData = {} - varData['name'] = childChild.value - # variable dimensionality - if 'dim' not in childChild.parameterValues: - dim = 1 - else: - dim = childChild.parameterValues['dim'] - varData['dim'] = dim - # set up mapping for variable to distribution - self.variables2distributionsMapping[varName] = varData - # flag distribution as needing to be sampled - self.toBeSampled[prefix + varName] = toBeSampled - elif childChild.getName() == 'function': - # can only have a function if doesn't already have a distribution or function - if not foundDistOrFunc: - foundDistOrFunc = True - else: - self.raiseAnError(IOError, 'A sampled variable cannot have both a distribution and a function!') - # function name - toBeSampled = childChild.value - # track variable as a functional sample - self.dependentSample[prefix + varName] = toBeSampled - - - else: - if tag not in childVals: - self.raiseAWarning('Unrecognized node in BasicStatistics "',tag,'" has been ignored!') - - - - - - - - - - def __runLocal(self, inputData): - """ - This method executes the postprocessor action. In this case, it computes all the requested statistical FOMs - @ In, inputData, tuple, (inputDataset, pbWeights), tuple, the dataset of inputs and the corresponding - variable probability weight - @ Out, outputSet or outputDict, xarray.Dataset or dict, dataset or dictionary containing the results - """ - inputDataset, pbWeights = inputData[0], inputData[1] - - - return outputDict - - - - + self.stat._handleInput(paramInput, ['subdomain']) + self.gridEntity._handleInput(subdomain, dimensionTags=["variable"]) def run(self, inputIn): """ This method executes the postprocessor action. In this case, it computes all the requested statistical FOMs @ In, inputIn, object, object contained the data to process. (inputToInternal output) - @ Out, outputSet, xarray.Dataset or dictionary, dataset or dictionary containing the results + @ Out, results, xarray.Dataset or dictionary, dataset or dictionary containing the results """ inputData = self.inputToInternal(inputIn) - outputSet = self.__runLocal(inputData) - return outputSet + results = {} + outputRealization = {} + midPoint = self.gridEntity.returnCellsMidPoints(returnDict=True) + firstPass = True + for i, (cellId, data) in enumerate(inputData.items()): + cellData = self.stat.inputToInternal(data) + res = self.stat._runLocal(cellData) + for k in res: + if firstPass: + results[k] = np.zeros(len(inputData), dtype=object if self.stat.dynamic else None) + results[k][i] = res[k][0] if not self.stat.dynamic else res[k] + for k in midPoint[cellId]: + #res[k] = np.atleast_1d(midPoint[cellId][k]) + if firstPass: + results[k] = np.zeros(len(inputData)) + results[k][i] = np.atleast_1d(midPoint[cellId][k]) + firstPass = False + outputRealization['data'] = results + if self.stat.dynamic: + dims = dict.fromkeys(results.keys(), inputIn[-1].indexes if type(inputIn) == list else inputIn.indexes) + for k in list(midPoint.values())[0]: + dims[k] = [] + outputRealization['dims'] = dims + + return outputRealization def collectOutput(self, finishedJob, output): """ diff --git a/ravenframework/Samplers/Grid.py b/ravenframework/Samplers/Grid.py index 88327bbca2..8aee4f482e 100644 --- a/ravenframework/Samplers/Grid.py +++ b/ravenframework/Samplers/Grid.py @@ -99,7 +99,6 @@ def localInputAndChecks(self, xmlNode, paramInput): if 'limit' in paramInput.parameterValues: self.raiseAnError(IOError,'limit is not used in Grid sampler') self.limit = 1 - # FIXME: THIS READ MORE XML MUST BE CONVERTED IN THE INPUTPARAMETER COLLECTOR!!!!!!! self.gridEntity._handleInput(paramInput, dimensionTags=["variable", "Distribution"], dimTagsPrefix={"Distribution": ""}) grdInfo = self.gridEntity.returnParameter("gridInfo") diff --git a/ravenframework/Samplers/LimitSurfaceSearch.py b/ravenframework/Samplers/LimitSurfaceSearch.py index 564cd061a8..d973bd00d4 100644 --- a/ravenframework/Samplers/LimitSurfaceSearch.py +++ b/ravenframework/Samplers/LimitSurfaceSearch.py @@ -165,7 +165,7 @@ def __init__(self): # cutoff (% of range space) self.sizeGrid = None # size of grid self.sizeSubGrid = None # size of subgrid - self.printTag = 'SAMPLER ADAPTIVE' + self.printTag = 'ADAPTIVE LIMIT SURFACE' self.acceptedScoringParam = ['distance','distancePersistence'] self.acceptedBatchParam = ['none','naive','maxV','maxP'] diff --git a/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainStats/subdomainSensitivity.csv b/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainStats/subdomainSensitivity.csv new file mode 100644 index 0000000000..a1e3b83234 --- /dev/null +++ b/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainStats/subdomainSensitivity.csv @@ -0,0 +1,5 @@ +skew_ans,skew_x1,skew_x2,skew_x3,skew_x4,skew_x5,vc_ans,vc_x1,vc_x2,vc_x3,vc_x4,vc_x5,percentile_5_ans,percentile_95_ans,percentile_5_x1,percentile_95_x1,percentile_5_x2,percentile_95_x2,percentile_5_x3,percentile_95_x3,percentile_5_x4,percentile_95_x4,percentile_5_x5,percentile_95_x5,mean_ans,mean_x1,mean_x2,mean_x3,mean_x4,mean_x5,kurt_ans,kurt_x1,kurt_x2,kurt_x3,kurt_x4,kurt_x5,median_ans,median_x1,median_x2,median_x3,median_x4,median_x5,max_ans,max_x1,max_x2,max_x3,max_x4,max_x5,min_ans,min_x1,min_x2,min_x3,min_x4,min_x5,samp_ans,samp_x1,samp_x2,samp_x3,samp_x4,samp_x5,var_ans,var_x1,var_x2,var_x3,var_x4,var_x5,sigma_ans,sigma_x1,sigma_x2,sigma_x3,sigma_x4,sigma_x5,nsen_ans_x1,nsen_ans_x2,nsen_ans_x3,nsen_ans_x4,nsen_ans_x5,nsen_x1_x1,nsen_x1_x2,nsen_x1_x3,nsen_x1_x4,nsen_x1_x5,nsen_x2_x1,nsen_x2_x2,nsen_x2_x3,nsen_x2_x4,nsen_x2_x5,nsen_x3_x1,nsen_x3_x2,nsen_x3_x3,nsen_x3_x4,nsen_x3_x5,nsen_x4_x1,nsen_x4_x2,nsen_x4_x3,nsen_x4_x4,nsen_x4_x5,nsen_x5_x1,nsen_x5_x2,nsen_x5_x3,nsen_x5_x4,nsen_x5_x5,sen_ans_x1,sen_ans_x2,sen_ans_x3,sen_ans_x4,sen_ans_x5,sen_x1_x1,sen_x1_x2,sen_x1_x3,sen_x1_x4,sen_x1_x5,sen_x2_x1,sen_x2_x2,sen_x2_x3,sen_x2_x4,sen_x2_x5,sen_x3_x1,sen_x3_x2,sen_x3_x3,sen_x3_x4,sen_x3_x5,sen_x4_x1,sen_x4_x2,sen_x4_x3,sen_x4_x4,sen_x4_x5,sen_x5_x1,sen_x5_x2,sen_x5_x3,sen_x5_x4,sen_x5_x5,pear_ans_x1,pear_ans_x2,pear_ans_x3,pear_ans_x4,pear_ans_x5,pear_x1_x1,pear_x1_x2,pear_x1_x3,pear_x1_x4,pear_x1_x5,pear_x2_x1,pear_x2_x2,pear_x2_x3,pear_x2_x4,pear_x2_x5,pear_x3_x1,pear_x3_x2,pear_x3_x3,pear_x3_x4,pear_x3_x5,pear_x4_x1,pear_x4_x2,pear_x4_x3,pear_x4_x4,pear_x4_x5,pear_x5_x1,pear_x5_x2,pear_x5_x3,pear_x5_x4,pear_x5_x5,cov_ans_x1,cov_ans_x2,cov_ans_x3,cov_ans_x4,cov_ans_x5,cov_x1_x1,cov_x1_x2,cov_x1_x3,cov_x1_x4,cov_x1_x5,cov_x2_x1,cov_x2_x2,cov_x2_x3,cov_x2_x4,cov_x2_x5,cov_x3_x1,cov_x3_x2,cov_x3_x3,cov_x3_x4,cov_x3_x5,cov_x4_x1,cov_x4_x2,cov_x4_x3,cov_x4_x4,cov_x4_x5,cov_x5_x1,cov_x5_x2,cov_x5_x3,cov_x5_x4,cov_x5_x5,vsen_ans_x1,vsen_ans_x2,vsen_ans_x3,vsen_ans_x4,vsen_ans_x5,vsen_x1_x1,vsen_x1_x2,vsen_x1_x3,vsen_x1_x4,vsen_x1_x5,vsen_x2_x1,vsen_x2_x2,vsen_x2_x3,vsen_x2_x4,vsen_x2_x5,vsen_x3_x1,vsen_x3_x2,vsen_x3_x3,vsen_x3_x4,vsen_x3_x5,vsen_x4_x1,vsen_x4_x2,vsen_x4_x3,vsen_x4_x4,vsen_x4_x5,vsen_x5_x1,vsen_x5_x2,vsen_x5_x3,vsen_x5_x4,vsen_x5_x5 +1.13562319798,1.18784679221,1.36862292907,1.68026867554,-0.856171580956,1.60736514212,6.96592176762,7.093032034,36.2916928212,16.7881473046,-89.6041251252,9.79144006142,693.369022927,713.216314561,0.295101283545,0.488279878488,-0.346539657264,0.820699832965,-1.51534435979,4.30601148378,-1.30344935706,0.91440289273,-0.0161693129385,0.216415508654,91.6234563368,0.0532264704085,0.0119515887199,0.122026631735,-0.00993950137461,0.012603305062,-2.35336080934,-2.2102337005,-0.273393970926,0.585576206834,-0.408970726581,-0.858084366364,705.987336016,0.42610796341,0.0482174514157,0.829831545839,0.0667250760007,0.0843267249239,714.817120398,0.494519153268,0.869971266171,4.31436874693,1.11162271183,0.231029931511,692.050662768,0.244892745738,-0.410541436108,-2.06891738295,-1.94142441906,-0.0521791665876,13.0,13.0,13.0,13.0,13.0,13.0,407352.632184,0.14253423142,0.188133325375,4.19676633852,0.793204563041,0.0152286721233,638.241828921,0.377537059665,0.433743386549,2.04860106866,0.890620324853,0.12340450609,0.114401640288,0.00441500308021,0.0186625706652,0.000441056131497,0.0579705991113,1.0,0.00479397271964,0.0195422499145,0.000262477699014,0.064507593418,0.128528508135,1.0,0.0629935426656,-0.0044066933008,-0.133662198003,0.11176065565,0.0134371295055,1.0,0.00489475504665,-0.00163823289325,0.0429708645621,-0.0269085687495,0.140119399135,1.0,0.253193786567,0.123987329986,-0.00958233077005,-0.000550589254462,0.00297260928573,1.0,1.0,2.2,-3.3,-0.2,11.0,1.0,0.0863187975092,0.0148720158683,0.0276083059406,0.833332940759,3.59887257493,1.0,-0.0599289752292,-0.114729191511,-4.94457649092,15.0282197079,-1.45249102893,1.0,-0.899281834097,-20.2620236215,7.12615942877,-0.710278420555,-0.229706825709,1.0,-9.14894619762,0.701125288236,-0.0997803662297,-0.0168702766301,-0.0298216988657,1.0,0.127670324275,0.0248095085664,0.0527687436795,-0.009930386082,0.0928742058462,1.0,0.0260183032322,0.0539431374313,-0.00803011562001,0.0995389320182,0.0260183032322,1.0,0.0311637456934,0.0107183135863,-0.0314646527131,0.0539431374313,0.0311637456934,1.0,-0.0266203250539,0.00784897924831,-0.00803011562001,0.0107183135863,-0.0266203250539,1.0,-0.02909491666,0.0995389320182,-0.0314646527131,0.00784897924831,-0.02909491666,1.0,30.7634341169,6.86809495997,68.9952850099,-5.64474073097,7.31495055396,0.14253423142,0.00426060650418,0.0417208536815,-0.00270006357147,0.00463749638461,0.00426060650418,0.188133325375,0.0276910811681,0.00414049178216,-0.00168417348925,0.0417208536815,0.0276910811681,4.19676633852,-0.0485694685147,0.0019842737812,-0.00270006357147,0.00414049178216,-0.0485694685147,0.793204563041,-0.00319772224147,0.00463749638461,-0.00168417348925,0.0019842737812,-0.00319772224147,0.0152286721233,196.929715859,33.8463656527,14.0127544631,-4.06570568113,421.434427746,1.0,0.0213499856028,0.00852408176804,-0.00140557971199,0.272429453648,0.0288600738735,1.0,0.00616974264753,0.00529875533602,-0.126750531712,0.25622169317,0.137194116378,1.0,-0.0600925991155,-0.0158615570264,-0.00802437140966,0.0223784270311,-0.0114132213641,1.0,-0.199679368012,0.0293585152583,-0.0101048522275,-5.68666383652e-05,-0.00376927375388,1.0 +1.24756816872,1.29212127903,1.37676123697,1.62139352565,-0.161134865941,1.51243092675,12.2951668821,12.4423399571,12.7280819462,31.1391034028,266.103676643,13.840841976,695.061931519,717.699218183,0.330199792239,0.485685780983,1.09365770285,2.11839899768,-1.38242457465,4.58106009969,-1.44425462477,1.32599293042,0.0337981183776,0.201293493705,56.6184758955,0.0319781673687,0.124661265927,0.090005856221,0.00414342348461,0.0103719048427,-2.64889119014,-2.50521265606,-2.21944366318,-0.766490093319,-1.28599606425,-1.85486520665,708.902866212,0.389498565698,1.6181753738,0.704148310794,0.0129682204659,0.144945299892,718.639090869,0.4958700578,2.3149614891,5.0811880962,1.38397023914,0.202365835728,694.058890028,0.317489474131,1.07041766172,-1.41316396069,-1.77978546462,0.0154569346876,8.0,8.0,8.0,8.0,8.0,8.0,484602.002618,0.158311064401,2.51761310804,7.85513661582,1.21568314847,0.0206082952528,696.133609746,0.397883229606,1.58669880823,2.80270166372,1.10258022314,0.143555895918,0.0681262705468,0.0647931091002,0.0104018158818,0.000115225301739,0.0549505381677,1.0,0.0645644829675,0.0104875649322,-1.8086498843e-06,0.0521963297537,0.0675677175287,1.0,0.00841685034141,-0.000222369638768,0.0568761432104,0.066274370612,0.0508247186943,1.0,0.00766125645329,0.0033697760906,-0.00083610022443,-0.0982276501308,0.560444854073,1.0,-0.386274880076,0.0647216631825,0.0673898160638,0.000661210586243,-0.0010361015757,1.0,1.0,2.2,-3.3,-0.2,11.0,1.0,-0.0209983261317,0.00379100566348,-0.0579451359603,-0.863726116848,-2.2259714188,1.0,0.032512340435,-0.290792490734,-1.31862364075,4.13815296626,0.334784570168,1.0,1.98229966457,-4.30599620273,-7.16621613145,-0.339251046675,0.224589845554,1.0,-6.8451341549,-0.666116638458,-0.00959312670238,-0.00304225606414,-0.0426857855373,1.0,0.0790469114661,0.0773774463594,0.0305831092996,0.00272639663796,0.0714118445158,1.0,0.076174482134,0.0302958511892,0.000426211100998,0.0678423202522,0.076174482134,1.0,0.0246737020315,-0.00432105210789,0.0713260952514,0.0302958511892,0.0246737020315,1.0,0.0654345974481,0.00537283835087,0.000426211100998,-0.00432105210789,0.0654345974481,1.0,-0.0198966002795,0.0678423202522,0.0713260952514,0.00537283835087,-0.0198966002795,1.0,21.8944047544,85.4675964348,59.6693229995,2.09262706575,7.13647727027,0.158311064401,0.048090538522,0.0337843575463,0.000186978034339,0.00387505044971,0.048090538522,2.51761310804,0.109725023517,-0.00755951960004,0.0162466557629,0.0337843575463,0.109725023517,7.85513661582,0.202206217199,0.00216173114494,0.000186978034339,-0.00755951960004,0.202206217199,1.21568314847,-0.00314927153173,0.00387505044971,0.0162466557629,0.00216173114494,-0.00314927153173,0.0206082952528,120.619970567,29.4276418461,6.54329603096,1.57451464792,299.965702335,1.0,0.0165621119524,0.00372612539643,-1.39588214737e-05,0.160929259785,0.263400873033,1.0,0.0116576327667,-0.00669033246914,0.683601722265,0.186536063924,0.0366956190382,1.0,0.166422271189,0.0292424185245,-0.000108333828685,-0.00326483731225,0.025800103101,1.0,-0.154311135119,0.0209920388511,0.00560688000707,7.61951896181e-05,-0.00259359126347,1.0 +1.11534307975,1.14697427582,1.6950498224,-1.21589760574,-0.777381466638,1.45331216359,5.8666326522,5.92039507558,9.8306075667,133.544093958,-15.4208189086,7.22868265057,697.09266223,728.620532506,0.527290348308,0.73395145956,-0.12689950564,0.781862090073,-5.42236270918,3.52980793877,-2.26064434654,1.52140190876,0.00875322644744,0.266499523255,106.305165571,0.0911985881881,0.04441982754,0.0247336423163,-0.0964378575282,0.0225978918057,-2.29822477352,-2.20596164537,-0.221769155426,1.60381083931,-0.854376051465,-1.36226527702,706.98529778,0.582637539753,0.268948453781,0.872878778555,-0.808375943761,0.195997397191,740.371277705,0.736948595145,0.991436733561,4.8202101499,2.25601859132,0.293632167942,690.376182739,0.502889102412,-0.218841823174,-8.645075201,-2.97940813934,-0.0155674690883,15.0,15.0,15.0,15.0,15.0,15.0,388943.507743,0.29152621087,0.190684088588,10.9100194247,2.21161731418,0.0266841988308,623.653355433,0.539931672409,0.436673892726,3.30303185342,1.48715073687,0.163352988435,0.127076356387,0.0419840580939,0.000514759228551,0.0175102059453,0.0812668569294,1.0,0.0406886908529,0.000726840025153,0.0128880580059,0.0794015455834,0.114104745985,1.0,-0.0040245417105,0.00917489816798,0.0906972613251,0.382443586758,-0.755118288795,1.0,-0.0976254180474,-0.789770104893,0.0903846081628,0.0229444820598,-0.00130119191554,1.0,0.0556366253012,0.119073231908,0.0485008519037,-0.00225090583928,0.0118970412461,1.0,1.0,2.2,-3.3,-0.2,11.0,1.0,-0.0390046382609,-0.00537911645997,0.0334207764657,-0.22226801547,-0.967741517234,1.0,-0.0689338383968,0.0738017076545,-0.826089562838,-7.9654773274,-4.11424549039,1.0,0.720549135868,-18.4905544438,10.9094204222,0.970975865385,0.15883561392,1.0,4.95015044063,-0.350998168771,-0.0525789483301,-0.0197186397977,0.0239475524725,1.0,0.148544127703,0.0907315712525,0.00372735277313,-0.0573255870123,0.122726758602,1.0,0.0875263269889,0.008281920253,-0.0460775324809,0.119117243031,0.0875263269889,1.0,-0.0569858744635,-0.0248927698036,0.0866556143507,0.008281920253,-0.0569858744635,1.0,0.0127968369262,-0.0444831818072,-0.0460775324809,-0.0248927698036,0.0127968369262,1.0,-0.0378347345827,0.119117243031,0.0866556143507,-0.0444831818072,-0.0378347345827,1.0,50.0192937116,24.7092135537,7.67814878458,-53.167564246,12.5028669989,0.29152621087,0.0206364379259,0.014770071926,-0.0369984055505,0.0105060755875,0.0206364379259,0.190684088588,-0.0821934493661,-0.0161653622546,0.00618131701322,0.014770071926,-0.0821934493661,10.9100194247,0.0628594227061,-0.0240013510987,-0.0369984055505,-0.0161653622546,0.0628594227061,2.21161731418,-0.00919121163736,0.0105060755875,0.00618131701322,-0.0240013510987,-0.00919121163736,0.0266841988308,148.125901664,100.475902185,2.21243455857,-19.3018114452,382.296134327,1.0,0.0835381712742,0.00268002517724,-0.0121878764704,0.320441788084,0.055576662302,1.0,-0.00722778499103,-0.00422601045652,0.17828020158,0.103721154779,-0.420461462728,1.0,0.0250382187334,-0.864412107753,-0.0955771151497,-0.0498137164074,0.00507342020084,1.0,-0.237432632667,0.0295048866992,0.0246740490566,-0.00205654007486,-0.00278778539651,1.0 +1.29675710115,1.31985540071,1.3994321439,1.55601985377,-0.444733446486,1.88600496,14.3525018896,14.4392762982,14.7882661314,-126.821415042,-129.33074584,21.372913521,692.541022029,721.254658924,0.528704721559,0.686620796382,1.14088482806,2.21836591872,-2.93918191421,5.57468244736,-1.26326299676,1.02284890766,-0.0223444368361,0.19105829425,49.9025748049,0.0426917728215,0.120313560552,-0.0299026943916,-0.00723785576764,0.00555152147673,-2.78369259604,-2.70640307532,-2.44660309341,0.917793514346,-1.57236044744,-0.657143549239,719.585695509,0.593073837655,1.85344540041,-2.35355098539,-0.299147291481,0.0948849803987,721.297404908,0.690784072691,2.31446826745,7.54417672423,1.19981248845,0.196922317121,685.972528649,0.521843100719,1.03016364958,-3.1145916485,-1.32923968002,-0.0364431326487,7.0,7.0,7.0,7.0,7.0,7.0,512980.827866,0.379996181936,3.16565566602,14.3815545829,0.876240682958,0.0140783418197,716.226799182,0.61643830343,1.77922895267,3.79230201631,0.936077284714,0.118652188432,0.0628801570384,0.0600826076937,0.000860340211323,0.000753928960821,0.0281728184018,1.0,0.060828342531,0.000238497797033,0.00105540279269,0.0271444421025,0.0638404511841,1.0,0.00116447745906,0.00178207361575,0.0219145629188,0.0185865425377,0.086467989174,1.0,0.0105775517202,0.0704967709646,0.0855124161755,0.137577246839,0.010997193283,1.0,-0.0923335705768,0.0597860793936,0.045989851027,0.00199238902775,-0.00250996604294,1.0,1.0,2.2,-3.3,-0.2,11.0,1.0,0.122278927312,0.0140052132051,-0.0171675930145,0.366754946272,3.49191718197,1.0,-0.066198461096,-0.0112089571678,-2.95710251791,42.7667332961,-7.07867433703,1.0,0.822298217691,-24.4139575434,-6.34039433333,-0.144963863461,0.0994534504795,1.0,2.89279381031,0.677366316434,-0.191250105416,-0.0147662238604,0.0144663256206,1.0,0.069492553565,0.0679902351299,-0.00915949891356,-0.00797465123317,0.0471373072574,1.0,0.0682790009919,-0.00399026801289,-0.0104486434157,0.0453831087256,0.0682790009919,1.0,-0.011369949677,-0.0164026012002,0.0373320696248,-0.00399026801289,-0.011369949677,1.0,0.0108718825907,-0.0125665238357,-0.0104486434157,-0.0164026012002,0.0108718825907,1.0,0.0137793082735,0.0453831087256,0.0373320696248,-0.0125665238357,0.0137793082735,1.0,30.6816318178,86.642095448,-24.8785577205,-5.3465541799,4.0058168539,0.379996181936,0.0748873757045,-0.00932813022106,-0.00602922138954,0.00331940016215,0.0748873757045,3.16565566602,-0.0767172976553,-0.0273184635212,0.0078811511536,-0.00932813022106,-0.0767172976553,14.3815545829,0.0385939630917,-0.00565449506122,-0.00602922138954,-0.0273184635212,0.0385939630917,0.876240682958,0.00153043495276,0.00331940016215,0.0078811511536,-0.00565449506122,0.00153043495276,0.0140783418197,73.5008535128,24.9205227669,-1.43576331922,-5.19808594877,253.245201996,1.0,0.0215841819369,-0.00034050088049,-0.0062251884684,0.208743559845,0.179914570925,1.0,-0.00468527776955,-0.029623085726,0.474936304176,-0.0130186137649,-0.0214907267564,1.0,0.0437004144121,-0.379723541841,-0.0144975599212,-0.00827640928384,0.00266183701672,1.0,0.120380884618,0.00777441839088,0.00212206873868,-0.000369892770621,0.00192517381398,1.0 diff --git a/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump.csv b/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump.csv new file mode 100644 index 0000000000..f946e050c4 --- /dev/null +++ b/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump.csv @@ -0,0 +1,3 @@ +x0,y0,filename +32.5,30.0,timedepSubdomainStats_dump_0.csv +77.5,30.0,timedepSubdomainStats_dump_1.csv diff --git a/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump_0.csv b/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump_0.csv new file mode 100644 index 0000000000..0b6d9f73b3 --- /dev/null +++ b/tests/framework/PostProcessors/SubdomainBasicStatistics/gold/subdomainTimeDepStats/timedepSubdomainStats_dump_0.csv @@ -0,0 +1,4 @@ +time,var_x0,var_y0,var_z0,var_x,var_y,var_z,mean_x0,mean_y0,mean_z0,mean_x,mean_y,mean_z,sigma_x0,sigma_y0,sigma_z0,sigma_x,sigma_y,sigma_z,vc_x0,vc_y0,vc_z0,vc_x,vc_y,vc_z,skew_x0,skew_y0,skew_z0,skew_x,skew_y,skew_z,kurt_x0,kurt_y0,kurt_z0,kurt_x,kurt_y,kurt_z,median_x0,median_y0,median_z0,median_x,median_y,median_z,samp_x0,samp_y0,samp_z0,samp_x,samp_y,samp_z,percentile_5_x0,percentile_5_y0,percentile_5_z0,percentile_5_x,percentile_5_y,percentile_5_z,percentile_95_x0,percentile_95_y0,percentile_95_z0,percentile_95_x,percentile_95_y,percentile_95_z,cov_x0_x0,cov_y0_x0,cov_z0_x0,cov_x_x0,cov_y_x0,cov_z_x0,cov_x0_y0,cov_y0_y0,cov_z0_y0,cov_x_y0,cov_y_y0,cov_z_y0,cov_x0_z0,cov_y0_z0,cov_z0_z0,cov_x_z0,cov_y_z0,cov_z_z0,cov_x0_x,cov_y0_x,cov_z0_x,cov_x_x,cov_y_x,cov_z_x,cov_x0_y,cov_y0_y,cov_z0_y,cov_x_y,cov_y_y,cov_z_y,cov_x0_z,cov_y0_z,cov_z0_z,cov_x_z,cov_y_z,cov_z_z,pearson_x0_x0,pearson_y0_x0,pearson_z0_x0,pearson_x_x0,pearson_y_x0,pearson_z_x0,pearson_x0_y0,pearson_y0_y0,pearson_z0_y0,pearson_x_y0,pearson_y_y0,pearson_z_y0,pearson_x0_z0,pearson_y0_z0,pearson_z0_z0,pearson_x_z0,pearson_y_z0,pearson_z_z0,pearson_x0_x,pearson_y0_x,pearson_z0_x,pearson_x_x,pearson_y_x,pearson_z_x,pearson_x0_y,pearson_y0_y,pearson_z0_y,pearson_x_y,pearson_y_y,pearson_z_y,pearson_x0_z,pearson_y0_z,pearson_z0_z,pearson_x_z,pearson_y_z,pearson_z_z +-0.05,848.0,212.0,53.0,891.92,234.32,64.52,12.0,6.0,3.0,12.4,6.4,3.4,29.1204395571,14.5602197786,7.28010988928,29.865029717,15.3075144945,8.03243425121,2.42670329643,2.42670329643,2.42670329643,2.40847013847,2.39179913977,2.36248066212,-4.19520360525e+15,-4.19520360525e+15,-4.19520360525e+15,-4.16025013377e+15,-4.12647562734e+15,-4.06264095907e+15,nan,nan,nan,nan,nan,nan,30.0,15.0,7.5,31.0,16.0,8.5,2.0,2.0,2.0,2.0,2.0,2.0,21.0,10.5,5.25,22.0,11.5,6.25,39.0,19.5,9.75,40.0,20.5,10.75,848.0,169.6,84.8,347.84,178.24,93.44,169.6,212.0,42.4,173.92,89.12,46.72,84.8,42.4,53.0,86.96,44.56,23.36,347.84,173.92,86.96,891.92,182.848,95.888,178.24,89.12,44.56,182.848,234.32,49.168,93.44,46.72,23.36,95.888,49.168,64.52,1.0,0.4,0.4,0.399961920485,0.399855034058,0.399473269683,0.4,1.0,0.4,0.399961920485,0.399855034058,0.399473269683,0.4,0.4,1.0,0.399961920485,0.399855034058,0.399473269683,0.399961920485,0.399961920485,0.399961920485,1.0,0.399965547881,0.399718390269,0.399855034058,0.399855034058,0.399855034058,0.399965547881,1.0,0.399880921478,0.399473269683,0.399473269683,0.399473269683,0.399718390269,0.399880921478,1.0 +0.0,848.0,212.0,53.0,883.2452,338.9117,201.8925,12.0,6.0,3.0,12.34,7.49,5.55,29.1204395571,14.5602197786,7.28010988928,29.7194414483,18.4095545845,14.2088880635,2.42670329643,2.42670329643,2.42670329643,2.40838261332,2.45788445721,2.56016001143,-4.19520360525e+15,-4.19520360525e+15,-4.19520360525e+15,-4.16007741029e+15,-4.25042411074e+15,-4.39581667348e+15,nan,nan,nan,nan,nan,nan,30.0,15.0,7.5,30.85,18.725,13.875,2.0,2.0,2.0,2.0,2.0,2.0,21.0,10.5,5.25,21.895,12.8075,8.8125,39.0,19.5,9.75,39.805,24.6425,18.9375,848.0,169.6,84.8,346.144,214.384,164.88,169.6,212.0,42.4,173.072,107.192,82.44,84.8,42.4,53.0,86.536,53.596,41.22,346.144,173.072,86.536,883.2452,218.70488,168.0516,214.384,107.192,53.596,218.70488,338.9117,104.4126,164.88,82.44,41.22,168.0516,104.4126,201.8925,1.0,0.4,0.4,0.399961546469,0.399899777672,0.39848314806,0.4,1.0,0.4,0.399961546469,0.399899777672,0.39848314806,0.4,0.4,1.0,0.399961546469,0.399899777672,0.39848314806,0.399961546469,0.399961546469,0.399961546469,1.0,0.399737184924,0.397962284983,0.399899777672,0.399899777672,0.399899777672,0.399737184924,1.0,0.399162318442,0.39848314806,0.39848314806,0.39848314806,0.397962284983,0.399162318442,1.0 +0.05,848.0,212.0,53.0,876.52840325,578.784090983,471.826916302,12.0,6.0,3.0,12.2915,9.35115,8.06685,29.1204395571,14.5602197786,7.28010988928,29.6062223739,24.0579319764,21.7215772057,2.42670329643,2.42670329643,2.42670329643,2.40867448024,2.57272442175,2.69269630719,-4.19520360525e+15,-4.19520360525e+15,-4.19520360525e+15,-4.16065319862e+15,-4.41032236144e+15,-4.51816126516e+15,nan,nan,nan,nan,nan,nan,30.0,15.0,7.5,30.72875,23.377875,20.167125,2.0,2.0,2.0,2.0,2.0,2.0,21.0,10.5,5.25,21.804125,14.7152625,11.6528775,39.0,19.5,9.75,39.653375,32.0404875,28.6813725,848.0,169.6,84.8,344.8264,278.98584,249.92616,169.6,212.0,42.4,172.4132,139.49292,124.96308,84.8,42.4,53.0,86.2066,69.74646,62.48154,344.8264,172.4132,86.2066,876.52840325,283.24750503,253.52499777,278.98584,139.49292,69.74646,283.24750503,578.784090983,208.626847407,249.92616,124.96308,62.48154,253.52499777,208.626847407,471.826916302,1.0,0.4,0.4,0.399962786288,0.398222635612,0.39511401449,0.4,1.0,0.4,0.399962786288,0.398222635612,0.39511401449,0.4,0.4,1.0,0.399962786288,0.398222635612,0.39511401449,0.399962786288,0.399962786288,0.399962786288,1.0,0.397671807893,0.394227062996,0.398222635612,0.398222635612,0.398222635612,0.397671807893,1.0,0.399227589572,0.39511401449,0.39511401449,0.39511401449,0.394227062996,0.399227589572,1.0 diff --git a/tests/framework/PostProcessors/SubdomainBasicStatistics/sensitivity_subdomain.xml b/tests/framework/PostProcessors/SubdomainBasicStatistics/sensitivity_subdomain.xml new file mode 100644 index 0000000000..6cb4862659 --- /dev/null +++ b/tests/framework/PostProcessors/SubdomainBasicStatistics/sensitivity_subdomain.xml @@ -0,0 +1,381 @@ + + + + subdomainStats + 1 + sample,PP + + + + framework/PostProcessors/SubdomainBasicStatistics.subdomainSensitivity + aalfonsi + 2023-05-17 + PostProcessors.SubdomainBasicStatistics + + This test checks the statistics and sensitivities (and other metrics) calculated + by subdomain basic statistics PP on + static data (not time-dependent). This test shows how to compute subdomain statistics + using a 2D grid defined in the input space. + + + added test + + + + + + x1,x2,x3,x4,x5 + OutputPlaceHolder + + + x1,x2,x3,x4,x5 + ans + + + InputOutput_vars + + + + + + 0.5 + 0.1 + + + -0.4 + 1.8 + + + 0.3 + 3 + + + -0.2 + 1.0 + + + 0.1 + 0.1 + + + + + + + 100 + 1234 + True + + + NDist1 + + + NDist2 + + + NDist3 + + + NDist4 + + + NDist5 + + + + + + + x1,x2,x3,x4,x5,ans + + + + + 0.0013 0.998 + + + -0.5 2.5 + + + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + ans,x1,x2,x3,x4,x5 + + ans,x1,x2,x3,x4,x5 + x1,x2,x3,x4,x5 + + + ans,x1,x2,x3,x4,x5 + x1,x2,x3,x4,x5 + + + ans,x1,x2,x3,x4,x5 + x1,x2,x3,x4,x5 + + + ans,x1,x2,x3,x4,x5 + x1,x2,x3,x4,x5 + + + ans,x1,x2,x3,x4,x5 + x1,x2,x3,x4,x5 + + + + + + + dummyIN + poly + MC_external + collset + + + collset + InputOutput + subdomainSensitivity + subdomainSensitivity + + + + + + csv + subdomainSensitivity + + + + + skew_ans, + skew_x1, + skew_x2, + skew_x3, + skew_x4, + skew_x5, + vc_ans, + vc_x1, + vc_x2, + vc_x3, + vc_x4, + vc_x5, + percentile_5_ans, + percentile_95_ans, + percentile_5_x1, + percentile_95_x1, + percentile_5_x2, + percentile_95_x2, + percentile_5_x3, + percentile_95_x3, + percentile_5_x4, + percentile_95_x4, + percentile_5_x5, + percentile_95_x5, + mean_ans, + mean_x1, + mean_x2, + mean_x3, + mean_x4, + mean_x5, + kurt_ans, + kurt_x1, + kurt_x2, + kurt_x3, + kurt_x4, + kurt_x5, + median_ans, + median_x1, + median_x2, + median_x3, + median_x4, + median_x5, + max_ans, + max_x1, + max_x2, + max_x3, + max_x4, + max_x5, + min_ans, + min_x1, + min_x2, + min_x3, + min_x4, + min_x5, + samp_ans, + samp_x1, + samp_x2, + samp_x3, + samp_x4, + samp_x5, + var_ans, + var_x1, + var_x2, + var_x3, + var_x4, + var_x5, + sigma_ans, + sigma_x1, + sigma_x2, + sigma_x3, + sigma_x4, + sigma_x5, + nsen_ans_x1, + nsen_ans_x2, + nsen_ans_x3, + nsen_ans_x4, + nsen_ans_x5, + nsen_x1_x1, + nsen_x1_x2, + nsen_x1_x3, + nsen_x1_x4, + nsen_x1_x5, + nsen_x2_x1, + nsen_x2_x2, + nsen_x2_x3, + nsen_x2_x4, + nsen_x2_x5, + nsen_x3_x1, + nsen_x3_x2, + nsen_x3_x3, + nsen_x3_x4, + nsen_x3_x5, + nsen_x4_x1, + nsen_x4_x2, + nsen_x4_x3, + nsen_x4_x4, + nsen_x4_x5, + nsen_x5_x1, + nsen_x5_x2, + nsen_x5_x3, + nsen_x5_x4, + nsen_x5_x5, + sen_ans_x1, + sen_ans_x2, + sen_ans_x3, + sen_ans_x4, + sen_ans_x5, + sen_x1_x1, + sen_x1_x2, + sen_x1_x3, + sen_x1_x4, + sen_x1_x5, + sen_x2_x1, + sen_x2_x2, + sen_x2_x3, + sen_x2_x4, + sen_x2_x5, + sen_x3_x1, + sen_x3_x2, + sen_x3_x3, + sen_x3_x4, + sen_x3_x5, + sen_x4_x1, + sen_x4_x2, + sen_x4_x3, + sen_x4_x4, + sen_x4_x5, + sen_x5_x1, + sen_x5_x2, + sen_x5_x3, + sen_x5_x4, + sen_x5_x5, + pear_ans_x1, + pear_ans_x2, + pear_ans_x3, + pear_ans_x4, + pear_ans_x5, + pear_x1_x1, + pear_x1_x2, + pear_x1_x3, + pear_x1_x4, + pear_x1_x5, + pear_x2_x1, + pear_x2_x2, + pear_x2_x3, + pear_x2_x4, + pear_x2_x5, + pear_x3_x1, + pear_x3_x2, + pear_x3_x3, + pear_x3_x4, + pear_x3_x5, + pear_x4_x1, + pear_x4_x2, + pear_x4_x3, + pear_x4_x4, + pear_x4_x5, + pear_x5_x1, + pear_x5_x2, + pear_x5_x3, + pear_x5_x4, + pear_x5_x5, + cov_ans_x1, + cov_ans_x2, + cov_ans_x3, + cov_ans_x4, + cov_ans_x5, + cov_x1_x1, + cov_x1_x2, + cov_x1_x3, + cov_x1_x4, + cov_x1_x5, + cov_x2_x1, + cov_x2_x2, + cov_x2_x3, + cov_x2_x4, + cov_x2_x5, + cov_x3_x1, + cov_x3_x2, + cov_x3_x3, + cov_x3_x4, + cov_x3_x5, + cov_x4_x1, + cov_x4_x2, + cov_x4_x3, + cov_x4_x4, + cov_x4_x5, + cov_x5_x1, + cov_x5_x2, + cov_x5_x3, + cov_x5_x4, + cov_x5_x5, + vsen_ans_x1, + vsen_ans_x2, + vsen_ans_x3, + vsen_ans_x4, + vsen_ans_x5, + vsen_x1_x1, + vsen_x1_x2, + vsen_x1_x3, + vsen_x1_x4, + vsen_x1_x5, + vsen_x2_x1, + vsen_x2_x2, + vsen_x2_x3, + vsen_x2_x4, + vsen_x2_x5, + vsen_x3_x1, + vsen_x3_x2, + vsen_x3_x3, + vsen_x3_x4, + vsen_x3_x5, + vsen_x4_x1, + vsen_x4_x2, + vsen_x4_x3, + vsen_x4_x4, + vsen_x4_x5, + vsen_x5_x1, + vsen_x5_x2, + vsen_x5_x3, + vsen_x5_x4, + vsen_x5_x5 + + + diff --git a/tests/framework/PostProcessors/SubdomainBasicStatistics/tests b/tests/framework/PostProcessors/SubdomainBasicStatistics/tests new file mode 100644 index 0000000000..8874543a0f --- /dev/null +++ b/tests/framework/PostProcessors/SubdomainBasicStatistics/tests @@ -0,0 +1,18 @@ +[Tests] + [./subdomainSensitivity] + type = 'RavenFramework' + input = 'sensitivity_subdomain.xml' + UnorderedCsv = 'subdomainStats/subdomainSensitivity.csv' + rel_err = 0.00001 + zero_threshold = 1e-9 + [../] + [./subdomainTimeDepStats] + type = 'RavenFramework' + input = 'time_dep_subdomain.xml' + UnorderedCsv = 'subdomainTimeDepStats/timedepSubdomainStats_dump.csv subdomainTimeDepStats/timedepSubdomainStats_dump_0.csv' + rel_err = 0.00001 + zero_threshold = 1e-9 + [../] + + +[] diff --git a/tests/framework/PostProcessors/SubdomainBasicStatistics/time_dep_subdomain.xml b/tests/framework/PostProcessors/SubdomainBasicStatistics/time_dep_subdomain.xml new file mode 100644 index 0000000000..a1161aaee6 --- /dev/null +++ b/tests/framework/PostProcessors/SubdomainBasicStatistics/time_dep_subdomain.xml @@ -0,0 +1,134 @@ + + + + subdomainTimeDepStats + FirstMRun,timeDepBasicStatPP + 1 + + + + framework/PostProcessors/SubdomainBasicStatistics.subdomainTimeDepStats + aalfonsi + 2023-05-17 + PostProcessors.SubdomainBasicStatistics + + This test checks the statistics and sensitivities (and other metrics) calculated + by subdomain basic statistics PP on + time-dependent data. This test shows how to compute subdomain statistics + using a 2D grid defined in the input space on time-dependent (HistorySet) data. + + + added test + + + + ../../BasicStatistics/basicStatisticsTimeDependent/samples.csv + + + + x,y,z,time,x0,x01,x02,y0,y02,y01,z0,z02,z01 + + + + + 10 100 + + + 10 50 + + + time + x0,y0,z0,x,y,z + + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + + + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + x0,y0,z0,x,y,z + + + + + + samples + + + + + + + + + inputPlaceholder + model + customSamplerFile + timedepSubdomainStats + + + timedepSubdomainStats + timeDepBasicStat + basicStatHistory + timedepSubdomainStats_dump + + + + + + csv + basicStatHistory + input,output + + + + + + x0,y0,z0 + OutputPlaceHolder + + + x0,y0,z0 + time,x,y,z + + + x0,y0 + + var_x0,var_y0,var_z0,var_x,var_y,var_z, + mean_x0, mean_y0, mean_z0, mean_x, mean_y, mean_z, + sigma_x0, sigma_y0, sigma_z0, sigma_x, sigma_y, sigma_z, + vc_x0, vc_y0, vc_z0, vc_x, vc_y, vc_z, + skew_x0, skew_y0, skew_z0, skew_x, skew_y, skew_z, + kurt_x0, kurt_y0, kurt_z0, kurt_x, kurt_y, kurt_z, + median_x0, median_y0, median_z0, median_x, median_y, median_z, + samp_x0, samp_y0, samp_z0, samp_x, samp_y, samp_z, + percentile_5_x0, percentile_5_y0, percentile_5_z0, percentile_5_x, percentile_5_y, percentile_5_z, + percentile_95_x0, percentile_95_y0, percentile_95_z0, percentile_95_x, percentile_95_y, percentile_95_z, + cov_x0_x0, cov_y0_x0, cov_z0_x0, cov_x_x0, cov_y_x0, cov_z_x0, + cov_x0_y0, cov_y0_y0, cov_z0_y0, cov_x_y0, cov_y_y0, cov_z_y0, + cov_x0_z0, cov_y0_z0, cov_z0_z0, cov_x_z0, cov_y_z0, cov_z_z0, + cov_x0_x, cov_y0_x, cov_z0_x, cov_x_x, cov_y_x, cov_z_x, + cov_x0_y, cov_y0_y, cov_z0_y, cov_x_y, cov_y_y, cov_z_y, + cov_x0_z, cov_y0_z, cov_z0_z, cov_x_z, cov_y_z, cov_z_z, + pearson_x0_x0, pearson_y0_x0, pearson_z0_x0, pearson_x_x0, pearson_y_x0, pearson_z_x0, + pearson_x0_y0, pearson_y0_y0, pearson_z0_y0, pearson_x_y0, pearson_y_y0, pearson_z_y0, + pearson_x0_z0, pearson_y0_z0, pearson_z0_z0, pearson_x_z0, pearson_y_z0, pearson_z_z0, + pearson_x0_x, pearson_y0_x, pearson_z0_x, pearson_x_x, pearson_y_x, pearson_z_x, + pearson_x0_y, pearson_y0_y, pearson_z0_y, pearson_x_y, pearson_y_y, pearson_z_y, + pearson_x0_z, pearson_y0_z, pearson_z0_z, pearson_x_z, pearson_y_z, pearson_z_z + + + time + + + + + From d3e4c2160b3a7e523bb9f67db387573e9b521f4c Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Mon, 22 May 2023 13:07:54 -0600 Subject: [PATCH 04/10] added documentations --- .../PostProcessors/BasicStatistics.tex | 14 ++- .../SubdomainBasicStatistics.tex | 112 ++++++++++++++++++ doc/user_manual/postprocessor.tex | 4 + doc/user_manual/raven_user_manual.tex | 13 ++ doc/user_manual/sampler.tex | 43 +++++++ 5 files changed, 183 insertions(+), 3 deletions(-) create mode 100644 doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex diff --git a/doc/user_manual/PostProcessors/BasicStatistics.tex b/doc/user_manual/PostProcessors/BasicStatistics.tex index 428ba493dd..1d61910442 100644 --- a/doc/user_manual/PostProcessors/BasicStatistics.tex +++ b/doc/user_manual/PostProcessors/BasicStatistics.tex @@ -12,7 +12,9 @@ \subsubsection{BasicStatistics} \textit{\textbf{Interfaced}} post-processor of type \textbf{HistorySetSync}). \end{itemize} % -\ppType{BasicStatistics post-processor}{BasicStatistics} + +\newcommand{\basicStatisticsBody} +{ \begin{itemize} \item \xmlNode{"metric"}, \xmlDesc{comma separated string or node list, required field}, specifications for the metric to be calculated. The name of each node is the requested metric. There are @@ -123,7 +125,7 @@ \subsubsection{BasicStatistics} computed and stored in the given \textbf{DataSet}, and the users do not need to specify these quantities in their RAVEN input files. % \item \xmlNode{pivotParameter}, \xmlDesc{string, optional field}, name of the parameter that needs - to be used for the computation of the Dynamic BasicStatistics (e.g. time). This node needs to + to be used for the computation of the Dynamic statistics (e.g. time). This node needs to be inputted just in case an \textbf{HistorySet} is used as Input. It represents the reference monotonic variable based on which the statistics is going to be computed (e.g. time-dependent statistical moments). @@ -146,7 +148,7 @@ \subsubsection{BasicStatistics} \item vector metrics, such as \xmlNode{covariance} and \xmlNode{sensitivity}, are requested, the index variables \xmlString{targets} and \xmlString{features} will be required. \item If \xmlNode{percentile} is requested, an additional index variable \xmlString{percent} should be added. - \item when dynamic BasicStatistics (e.g. time) is requested, the index variable \xmlString{time} will be required. + \item when dynamic statistics (e.g. time) is requested, the index variable \xmlString{time} will be required. \end{itemize} \default{False} % @@ -164,6 +166,12 @@ \subsubsection{BasicStatistics} \textbf{Example (Static Statistics):} This example demonstrates how to request the expected value of \xmlString{x01} and \xmlString{x02}, along with the sensitivity of both \xmlString{x01} and \xmlString{x02} to \xmlString{a} and \xmlString{b}. +} + +% in the following we have the description of the statistics +\ppType{BasicStatistics post-processor}{BasicStatistics} +\basicStatisticsBody + \begin{lstlisting}[style=XML,morekeywords={name,subType,debug}] ... diff --git a/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex b/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex new file mode 100644 index 0000000000..9bd682d0a2 --- /dev/null +++ b/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex @@ -0,0 +1,112 @@ +\subsubsection{SubdomainBasicStatistics} +\label{SubdomainBasicStatistics} +The \textbf{SubdomainBasicStatistics} post-processor is aimed to allow the \textbf{BasicStatistics} +post-processor to be used on subdomains of the sampling space. This post-processor fully leverages +the \textbf{BasicStatistics} post-processor and, consequentially, computes many of the most important statistical quantities on +a mesh/grid (subdomains). +Similarly to the \textbf{BasicStatistics}, the post-processor can accept as input both \textit{\textbf{PointSet}} +and \textit{\textbf{HistorySet}}. + +\ppType{SubdomainBasicStatistics post-processor}{SubdomainBasicStatistics} + +\begin{itemize} + \item \xmlNode{"subdomain"}, \xmlDesc{XML node, required parameter}, definition of the subdomain grid. This node + can specify the following XML node: + \begin{itemize} + \item \xmlNode{variable}, \xmlDesc{XML node, required parameter} can specify the \xmlAttr{name} attribute, which is \xmlDesc{required string attribute} and specifies the user-defined name of this variable. + \variableChildrenIntro + \begin{itemize} + \item \gridDescriptionNoCFD + \end{itemize} + \end{itemize} +\end{itemize} + +Based on the definition above, the post-processor will compute \textbf{BasicStatistics} on each ``cell'' defined +by the discretization grid. The results will be exported as function of the variables used +in the definition of the grid (see below for a detailed example). + +In addition to the subdomain grid defined above, the same nodes/settings available for \textbf{BasicStatistics} +are available: + +\basicStatisticsBody + +\textbf{Example (Static Subdomain Statistics):} +\begin{lstlisting}[style=XML,morekeywords={name,subType,debug}] + + ... + + ... + + + + 10 100 + + + 10 50 + + + x01,x02 + + x01,x02 + a,b + + + ... + + ... + + + a,b + + + mean_x01,mean_x02, + sen_x01_a, sen_x01_b, + sen_x02_a, sen_x02_b + + + +\end{lstlisting} + +In this case, the RAVEN variables ``mean\_x01, mean\_x02, sen\_x01\_a, sen\_x02\_a, sen\_x01\_b, sen\_x02\_b'' +will be computed for each cell defined above (2 cells in this case). +The results will be associated to the mid-point of each point. For the example above: +\begin{itemize} + \item $(a=32.5,b=30.0)$ + \item $(a=77.5,b=30.0)$ +\end{itemize} +The results will be then accessible for the RAVEN entities \textbf{DataObjects} and \textbf{OutStreams}. + +\textbf{Example (Dynamic Subdomain Statistics):} + +\begin{lstlisting}[style=XML,morekeywords={name,subType,debug}] + + ... + + ... + + x01,x02 + + x01,x02 + a,b + + time + + ... + + ... + + + a,b + + + mean_x01,mean_x02, + sen_x01_a, sen_x01_b, + sen_x02_a, sen_x02_b + + + time + + + +\end{lstlisting} + diff --git a/doc/user_manual/postprocessor.tex b/doc/user_manual/postprocessor.tex index 2f533a96fd..e884324241 100644 --- a/doc/user_manual/postprocessor.tex +++ b/doc/user_manual/postprocessor.tex @@ -12,6 +12,7 @@ \subsection{PostProcessor} \begin{itemize} \itemsep0em \item \textbf{BasicStatistics} + \item \textbf{SubdomainBasicStatistics} \item \textbf{ComparisonStatistics} \item \textbf{ImportanceRank} \item \textbf{SafestPoint} @@ -69,6 +70,9 @@ \subsection{PostProcessor} %%%%% PP BasicStatistics %%%%%%% \input{PostProcessors/BasicStatistics.tex} +%%%%% PP SubdomainBasicStatistics %%%%%%% +\input{PostProcessors/SubdomainBasicStatistics.tex} + %%%%% PP ComparisonStatistics %%%%%%% \input{PostProcessors/ComparisonStatistics.tex} diff --git a/doc/user_manual/raven_user_manual.tex b/doc/user_manual/raven_user_manual.tex index f63af1ec67..a16c63826c 100644 --- a/doc/user_manual/raven_user_manual.tex +++ b/doc/user_manual/raven_user_manual.tex @@ -90,6 +90,19 @@ \usepackage{listings} +% the following allows for 8 levels list nesting +\usepackage{enumitem} +\setlistdepth{8} +\setlist[itemize,1]{label=$\bullet$} +\setlist[itemize,2]{label=$\bullet$} +\setlist[itemize,3]{label=$\bullet$} +\setlist[itemize,4]{label=$\bullet$} +\setlist[itemize,5]{label=$\bullet$} +\setlist[itemize,6]{label=$\bullet$} +\setlist[itemize,7]{label=$\bullet$} +\setlist[itemize,8]{label=$\bullet$} +\renewlist{itemize}{itemize}{8} + % Python style for highlighting \newcommand\pythonstyle{\lstset{ language=Python, diff --git a/doc/user_manual/sampler.tex b/doc/user_manual/sampler.tex index 0f696adebb..ab23073c20 100644 --- a/doc/user_manual/sampler.tex +++ b/doc/user_manual/sampler.tex @@ -118,6 +118,34 @@ \section{Samplers} code will raise an error. \end{itemize} } +\newcommand{\constructionGridDescriptionNoCFD} +{ +Based on the \xmlAttr{construction} type, the content of the \xmlNode{grid} +XML node and the requirements for other attributes change: +\begin{itemize} + \item \xmlAttr{construction}\textbf{\texttt{=}}\xmlString{equal}. + The grid is going to be constructed equally-spaced + (\xmlAttr{type}\textbf{\texttt{=}}\xmlString{value}). + This construction type requires the definition of additional attributes: + \begin{itemize} + \item \xmlAttr{steps}, \xmlDesc{required integer attribute}, number + of equally spaced/probable discretization steps. + \end{itemize} + This construction type requires that the content of the \xmlNode{grid} + node represents the lower and upper bounds (either + in probability or value). Two values need to be specified; the lowest one + will be considered as the $lowerBound$, the largest, the $upperBound$. + The $stepSize$ is determined as follows: + \\ $stepSize=(upperBound - lowerBound)/steps$ + \item \xmlAttr{construction}\textbf{\texttt{=}}\xmlString{custom}. + The grid will be directly specified by the user. + No additional attributes are needed. + This construction type requires that the \xmlNode{grid} node contains + the actual mesh bins. + For example, if the grid \xmlAttr{type} is \xmlString{value}, in the body + of \xmlNode{grid}, the user will specify the values representing the nodes of the grid. +\end{itemize} +} \newcommand{\variableDescription} { \xmlNode{variable}, \xmlDesc{XML node, required parameter} can specify the following attribute: @@ -234,6 +262,21 @@ \section{Samplers} information is requested.} } +\newcommand{\gridDescriptionNoCFD} +{ + \xmlNode{grid}, \xmlDesc{space separated floats, required + field}, the content of this XML node allows for a type of the grid only: + \begin{itemize} + \itemsep0em + \item \xmlAttr{type}, \xmlDesc{required string attribute}, user-defined + discretization metric type: 1) \xmlString{value}, the grid will be provided + using variable values. + \item \xmlAttr{construction}, \xmlDesc{required string attribute}, how + the grid needs to be constructed, independent of its type. + \end{itemize} + \constructionGridDescriptionNoCFD +} + \newcommand{\convergenceDescription} { \xmlNode{Convergence}, \xmlDesc{float, required field}, Convergence From 169c687437ddf36b50ec6e426c43209526a648fa Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Mon, 22 May 2023 13:14:22 -0600 Subject: [PATCH 05/10] removed unused imports --- .../Models/PostProcessors/SubdomainBasicStatistics.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py index d5ffc7fd85..c6d1e7e384 100644 --- a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py +++ b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py @@ -18,12 +18,6 @@ """ #External Modules--------------------------------------------------------------- import numpy as np -import os -import copy -from collections import OrderedDict, defaultdict -import six -import xarray as xr -import scipy.stats as stats #External Modules End----------------------------------------------------------- #Internal Modules--------------------------------------------------------------- @@ -31,7 +25,6 @@ from .BasicStatistics import BasicStatistics from ...utils import utils from ...utils import InputData, InputTypes -from ...utils import mathUtils #Internal Modules End----------------------------------------------------------- class SubdomainBasicStatistics(PostProcessorInterface): From a262164aacceec207a09c39d1df2decf15776ca2 Mon Sep 17 00:00:00 2001 From: Andrea Alfonsi Date: Mon, 22 May 2023 13:32:31 -0600 Subject: [PATCH 06/10] Apply suggestions from code review --- ravenframework/Models/PostProcessors/BasicStatistics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ravenframework/Models/PostProcessors/BasicStatistics.py b/ravenframework/Models/PostProcessors/BasicStatistics.py index 25f29dab6e..5c5b07aa35 100644 --- a/ravenframework/Models/PostProcessors/BasicStatistics.py +++ b/ravenframework/Models/PostProcessors/BasicStatistics.py @@ -1563,7 +1563,6 @@ def run(self, inputIn): """ inputData = self.inputToInternal(inputIn) outputSet = self._runLocal(inputData) - print(outputSet) return outputSet def collectOutput(self, finishedJob, output): From 12024fd20815539236a3a5a086949a93de7d8121 Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Mon, 22 May 2023 14:32:48 -0600 Subject: [PATCH 07/10] readded self.outputMultipleRealizations --- ravenframework/Models/PostProcessors/PostProcessorInterface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ravenframework/Models/PostProcessors/PostProcessorInterface.py b/ravenframework/Models/PostProcessors/PostProcessorInterface.py index f540d2177a..ed21f8cd3f 100644 --- a/ravenframework/Models/PostProcessors/PostProcessorInterface.py +++ b/ravenframework/Models/PostProcessors/PostProcessorInterface.py @@ -66,7 +66,7 @@ def __init__(self): ## However, the DataObject.load can not be directly used to collect single realization ## One possible solution is all postpocessors return a list of realizations, and we only ## use addRealization method to add the collections into the DataObjects - + self.outputMultipleRealizations = False def _handleInput(self, paramInput): """ From 138c038d297a3bb7dd44726dd15f7c5f5bc91030 Mon Sep 17 00:00:00 2001 From: Andrea Alfonsi Date: Wed, 24 May 2023 08:32:52 -0600 Subject: [PATCH 08/10] Update doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex --- doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex b/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex index 9bd682d0a2..0ba0a775eb 100644 --- a/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex +++ b/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex @@ -69,7 +69,7 @@ \subsubsection{SubdomainBasicStatistics} In this case, the RAVEN variables ``mean\_x01, mean\_x02, sen\_x01\_a, sen\_x02\_a, sen\_x01\_b, sen\_x02\_b'' will be computed for each cell defined above (2 cells in this case). -The results will be associated to the mid-point of each point. For the example above: +The results will be associated to the mid-point of each cell. For the example above: \begin{itemize} \item $(a=32.5,b=30.0)$ \item $(a=77.5,b=30.0)$ From 36bf939b9ac1e2834b52eb89f3cbd67542c34e2a Mon Sep 17 00:00:00 2001 From: Andrea Alfonsi Date: Wed, 24 May 2023 08:41:20 -0600 Subject: [PATCH 09/10] Update ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py --- .../Models/PostProcessors/SubdomainBasicStatistics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py index c6d1e7e384..646294e2e1 100644 --- a/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py +++ b/ravenframework/Models/PostProcessors/SubdomainBasicStatistics.py @@ -193,7 +193,7 @@ def run(self, inputIn): def collectOutput(self, finishedJob, output): """ Function to place all of the computed data into the output object - @ In, finishedJob, JobHandler External or Internal instance, A JobHandler object that is in charge of running this post-processor + @ In, finishedJob, JobHandler.Runner instance, the instance containing the completed job @ In, output, dataObjects, The object where we want to place our computed results @ Out, None """ From b0a2f97301a155c70f95bad28a75c97394ae5522 Mon Sep 17 00:00:00 2001 From: aalfonsi Date: Wed, 24 May 2023 09:58:06 -0600 Subject: [PATCH 10/10] moved latex outside sampler --- .../SubdomainBasicStatistics.tex | 2 +- doc/user_manual/postprocessor.tex | 45 +++++++++++++++++ doc/user_manual/sampler.tex | 50 +------------------ 3 files changed, 48 insertions(+), 49 deletions(-) diff --git a/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex b/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex index 9bd682d0a2..c4f279f203 100644 --- a/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex +++ b/doc/user_manual/PostProcessors/SubdomainBasicStatistics.tex @@ -16,7 +16,7 @@ \subsubsection{SubdomainBasicStatistics} \item \xmlNode{variable}, \xmlDesc{XML node, required parameter} can specify the \xmlAttr{name} attribute, which is \xmlDesc{required string attribute} and specifies the user-defined name of this variable. \variableChildrenIntro \begin{itemize} - \item \gridDescriptionNoCFD + \item \gridDescriptionNoCDF \end{itemize} \end{itemize} \end{itemize} diff --git a/doc/user_manual/postprocessor.tex b/doc/user_manual/postprocessor.tex index e884324241..0007eec2fd 100644 --- a/doc/user_manual/postprocessor.tex +++ b/doc/user_manual/postprocessor.tex @@ -1,5 +1,50 @@ \subsection{PostProcessor} \label{sec:models_postProcessor} + +\newcommand{\constructionNoCDF} +{ +Based on the \xmlAttr{construction} type, the content of the \xmlNode{grid} +XML node and the requirements for other attributes change: +\begin{itemize} + \item \xmlAttr{construction}\textbf{\texttt{=}}\xmlString{equal}. + The grid is going to be constructed equally-spaced + (\xmlAttr{type}\textbf{\texttt{=}}\xmlString{value}). + This construction type requires the definition of additional attributes: + \begin{itemize} + \item \xmlAttr{steps}, \xmlDesc{required integer attribute}, number + of equally spaced/probable discretization steps. + \end{itemize} + This construction type requires that the content of the \xmlNode{grid} + node represents the lower and upper bounds (either + in probability or value). Two values need to be specified; the lowest one + will be considered as the $lowerBound$, the largest, the $upperBound$. + The $stepSize$ is determined as follows: + \\ $stepSize=(upperBound - lowerBound)/steps$ + \item \xmlAttr{construction}\textbf{\texttt{=}}\xmlString{custom}. + The grid will be directly specified by the user. + No additional attributes are needed. + This construction type requires that the \xmlNode{grid} node contains + the actual mesh bins. + For example, if the grid \xmlAttr{type} is \xmlString{value}, in the body + of \xmlNode{grid}, the user will specify the values representing the nodes of the grid. +\end{itemize} +} + +\newcommand{\gridDescriptionNoCDF} +{ + \xmlNode{grid}, \xmlDesc{space separated floats, required + field}, the content of this XML node allows for a type of the grid only: + \begin{itemize} + \itemsep0em + \item \xmlAttr{type}, \xmlDesc{required string attribute}, user-defined + discretization metric type: 1) \xmlString{value}, the grid will be provided + using variable values. + \item \xmlAttr{construction}, \xmlDesc{required string attribute}, how + the grid needs to be constructed, independent of its type. + \end{itemize} + \constructionNoCDF +} + A Post-Processor (PP) can be considered as an action performed on a set of data or other type of objects. % diff --git a/doc/user_manual/sampler.tex b/doc/user_manual/sampler.tex index ab23073c20..b406853895 100644 --- a/doc/user_manual/sampler.tex +++ b/doc/user_manual/sampler.tex @@ -118,34 +118,7 @@ \section{Samplers} code will raise an error. \end{itemize} } -\newcommand{\constructionGridDescriptionNoCFD} -{ -Based on the \xmlAttr{construction} type, the content of the \xmlNode{grid} -XML node and the requirements for other attributes change: -\begin{itemize} - \item \xmlAttr{construction}\textbf{\texttt{=}}\xmlString{equal}. - The grid is going to be constructed equally-spaced - (\xmlAttr{type}\textbf{\texttt{=}}\xmlString{value}). - This construction type requires the definition of additional attributes: - \begin{itemize} - \item \xmlAttr{steps}, \xmlDesc{required integer attribute}, number - of equally spaced/probable discretization steps. - \end{itemize} - This construction type requires that the content of the \xmlNode{grid} - node represents the lower and upper bounds (either - in probability or value). Two values need to be specified; the lowest one - will be considered as the $lowerBound$, the largest, the $upperBound$. - The $stepSize$ is determined as follows: - \\ $stepSize=(upperBound - lowerBound)/steps$ - \item \xmlAttr{construction}\textbf{\texttt{=}}\xmlString{custom}. - The grid will be directly specified by the user. - No additional attributes are needed. - This construction type requires that the \xmlNode{grid} node contains - the actual mesh bins. - For example, if the grid \xmlAttr{type} is \xmlString{value}, in the body - of \xmlNode{grid}, the user will specify the values representing the nodes of the grid. -\end{itemize} -} + \newcommand{\variableDescription} { \xmlNode{variable}, \xmlDesc{XML node, required parameter} can specify the following attribute: @@ -203,6 +176,7 @@ \section{Samplers} the attribute \xmlAttr{dim} is required. \nb{Alternatively, this node must be omitted if the \xmlNode{function} node is supplied.} } + \newcommand{\functionDescription} { \xmlNode{function}, \xmlDesc{string, required field}, name of the function that @@ -213,9 +187,6 @@ \section{Samplers} if the \xmlNode{distribution} node is supplied.} } - - - \newcommand{\gridDescriptionOnlyCustom} { \xmlNode{grid}, \xmlDesc{space separated floats, required @@ -238,8 +209,6 @@ \section{Samplers} information is requested.} } - - \newcommand{\gridDescription} { \xmlNode{grid}, \xmlDesc{space separated floats, required @@ -262,21 +231,6 @@ \section{Samplers} information is requested.} } -\newcommand{\gridDescriptionNoCFD} -{ - \xmlNode{grid}, \xmlDesc{space separated floats, required - field}, the content of this XML node allows for a type of the grid only: - \begin{itemize} - \itemsep0em - \item \xmlAttr{type}, \xmlDesc{required string attribute}, user-defined - discretization metric type: 1) \xmlString{value}, the grid will be provided - using variable values. - \item \xmlAttr{construction}, \xmlDesc{required string attribute}, how - the grid needs to be constructed, independent of its type. - \end{itemize} - \constructionGridDescriptionNoCFD -} - \newcommand{\convergenceDescription} { \xmlNode{Convergence}, \xmlDesc{float, required field}, Convergence