Mandd/margin solver enh (#49)

* edit src file and new test

* edit in user manual

* delete whitespaces

* edit files for review

* added BElist capability

* edits

* addressing comemnts

* fixed test

* delete whitespaces

* fixed tests

* Update MCSSolver.py

* delete whitespaces

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

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,16 +145,22 @@ 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()
       elif input.type == 'PointSet':
         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,29 +248,62 @@ 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
 
     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

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

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

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
     </description>
   </TestInfo>
 
@@ -46,6 +46,7 @@
       <variables>statusA,statusB,statusC,statusD,statusE,TOP</variables>
       <solver type='margin'>
         <metric>2</metric>
+        <setType>cut</setType>
       </solver>
       <topEventID>TOP</topEventID>
       <map var='statusA'>A</map>

tests/test_MCSSolver_margin_pathsets.xml

@@ -0,0 +1,85 @@
+<Simulation verbosity="debug">
+  <TestInfo>
+    <name>SR2ML/tests.MCSSolver_margin_pathsets</name>
+    <author>mandd</author>
+    <created>2021-04-27</created>
+    <classesTested>SR2ML.MCSSolver</classesTested>
+    <description>
+      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
+    </description>
+  </TestInfo>
+
+  <RunInfo>
+    <WorkingDir>MCSSolverMargin</WorkingDir>
+    <Sequence>simRun</Sequence>
+    <batchSize>1</batchSize>
+  </RunInfo>
+
+  <Files>
+    <Input name="MCSlistFile" type="MCSlist">MCSlist.csv</Input>
+  </Files>
+
+  <Models>
+    <ExternalModel name="MCSmodel" subType="SR2ML.MCSSolver">
+      <inputs>statusA,statusB,statusC,statusD,statusE</inputs>
+      <outputs>TOP,sens_statusA,sens_statusB,sens_statusC,sens_statusD,sens_statusE</outputs>
+      <solver type='margin'>
+        <metric>2</metric>
+        <setType>path</setType>
+      </solver>
+      <topEventID>TOP</topEventID>
+      <map var='statusA'>A</map>
+      <map var='statusB'>B</map>
+      <map var='statusC'>C</map>
+      <map var='statusD'>D</map>
+      <map var='statusE'>E</map>
+    </ExternalModel>
+  </Models>
+
+  <Samplers>
+    <MonteCarlo name="MC_external">
+      <samplerInit>
+        <limit>1</limit>
+      </samplerInit>
+      <constant name="statusA">0.35</constant>
+      <constant name="statusB">0.2</constant>
+      <constant name="statusC">0.1</constant>
+      <constant name="statusD">0.3</constant>
+      <constant name="statusE">0.4</constant>
+    </MonteCarlo>
+  </Samplers>
+
+  <Steps>
+    <MultiRun name="simRun">
+      <Input   class="Files"        type="MCSlist"         >MCSlistFile</Input>
+      <Model   class="Models"       type="ExternalModel"   >MCSmodel</Model>
+      <Sampler class="Samplers"     type="MonteCarlo"      >MC_external</Sampler>
+      <Output  class="DataObjects"  type="PointSet"        >sim_PS</Output>
+      <Output  class="OutStreams"   type="Print"           >Print_sim_PS_path</Output>
+    </MultiRun>
+  </Steps>
+
+  <OutStreams>
+    <Print name="Print_sim_PS_path">
+      <type>csv</type>
+      <source>sim_PS</source>
+      <what>input,output</what>
+    </Print>
+  </OutStreams>
+
+  <DataObjects>
+    <PointSet name="inputPlaceHolder">
+      <Input>statusA,statusB,statusC,statusD,statusE</Input>
+      <Output>OutputPlaceHolder</Output>
+    </PointSet>
+    <PointSet name="sim_PS">
+      <Input>statusA,statusB,statusC,statusD,statusE</Input>
+      <Output>TOP,sens_statusA,sens_statusB,sens_statusC,sens_statusD,sens_statusE</Output>
+    </PointSet>
+  </DataObjects>
+
+</Simulation>

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'
+ [../]
 []