diff --git a/doc/user_manual/include/MCSsolver.tex b/doc/user_manual/include/MCSsolver.tex index e4de1a5..ec9f812 100644 --- a/doc/user_manual/include/MCSsolver.tex +++ b/doc/user_manual/include/MCSsolver.tex @@ -41,7 +41,8 @@ \section{MCSSolver} to be considered when evaluating the probability of their union \item \xmlNode{metric},\xmlDesc{string, required parameter}, metric used for margin calculations of the top event of $M(TE)$; it specifies the desired distance metric in terms of p values for the $L_p$ norm (allowed values - are: 0, 1, 2, inf). + are: 0, 1, 2, inf). + \item \xmlNode{setType},\xmlDesc{string, required parameter}, type of provided sets : minimal path sets (path) or minimal cut sets (cut) \end{itemize} \item \xmlNode{topEventID},\xmlDesc{string, required parameter}, the name of the alias variable for the Top Event \item \xmlNode{map},\xmlDesc{string, required parameter}, the name ID of the ET branching variable @@ -152,7 +153,7 @@ \subsection{Time dependent calculation} \end{itemize} \subsection{MCSSolver model reference tests} -The following is the provided analytic test: +The following is the provided analytic tests: \begin{itemize} \item test\_MCSSolver.xml \item test\_MCSSolver\_TD.xml diff --git a/src/MCSSolver.py b/src/MCSSolver.py index 65a9eb6..4640869 100644 --- a/src/MCSSolver.py +++ b/src/MCSSolver.py @@ -22,6 +22,7 @@ import itertools import math import xarray as xr +import pandas as pd import copy #External Modules End----------------------------------------------------------- @@ -45,6 +46,8 @@ def __init__(self): self.timeDepData = None # This variable contains the basic event temporal profiles as xr.Dataset self.topEventTerms = {} # Dictionary containing, for each order, a list of terms containing the union of MCSs self.mcsList = None # List containing all the MCSs; each MCS is a list of basic events + self.solver['setType'] = None # Type of sets provided path sets (path) or cut sets (cut) + self.fullListBE = None def initialize(self, container, runInfoDict, inputFiles): """ @@ -75,7 +78,12 @@ def _readMoreXML(self, container, xmlNode): container.beId = None # ID of the variable containing the IDs of the BEs container.tdFromPS = False # boolean variable which flags when TD calculation is generated from PS + self.solver['setType'] = None # type of set provided by the user: path sets for cut sets + self.fullListBE = None # list of BEs + metricOrder = {'0':0, '1':1, '2':2, 'inf':np.inf} + setTypes = ['path','cut'] + externalModelNodes = ['inputs','outputs'] for child in xmlNode: if child.tag == 'topEventID': @@ -103,14 +111,27 @@ def _readMoreXML(self, container, xmlNode): if metricValue not in metricOrder.keys(): raise IOError("MCSSolver: value in xml node metric is not allowed (0,1,2,inf)") self.solver['metric'] = metricOrder[metricValue] + elif childChild.tag == 'setType': + setType = childChild.text.strip() + if setType not in setTypes: + raise IOError("MCSSolver: set type in xml node setType is not allowed (cut,or path)") + self.solver['setType'] = setType elif child.tag == 'variables': variables = [str(var.strip()) for var in child.text.split(",")] elif child.tag == 'map': container.mapping[child.get('var')] = child.text.strip() container.invMapping[child.text.strip()] = child.get('var') - else: + elif child.tag not in externalModelNodes: raise IOError("MCSSolver: xml node " + str(child.tag) + " is not allowed") + if self.solver['setType'] == None: + ("MCSSolver: set type in xml node setType has not been specified (cut,or path)") + + if not bool(container.mapping): + variables.remove(container.topEventID) + for variable in (variables): + container.mapping[variable] = variable + container.invMapping[variable] = variable def createNewInput(self, container, inputs, samplerType, **kwargs): """ @@ -124,7 +145,6 @@ def createNewInput(self, container, inputs, samplerType, **kwargs): """ if len(inputs) > 2: raise IOError("MCSSolver: More than one file has been passed to the MCS solver") - for input in inputs: if input.type == 'HistorySet': self.timeDepData = input.asDataset() @@ -132,8 +152,15 @@ def createNewInput(self, container, inputs, samplerType, **kwargs): self.timeDepData = self.generateHistorySetFromSchedule(container,input.asDataset()) container.tdFromPS = True else: - mcsIDs, probability, self.mcsList, self.beList = mcsReader(input.getFilename(), type=container.fileFrom) - + if input.__getstate__()['type'] == 'MCSlist': + mcsIDs, probability, self.mcsList, self.beList = mcsReader(input.getFilename(), type=container.fileFrom) + elif input.__getstate__()['type'] == 'BElist': + self.fullListBE = beReader(input.getFilename()) + # Check list of BEs + if self.fullListBE is not None: + missingBEs = set(self.beList) - set(self.fullListBE) + if len(missingBEs)>0: + raise IOError("MCSSolver: there are BEs in the MCSs not defined in the BE file: " +str(missingBEs)) # mcsList is supposed to be a list of lists # E.g., if the MCS are ABC CD and AE --> MCS1=['A','B','C'], MCS2=['D','C'], MCS3=['A','E'] # then mcsList = [MCS1,MCS2,MCS3] = @@ -152,7 +179,7 @@ def createNewInput(self, container, inputs, samplerType, **kwargs): # (['D', 'C'], ['A', 'E']) ] basicEventCombined = list(set(itertools.chain.from_iterable(term)) for term in terms) - self.topEventTerms[order]=basicEventCombined + self.topEventTerms[order] = basicEventCombined return kwargs @@ -164,6 +191,7 @@ def run(self, container, inputs): @ In, inputs, dict, dictionary of inputs from RAVEN @ Out, None """ + if self.timeDepData is None: self.runStatic(container, inputs) else: @@ -185,6 +213,10 @@ def runStatic(self, container, inputs): topEventValue = self.mcsSolverProbability(inputForSolver) else: topEventValue = self.mcsSolverMargin(inputForSolver) + sensitivities = self.marginSensitivities(topEventValue, inputForSolver) + for key in sensitivities: + keyID = "sens_" + str(container.invMapping[key]) + container.__dict__[keyID] = sensitivities[key] container.__dict__[container.topEventID] = np.asarray(float(topEventValue)) @@ -198,6 +230,10 @@ def runDynamic(self, container, inputs): @ Out, None """ topEventValue = np.zeros([self.timeDepData[container.timeID].shape[0]]) + if self.solver['type'] == 'margin': + sensitivities = {} + for key in inputs: + sensitivities[key] = np.zeros([self.timeDepData[container.timeID].shape[0]]) for index,t in enumerate(self.timeDepData[container.timeID]): inputForSolver = {} @@ -212,6 +248,10 @@ def runDynamic(self, container, inputs): else: topEventValue[index] = self.mcsSolverMargin(inputForSolver) + sensValues = self.marginSensitivities(topEventValue, inputForSolver) + for key in sensitivities: + sensitivities[key][index] = sensValues[key] + if container.tdFromPS: for key in container.invMapping.keys(): container.__dict__[key] = self.timeDepData[key][0].values @@ -219,22 +259,51 @@ def runDynamic(self, container, inputs): container.__dict__[container.timeID] = self.timeDepData[container.timeID].values container.__dict__[container.topEventID] = topEventValue + if self.solver['type'] == 'margin': + for key in sensitivities: + keyID = "sens_" + str(container.invMapping[key]) + container.__dict__[keyID] = sensitivities[key] + + + def marginSensitivities(self, MsysBase, inputDict, epsilon = 0.01): + """ + This method calculates the sensitivity (derivative based) of the top event margin vs. basic event margin + @ In, inputDict, dictionary containing the probability value of all basic events + @ In, MsysBase, float, base top event margin + @ Out, sensDict, dict, dictionary containing the sensitivity values of all basic events + """ + sensDict={} + for key in inputDict: + tempDict = copy.deepcopy(inputDict) + # The sensitivitiy of each basic event is calculated in a derivative form as d M_sys/d M_be + # Thus it is needed to calculate: + # * M_sys with nominal value of M_be + # * M_sys with modified value of M_be as M_be*(1.-epsilon) + tempDict[key] = tempDict[key] * (1.-epsilon) + deltaMsys = self.mcsSolverMargin(tempDict) + sensDict[key] = (MsysBase - deltaMsys) / (inputDict[key] - tempDict[key]) + + return sensDict + def mcsSolverProbability(self, inputDict): """ This method determines the probability of the TopEvent of the FT provided the probability of its Basic Events - @ In, inputs, inputDict, dictionary containing the probability value of all basic events + @ In, inputDict, dictionary containing the probability value of all basic events @ Out, teProbability, float, probability value of the top event """ teProbability = 0.0 multiplier = 1.0 + if self.fullListBE is not None: + inputDict = {**self.fullListBE, **inputDict} # perform probability calculation for each order level for order in range(1,self.solver['solverOrder']+1): orderProbability=0 for term in self.topEventTerms[order]: # map the sampled values of the basic event probabilities to the MCS basic events ID termValues = list(map(inputDict.get,term)) + #print(term, termValues) orderProbability = orderProbability + np.prod(termValues) teProbability = teProbability + multiplier * orderProbability multiplier = -1.0 * multiplier @@ -249,12 +318,16 @@ def mcsSolverMargin(self, inputDict): @ Out, teMargin, float, margin value of the top event """ mcsMargins = np.zeros(len(self.mcsList)) - for index,mcs in enumerate(self.mcsList): termValues = list(map(inputDict.get,mcs)) - mcsMargins[index] = np.linalg.norm(termValues, ord=self.solver['metric']) - - teMargin = np.amin(mcsMargins) + if self.solver['setType']=='cut': + mcsMargins[index] = np.linalg.norm(termValues, ord=self.solver['metric']) + else: + mcsMargins[index] = np.amin(termValues) + if self.solver['setType']=='cut': + teMargin = np.amin(mcsMargins) + else: + teMargin = np.linalg.norm(mcsMargins, ord=self.solver['metric']) return teMargin @@ -286,3 +359,24 @@ def generateHistorySetFromSchedule(self, container, inputDataset): basicEventHistorySet[key][0][indexes] = 1.0 return basicEventHistorySet + + +def beReader(fileID): + """ + This method is designed to read the file containing the set of basic events and their associated probability value + @ In, fileID, file, file which is structured in two columns; the first containing the ID of the basic events, + the second one containing the corresponding probability values + @ Out, BEdict, dict, dictionary containing the ID of the basic events and the corresponding probability values + """ + data = pd.read_csv(fileID) + fromList = ['FALSE','TRUE'] + toList = [0.0,1.0] + data = data.replace(to_replace=fromList, value=toList) + pd.set_option('display.max_rows', None) + finalColumns = ['Name','Value'] + data = data[data.columns.intersection(finalColumns)] + data = data[finalColumns] + data['Value'] = data['Value'].astype(float) + data['Name'] = data['Name'].str.replace(' ','') + BEdict = dict(data.values) + return BEdict diff --git a/src/PostProcessors/MCSimporter.py b/src/PostProcessors/MCSimporter.py index 531f94d..2921680 100644 --- a/src/PostProcessors/MCSimporter.py +++ b/src/PostProcessors/MCSimporter.py @@ -162,7 +162,7 @@ def mcsReader(mcsListFile, type=None): elementsList.pop(0) probability=np.append(probability,float(elementsList[0])) elementsList.pop(0) - mcsList.append(list(element.rstrip('\n') for element in elementsList)) + mcsList.append(list(element.rstrip('\n').strip() for element in elementsList)) beList.update(elementsList) elif type.lower() == 'saphire': lines = lines[1:] # skip additional description line diff --git a/tests/gold/MCSSolverMargin/Print_sim_PS_path.csv b/tests/gold/MCSSolverMargin/Print_sim_PS_path.csv new file mode 100644 index 0000000..042e4b8 --- /dev/null +++ b/tests/gold/MCSSolverMargin/Print_sim_PS_path.csv @@ -0,0 +1,2 @@ +statusA,statusB,statusC,statusD,statusE,TOP,sens_statusA,sens_statusB,sens_statusC,sens_statusD,sens_statusE +0.35,0.2,0.1,0.3,0.4,0.471699056603,0.740321883138,0.845273020793,0.422067985733,0.0,0.0 diff --git a/tests/test_MCSSolver_margin.xml b/tests/test_MCSSolver_margin.xml index dad4d75..0593aec 100644 --- a/tests/test_MCSSolver_margin.xml +++ b/tests/test_MCSSolver_margin.xml @@ -27,7 +27,7 @@ This model is designed to read from file a list of Minimal Cut Sets (MCSs) and to import such Boolean logic structure as a RAVEN model. Provided the sampled values of Basic Events (BEs) margin values, the MCSSolver determines the - margin of Top Event (TE) using the euclidean metrics + margin of Top Event (TE) using the Euclidean metrics @@ -46,6 +46,7 @@ statusA,statusB,statusC,statusD,statusE,TOP 2 + cut TOP A diff --git a/tests/test_MCSSolver_margin_pathsets.xml b/tests/test_MCSSolver_margin_pathsets.xml new file mode 100644 index 0000000..34c76f2 --- /dev/null +++ b/tests/test_MCSSolver_margin_pathsets.xml @@ -0,0 +1,85 @@ + + + SR2ML/tests.MCSSolver_margin_pathsets + mandd + 2021-04-27 + SR2ML.MCSSolver + + This model is designed to read from file a list of Minimal Path Sets (MPSs) and + to import such Boolean logic structure as a RAVEN model. Provided the sampled + margin values of all Basic Events (BEs), the MCSSolver determines the + margin of Top Event (TE) using the Euclidean metrics and it calculates the + derivative based importance measures for each basic event + + + + + MCSSolverMargin + simRun + 1 + + + + MCSlist.csv + + + + + statusA,statusB,statusC,statusD,statusE + TOP,sens_statusA,sens_statusB,sens_statusC,sens_statusD,sens_statusE + + 2 + path + + TOP + A + B + C + D + E + + + + + + + 1 + + 0.35 + 0.2 + 0.1 + 0.3 + 0.4 + + + + + + MCSlistFile + MCSmodel + MC_external + sim_PS + Print_sim_PS_path + + + + + + csv + sim_PS + input,output + + + + + + statusA,statusB,statusC,statusD,statusE + OutputPlaceHolder + + + statusA,statusB,statusC,statusD,statusE + TOP,sens_statusA,sens_statusB,sens_statusC,sens_statusD,sens_statusE + + + + diff --git a/tests/tests b/tests/tests index 4c15ee4..b21cc26 100644 --- a/tests/tests +++ b/tests/tests @@ -113,4 +113,10 @@ input = 'test_MCSSolver_Saphire.xml' UnorderedCsv = 'MCSSolverSaphire/Print_sim_PS.csv' [../] + + [./TestMCSSolver_margin_paths] + type = 'RavenFramework' + input = 'test_MCSSolver_margin_pathsets.xml' + UnorderedCsv = 'MCSSolverMargin/Print_sim_PS_path.csv' + [../] []