From 0510b9b0eba7235443dec10435a076b0bfed5656 Mon Sep 17 00:00:00 2001 From: HugoNip Date: Fri, 27 Mar 2020 09:48:48 -0400 Subject: [PATCH 01/10] add new test code --- include/utils/FaultTreeUtils.h | 106 +++++ include/utils/QuantificationUtils.h | 189 ++++++++ include/utils/Utils.h | 19 + src/utils/FaultTreeUtils.C | 519 +++++++++++++++++++++ src/utils/QuantificationUtils.C | 679 ++++++++++++++++++++++++++++ unit/hazard.txt | 6 + unit/logic1.txt | 5 + unit/logic1_bas_events_LNORM.txt | 4 + unit/logic2.txt | 7 + unit/logic2_bas_events_PE.txt | 5 + unit/src/TestFaultTreeUtils.C | 176 +++++++ unit/src/TestQuantificationUtils.C | 456 +++++++++++++++++++ 12 files changed, 2171 insertions(+) create mode 100644 include/utils/FaultTreeUtils.h create mode 100644 include/utils/QuantificationUtils.h create mode 100644 include/utils/Utils.h create mode 100644 src/utils/FaultTreeUtils.C create mode 100644 src/utils/QuantificationUtils.C create mode 100644 unit/hazard.txt create mode 100644 unit/logic1.txt create mode 100644 unit/logic1_bas_events_LNORM.txt create mode 100644 unit/logic2.txt create mode 100644 unit/logic2_bas_events_PE.txt create mode 100644 unit/src/TestFaultTreeUtils.C create mode 100644 unit/src/TestQuantificationUtils.C diff --git a/include/utils/FaultTreeUtils.h b/include/utils/FaultTreeUtils.h new file mode 100644 index 0000000000..878b9591c0 --- /dev/null +++ b/include/utils/FaultTreeUtils.h @@ -0,0 +1,106 @@ +#ifndef _FAULT_TREE_H +#define _FAULT_TREE_H + +#include +#include +#include +#include + +#include "Utils.h" + +namespace FTAUtils +{ +class Parser; +class Quantification; +class FaultTree; +class CException; +std::string trim(const std::string & str); +std::string str2Upper(const std::string & str_in, bool trim_input = false); +double interpolate(std::vector> data, double x, bool extrapolate); +double normalCDF(double x); // Phi(-∞, x) aka N(x) +double normalCDF(double x, double mu, double sigma); +double lnCDF(double x, double mu, double sigma); +} // namespace FTAUtils + +// Fault Tree class +class FTAUtils::FaultTree +{ +public: + enum _operator_t + { + AND, + OR, + UNDEF + }; + + struct _node + { + std::string _name; + _operator_t _op; + std::vector _child; + _node(std::string name, _operator_t op) + { + this->_name = name; + this->_op = op; + } + }; + + FaultTree(std::string file_name, std::string root = ""); + FaultTree(std::set> & sets_link, std::map _node_base); + ~FaultTree(); + std::string getRoot(); + void printSets(); + static void printSets(std::set> sets); + static void printRow(std::set row); + std::set> getCutSets(); + std::map buildTree(Parser parser); + std::set> computeMinimumCutSets(); + std::vector event(std::map); + +private: + // Hash map for operators + std::map _opDict = {{"AND", AND}, {"OR", OR}}; + // Inverse mapping for printing purpose only + std::map<_operator_t, std::string> _opDictInv = {{AND, "AND"}, {OR, "OR"}}; + std::map _node_d_b; + std::string _root; + + // Cut sets container for in-place computations + std::set> _sets; + + // Member functions + _operator_t str2Operator(std::string op); + void rmSets(); + _node * getNode(std::string name); + void cutSetsExpand(_node * node); + void removeSubsets(); +}; + +class FTAUtils::Parser +{ +public: + // Supported formats for parsing + enum parseFormatT + { + FORMAT_CSV, + FORMAT_UNDEF + }; + + Parser(std::string fileName, parseFormatT format); + ~Parser(); + std::vector yieldLine(); + std::vector> yieldLines(); + +private: + // Handle to file + std::ifstream * fileP; +}; + +class FTAUtils::CException +{ +public: + std::string msg; + CException(std::string s) : msg(s) {} +}; + +#endif // _FAULT_TREE_H diff --git a/include/utils/QuantificationUtils.h b/include/utils/QuantificationUtils.h new file mode 100644 index 0000000000..7a004d36b7 --- /dev/null +++ b/include/utils/QuantificationUtils.h @@ -0,0 +1,189 @@ +#ifndef _QUANTIFICATION_H +#define _QUANTIFICATION_H + +#include +#include "FaultTreeUtils.h" +#include "MastodonUtils.h" + +#define CLIP(A, MIN, MAX) (((A) < MIN) ? MIN : ((A) > MAX ? MAX : (A))) + +#define GEN_QUALIFICATION_R_VEC(gen, d, n, rv) \ + for (int index = 0; index < n; index++) \ + { \ + double dataPoint = d(gen); \ + rv.push_back(CLIP(dataPoint, 0.0, 1.0)); \ + } + +// Quantification class +class FTAUtils::Quantification +{ +public: + enum _dist_t + { + PE, // = 0 + EXP, // = 1 + GAMMA, // = 2 + WEIBULL, + EXTREME, + NORM, + LNORM, + CHI_SQ, + CAUCHY, + FISHER_F, + STUDENT_T + }; + + enum _analysis_t + { + FRAGILITY, // = 0 + RISK // = 1 + }; + + Quantification(std::map>> & params_double, + std::map>> & params_string, + std::map & params_int, + std::map & params_bool, + std::map & params_analysis_t, + std::string events_file, + std::string events_prob_file, + _analysis_t analysis = FRAGILITY, + std::string hazard_file = "", + double im_lower = 0.1, + double im_upper = 4, + int n_bins = 15, + bool uncertainty = false, + std::string root = "", + int n_sample = 1, + int seed = 0); + +private: + struct Stats + { + double _pe; + double _mean; + double _median; + double _sd; + double _variance; + double _p5; + double _p95; + Stats(std::vector vec_in) + { + + std::vector vec; + vec.clear(); + + // Remove all INF from vector + for (int index = 0; index < vec_in.size(); index++) + { + if (!std::isinf(vec_in[index])) + { + vec.push_back(vec_in[index]); + } + } + + int size = vec.size(); + ASSERT(size > 0, "Size 0 vector passed for stats computation"); + + _pe = vec[0]; + + // TODO: Remove 0th element from further calculations + // vec.erase( vec.begin() ); + + sort(vec.begin(), vec.end()); + _mean = accumulate(vec.begin(), vec.end(), 0.0) / size; + + _variance = 0; + for (int i = 0; i < vec.size(); i++) + { + _variance += pow(vec[i] - _mean, 2); + } + _variance /= vec.size(); + _sd = sqrt(_variance); + + _median = (size % 2) ? vec[size / 2] : ((vec[size / 2 - 1] + vec[size / 2]) / 2); + + _p5 = getPercentile(vec, 5); + _p95 = getPercentile(vec, 95); + } + double getPercentile(std::vector vec, int percentile) + { + double p_i = (vec.size() + 1) * (percentile / 100.0); + int p_min = floor(p_i); + p_min = (p_min < 1) ? 1 : p_min; + int p_max = ceil(p_i); + p_max = (p_max > (vec.size() - 1)) ? vec.size() : p_max; + return vec[p_min - 1] + ((p_i - p_min) * (vec[p_max - 1] - vec[p_min - 1])); + } + void printStats() + { + std::cout << "Mean: " << _mean << ", Median: " << _median << ", SD: " << _sd + << ", 5th: " << _p5 << ", 95th: " << _p95 << "\n\n"; + } + }; + + // TODO: 3 of the following are not supported as they need 1 arg rather than 2 + // for their distribution. Remove these 3 if not needed else add support + std::map _str2dist = { + /*{"EXP" , EXP },*/ + {"GAMMA", GAMMA}, + {"WEIBULL", WEIBULL}, + {"EXTREME", EXTREME}, + {"NORM", NORM}, + {"LNORM", LNORM}, + {"PE", PE}, + /*{"CHI_SQ" , CHI_SQ },*/ + {"CAUCHY", CAUCHY}, + {"FISHER_F", FISHER_F}, + /*{"STUDENT_T", STUDENT_T}*/ + }; + + std::vector> _cut_set_prob; + std::map> _b_nodes; + std::vector getProbVector(_dist_t dist, + double a, + double b, + int n, + int seed, + std::vector im, + _analysis_t analysis, + bool uc); + std::vector> beProb(FTAUtils::Parser parser, + int n_sample, + int seed, + _analysis_t analysis, + std::vector intmes, + bool uncert); + + std::vector> computeCutSetProb(std::set> cut_sets, + int n, + bool bypass = false, + std::string bypass_key = "", + double bypass_value = 0); + std::vector cutSetRowProb(std::set row, + int n, + bool sign_positive = true, + bool bypass = false, + std::string bypass_key = "", + double bypass_value = 0); + std::vector * + cutSetProbWDigest(std::set> cut_sets, int n, bool only_min_max = false); + double getGateProb(double a, double b, bool is_and); + std::vector> minCutIM(std::vector upperBound); + std::map>> + beIM(std::set> cut_sets, + int n, + std::vector upper_bound, + std::vector & count_v); + std::vector> linesToDouble(std::vector> lines); + std::vector getBinMeans(double im_lower, double im_upper, int n_bins); + std::vector hazInterp(std::vector> hazard, + std::vector im_bins); + std::vector fragility(std::set> cut_sets, + int n, + std::vector im_bins, + double & mu, + double & sigma); + void computeRisk(int n, std::vector hazard); +}; + +#endif // _QUANTIFICATION_H diff --git a/include/utils/Utils.h b/include/utils/Utils.h new file mode 100644 index 0000000000..e58a3dde50 --- /dev/null +++ b/include/utils/Utils.h @@ -0,0 +1,19 @@ +#ifndef _UTILS_H +#define _UTILS_H + +#include +#include + +//-------------------------- MACRO UTILITIES BEGIN --------------------- +// A custom assert block for assertions +#define ASSERT(condition, statement, ...) \ + if (!(condition)) \ + { \ + fprintf(stderr, \ + "[ASSERT] In File: %s, Line: %d => " #statement "\n", \ + __FILE__, \ + __LINE__, \ + ##__VA_ARGS__); \ + abort(); \ + } +#endif // _UTILS_H diff --git a/src/utils/FaultTreeUtils.C b/src/utils/FaultTreeUtils.C new file mode 100644 index 0000000000..171b48850e --- /dev/null +++ b/src/utils/FaultTreeUtils.C @@ -0,0 +1,519 @@ +#include "FaultTreeUtils.h" +#include "MastodonUtils.h" +// #include "DelimitedFileReader.h" + +/* + * Constructor for fault tree class + * NOTE: Default value of root is the first node in file + */ +/*!public*/ +FTAUtils::FaultTree::FaultTree(std::string file_name, std::string root) +/*!endpublic*/ +{ + FTAUtils::Parser parser = FTAUtils::Parser(file_name, FTAUtils::Parser::FORMAT_CSV); + buildTree(parser); + + // Override root node if not default + if (root != "") + { + ASSERT(getNode(root), + "Root node requested (%s) not found in heirarchy. Please the check " + "file %s", + root.c_str(), + file_name.c_str()); + this->_root = root; + } + + computeMinimumCutSets(); +} + +/*!public*/ +FTAUtils::FaultTree::FaultTree(std::set> & sets_link, + std::map _node_base) +/*!endpublic*/ +{ + std::string root = ""; + _node_d_b = _node_base; + + // Clearing up sets vector for safety + rmSets(); + + // Add root to set + std::set root_set, sub_sets; + root_set.insert("TOP"); // "TOP", size: 1 + + _sets.insert(root_set); + + // Start with root node to expand on heirarchy + cutSetsExpand(getNode("TOP")); + + // Check if first row is empty + // NOTE: Since its std::set, only first row could + // be empty + if (_sets.begin()->size() == 0) + { + _sets.erase(_sets.begin()); + } + + removeSubsets(); + + sets_link = _sets; +} + +/* + * Destructor + */ +/*!public*/ +FTAUtils::FaultTree::~FaultTree() +/*!endpublic*/ {} + +/* + * Builds m-ary fault tree + */ +/*!public*/ +std::map +FTAUtils::FaultTree::buildTree(FTAUtils::Parser parser) +/*!endpublic*/ +{ + std::vector line; + while (true) + { + line = parser.yieldLine(); + + // Stop if no new line to process + if (line.size() == 0) + break; + + // Stash name, operator + _node * node = new _node(line[0], str2Operator(line[1])); + + // Add children + for (int i = 2; i < line.size(); i++) + node->_child.push_back(line[i]); + + // Stash the first entry as ROOT of tree + if (_node_d_b.size() == 0) + _root = line[0]; + + // Add the newly created node to node lookup hashmap + _node_d_b[line[0]] = node; + } + + return _node_d_b; +} + +/* + * Translates string to opeartor + */ +/*!private*/ +FTAUtils::FaultTree::_operator_t +FTAUtils::FaultTree::str2Operator(std::string op) +/*!endprivate*/ +{ + std::string op_s = FTAUtils::str2Upper(op, true); + ASSERT(_opDict.count(op_s) != 0, "Illegal Operator found: %s", op.c_str()); + + return _opDict[op_s]; +} + +/* + * Computes minimum cut sets based on MOCUS Algorithm + */ +/*!public*/ +std::set> +FTAUtils::FaultTree::computeMinimumCutSets() +/*!endpublic*/ +{ + // Clearing up sets vector for safety + rmSets(); + + // Add root to set + std::set root_set, sub_sets; + root_set.insert(getRoot()); + + _sets.insert(root_set); + + // Start with root node to expand on heirarchy + cutSetsExpand(getNode(getRoot())); + + // Check if first row is empty + // NOTE: Since its std::set, only first row could + // be empty + if (_sets.begin()->size() == 0) + { + _sets.erase(_sets.begin()); + } + + removeSubsets(); + + return _sets; +} + +/*!private*/ +void +FTAUtils::FaultTree::removeSubsets() +/*!endprivate*/ +{ + std::set rm_its; + + uint64_t index1 = 0; + for (std::set>::iterator row1_it = _sets.begin(); row1_it != _sets.end(); + ++row1_it) + { + uint64_t index2 = index1 + 1; + for (std::set>::iterator row2_it = next(row1_it, 1); + row2_it != _sets.end(); + ++row2_it) + { + bool row1_is_subset = + includes(row2_it->begin(), row2_it->end(), row1_it->begin(), row1_it->end()); + bool row2_is_subset = + includes(row1_it->begin(), row1_it->end(), row2_it->begin(), row2_it->end()); + // Remove row1 if row2 is a subset of row1 by marking it + if (row2_is_subset) + rm_its.insert(index1); + // Remove row2 if row1 is a subset of row1 by marking it + else if (row1_is_subset) + rm_its.insert(index2); + index2++; + } + index1++; + } + + uint64_t index = 0; + // Remove rows marked to be removed + for (std::set::iterator it = rm_its.begin(); it != rm_its.end(); ++it) + { + _sets.erase(next(_sets.begin(), *it - index)); + index++; + } +} + +/* + * Recursive call function to flood fill sets by expanding them based on + * operation ALGO: + * 1. Iterate through entire list and match for own node's name + * 2. Replace self with children based on operation + * (i) . Replace self with children in same row if AND + * (ii). Replace self with child one per row if OR + * 3. Recurse on each non leaf child + * + * NOTE: 1. Updates "sets" variable + * 2. Uses std::set 2d array, hence absorption and idempotence properties + * are implicit + */ +/*!private*/ +void +FTAUtils::FaultTree::cutSetsExpand(_node * node) +/*!endprivate*/ +{ + ASSERT(node, "NULL node encountered"); + ASSERT(node->_child.size() != 0, "Empty child list found for node '%s'", node->_name.c_str()); + + // New rows which have to be appended at end. Adding them might + // result in data hazards + std::set> mk_rows; + + // 1. Iterate through entire list and match for own node's name + for (std::set>::iterator row_it = _sets.begin(); row_it != _sets.end(); + ++row_it) + { + std::set row(*row_it); + // Search for the element in row + std::set::iterator it = find(row.begin(), row.end(), node->_name); + if (it != row.end()) + { + // Erase self to add children + row.erase(it); + + // 2. Replace self with children based on operation + // (i) . Replace self with children in same row if AND + // (ii). Replace self with child one per row if OR + switch (node->_op) + { + case AND: + // Append all children in same row + for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) + { + row.insert(node->_child[c_id]); + } + break; + + case OR: + // Replicate line and append one child per row + for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) + { + std::set row_set(row); + row_set.insert(node->_child[c_id]); + mk_rows.insert(row_set); + } + + // Clear this row to be removed at end (postprocessing) + row.clear(); + break; + + default: + ASSERT(false, "Unknown Operator found!"); + break; + } + } + + // The following 2 lines might result in an empty row to be pushed. + // If we gate insertion with count, we might have a data hazard due + // to wrong increment of iterators. + // A quick fix to this is to check for first row and delete it at end + _sets.erase(row_it); + _sets.insert(row); + } + + // Add newly created rows + for (std::set>::iterator it = mk_rows.begin(); it != mk_rows.end(); ++it) + { + std::set row(*it); + _sets.insert(*it); + } + + // 3. Recurse on each non leaf child + for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) + { + _node * child = getNode(node->_child[c_id]); + if (child) + cutSetsExpand(child); + } +} + +/*!private*/ +void +FTAUtils::FaultTree::rmSets() +/*!endprivate*/ +{ + for (std::set>::iterator row_it = _sets.begin(); row_it != _sets.end(); + ++row_it) + { + _sets.erase(row_it); + } +} + +/*!private*/ +FTAUtils::FaultTree::_node * +FTAUtils::FaultTree::getNode(std::string name) +/*!endprivate*/ { return (_node_d_b.count(name) != 0) ? _node_d_b[name] : NULL; } + +/*!public*/ +std::string +FTAUtils::FaultTree::getRoot() +/*!endpublic*/ { return _root; } + +/*!public*/ +void +FTAUtils::FaultTree::printSets(std::set> sets) +/*!endpublic*/ +{ + std::cout << "------------- SETS BEGIN ----------------- " << std::endl; + for (std::set>::iterator row = sets.begin(); row != sets.end(); ++row) + { + printRow(*row); + std::cout << std::endl; + } + std::cout << "-------------- SETS END ------------------ " << std::endl; +} + +/*!public*/ +void +FTAUtils::FaultTree::printSets() +/*!endpublic*/ { printSets(_sets); } + +/*!public*/ +void +FTAUtils::FaultTree::printRow(std::set row) +/*!endpublic*/ +{ + for (std::set::iterator col = row.begin(); col != row.end(); ++col) + { + std::cout << *col << ", "; + } +} + +/* + * Function returns cut sets at the given point + * NOTE: If MOCUS ran before this function call, min cut sets will be returned + */ +/*!public*/ +std::set> +FTAUtils::FaultTree::getCutSets() +/*!endpublic*/ { return _sets; } + +// **************************Parser Definition************************** +/* + * Constructor for parser class + */ +/*!public*/ +FTAUtils::Parser::Parser(std::string fileName, FTAUtils::Parser::parseFormatT format) +/*!endpublic*/ +{ + // Assertion on supported parsing formats + ASSERT(format == FORMAT_CSV, "Unsupported parse format"); + fileP = new std::ifstream; + fileP->open(fileName, std::ifstream::in); + // ASSERT( fileP->is_open(), "Unable to open file: %s", fileName.c_str()); + + if (!fileP->is_open()) + throw FTAUtils::CException("Unable to open file."); +} + +/* + * Destructor for parser class + */ +/*!public*/ +FTAUtils::Parser::~Parser() +/*!endpublic*/ {} + +/* + * Yields all records, populates the standard structure and returns + * This function acts as an abstract layer to hide different formats + * that might be supported in future + * + * Returns: Array of strings + */ +/*!public*/ +std::vector> +FTAUtils::Parser::yieldLines() +/*!endpublic*/ +{ + std::vector> lines; + std::vector line; + while (true) + { + line = FTAUtils::Parser::yieldLine(); + + // Stop if no new line to process + if (line.size() == 0) + break; + lines.push_back(line); + } + return lines; +} + +/* + * Yields a single record, populates the standard structure and returns + * This function acts as an abstract layer to hide different formats + * that might be supported in future + * + * Returns: Array of strings + */ +/*!public*/ +std::vector +FTAUtils::Parser::yieldLine() +/*!endpublic*/ +{ + ASSERT(fileP->is_open(), "Illegal call to yieldLine. File not open yet!"); + + std::string buffer; + std::vector line; + + // Get a line that has something except \n + if (getline(*fileP, buffer)) + { + std::string token; + std::istringstream iss(buffer); + while (getline(iss, token, ',')) + { + line.push_back(FTAUtils::trim(token)); + } + } + return line; +} + +/* + * String trim method which removes all leading and lagging whitespace in a string + */ +std::string +FTAUtils::trim(const std::string & str) +{ + std::string whiteSpace = " \t\n"; + + size_t first = str.find_first_not_of(whiteSpace); + if (std::string::npos == first) + { + return str; + } + size_t last = str.find_last_not_of(whiteSpace); + return str.substr(first, (last - first + 1)); +} + +// defined in Utils.cpp in FTA code + +/* + * Converts ASCII string to upper case + */ +std::string +FTAUtils::str2Upper(const std::string & str_in, bool trim_input) +{ + std::string str = str_in; + + if (trim_input) + str = trim(str_in); + + transform(str.begin(), str.end(), str.begin(), ::toupper); + return str; +} + +/* + * Returns interpolated value at x from parallel arrays ( xData, yData ) + * Assumes that xData has at least two elements, is sorted and is strictly monotonic increasing + * boolean argument extrapolate determines behaviour beyond ends of array (if needed) + */ +double +FTAUtils::interpolate(std::vector> data, double x, bool extrapolate) +{ + int size = data.size(); + + // find left end of interval for interpolation + int i = 0; + + // special case: beyond right end + if (x >= data[size - 2][0]) + { + i = size - 2; + } + else + { + while (x > data[i + 1][0]) + i++; + } + + // points on either side (unless beyond ends) + double xL = data[i][0], yL = data[i][1], xR = data[i + 1][0], yR = data[i + 1][1]; + + // if beyond ends of array and not extrapolating + if (!extrapolate) + { + if (x < xL) + yR = yL; + if (x > xR) + yL = yR; + } + + // gradient + double dydx = (yR - yL) / (xR - xL); + + // linear interpolation + return yL + dydx * (x - xL); +} + +// Phi(-∞, x) aka N(x) +double +FTAUtils::normalCDF(double x) +{ + return std::erfc(-x / std::sqrt(2)) / 2; +} + +double +FTAUtils::normalCDF(double x, double mu, double sigma) +{ + double y; + + y = (x - mu) / sigma; + + return FTAUtils::normalCDF(y); +} diff --git a/src/utils/QuantificationUtils.C b/src/utils/QuantificationUtils.C new file mode 100644 index 0000000000..5e9e29603a --- /dev/null +++ b/src/utils/QuantificationUtils.C @@ -0,0 +1,679 @@ +#include "QuantificationUtils.h" +// #include "DelimitedFileReader.h" + +// A macro to reduce typing while printing stats +#define QUANT_BE_STAT_PRINT_HELPER(index, i) \ + "PE: " << be_stats[index][i][0] << ", Mean: " << be_stats[index][i][1] \ + << ", Median: " << be_stats[index][i][2] << ", SD: " << be_stats[index][i][3] \ + << ", 5th: " << be_stats[index][i][4] << ", 95th: " << be_stats[index][i][5] << "\n" + +/* + * Constructor for qualifications class + */ +/*!public*/ +FTAUtils::Quantification::Quantification( + std::map>> & params_double, + std::map>> & params_string, + std::map & params_int, + std::map & params_bool, + std::map & params_analysis_t, + std::string events_file, + std::string events_prob_file, + _analysis_t analysis, + std::string hazard_file, + double im_lower, + double im_upper, + int n_bins, + bool uncertainty, + std::string root, + int n_sample, + int seed) +/*!endpublic*/ +{ + //---------------- SAVE PARAMETERS ------------------------------------------ + std::vector params_double1; + params_double1.push_back(im_lower); + params_double1.push_back(im_upper); + + std::vector> params_double2; + params_double2.push_back(params_double1); + params_double.insert( + std::pair>>("IM", params_double2)); + + params_int.insert(std::pair("n_bins", n_bins)); + params_int.insert(std::pair("nsamp", n_sample)); + params_int.insert(std::pair("seed", seed)); + params_bool.insert(std::pair("uncertainty", uncertainty)); + params_analysis_t.insert(std::pair("analysis", analysis)); + + //---------------- COMPUTE -------------------------------------------------- + // Construct the fault tree with MOCUS algo for minimal cut sets + FTAUtils::FaultTree ft = FTAUtils::FaultTree(events_file, root); + Parser parser_events = Parser(events_file, Parser::FORMAT_CSV); + std::vector> lines_events = parser_events.yieldLines(); + std::set> cut_sets = ft.getCutSets(); + // ft.printSets(); + params_string.insert( + std::pair>>("events_files", lines_events)); + + // Read basic events file + FTAUtils::Parser parser_event_prob = + FTAUtils::Parser(events_prob_file, FTAUtils::Parser::FORMAT_CSV); + std::vector> lines_events_prob; + + if (analysis == FRAGILITY) + { + // Read and interpolate hazard curve + FTAUtils::Parser parser_hazard = FTAUtils::Parser(hazard_file, FTAUtils::Parser::FORMAT_CSV); + std::vector> lines_hazard = parser_hazard.yieldLines(); + std::vector> hazard = FTAUtils::Quantification::linesToDouble(lines_hazard); + params_double.insert( + std::pair>>("hazard", hazard)); + + // List of Intensity Measure bin values + std::vector im_bins = getBinMeans(im_lower, im_upper, n_bins); + + // List of Intensity Measure bin values + std::vector hazard_freq = hazInterp(hazard, im_bins); + + // Dictionary of basic events for fragility input + lines_events_prob = beProb(parser_event_prob, n_sample, seed, analysis, im_bins, uncertainty); + + // Top event fragility (lognormal parameters) + double mu, sigma; + fragility(cut_sets, n_bins, im_bins, mu, sigma); + /* + std::cout << "----------- LN PARAMS BEGIN --------------" << std::endl; + std::cout << "mu: " << mu << std::endl; + std::cout << "sigma: " << sigma << std::endl; + std::cout << "------------ LN PARAMS END ---------------" << std::endl; + */ + + // Dictionary of basic events risk (convoluting fragility and hazard) + computeRisk(n_bins, hazard_freq); + } + else + { + // Dictionary of basic events risk (risk inputs) + lines_events_prob = + beProb(parser_event_prob, n_sample, seed, analysis, std::vector(), uncertainty); + /* + std::cout << "------------ LN PARAMS END ---------------" << std::endl; + */ + } + params_string.insert(std::pair>>( + "basic_events", lines_events_prob)); + + // Calculate minimal cut sets probability + // Digest cut set probability to unified probability + // Adding 1 to accomodate mean as 0th element + std::vector * digest = cutSetProbWDigest(cut_sets, n_sample + 1); + std::vector> mc_im = minCutIM(digest[1]); + Stats s0(digest[0]); + Stats s1(digest[1]); + Stats s2(digest[2]); + std::vector count_v; + std::map>> be_stats = + beIM(cut_sets, n_sample + 1, digest[1], count_v); + + //---------------- PRINTING -------------------------------------------------- + /* + std::cout << "---------- PROBABILITY BEGIN -------------" << std::endl; + std::cout << "1. Exact solution: " << s0._pe << std::endl; + s0.printStats(); + std::cout << "2. Upper Bound : " << s1._pe << std::endl; + s1.printStats(); + std::cout << "3. Rare Event : " << s2._pe << std::endl; + s2.printStats(); + */ + + // return the FTA top event risk + std::vector fta_0; + std::vector> fta_1; + fta_0.push_back(s0._pe); + fta_0.push_back(s1._pe); + fta_0.push_back(s2._pe); + fta_1.push_back(fta_0); + params_double.insert(std::make_pair("fta", fta_1)); + + /* + std::cout << "----------- PROBABILITY END --------------" << std::endl; + std::cout << "-------- CUT SET DETAILS BEGIN -----------" << std::endl; + + int i = 0; + for (std::set>::iterator row = cut_sets.begin(); row != cut_sets.end(); + ++row) { + ft.printRow(*row); + std::cout << "\t(PE: " << _cut_set_prob[i][0] << ", IM: " << mc_im[i][0] << ")" + << std::endl; + i++; + } + + std::cout << "--------- CUT SET DETAILS END ------------" << std::endl; + i = 0; + for (std::map>::iterator bn_it = _b_nodes.begin(); + bn_it != _b_nodes.end(); ++bn_it) { + std::cout << bn_it->first << " (" << count_v[i] << ")" << std::endl; + std::cout << "\t[FV] " QUANT_BE_STAT_PRINT_HELPER("fv", i); + std::cout << "\t[RRR] " QUANT_BE_STAT_PRINT_HELPER("rrr", i); + std::cout << "\t[RIR] " QUANT_BE_STAT_PRINT_HELPER("rir", i); + std::cout << "\t[RRI] " QUANT_BE_STAT_PRINT_HELPER("rri", i); + std::cout << "\t[RII] " QUANT_BE_STAT_PRINT_HELPER("rii", i); + std::cout << "\t[BI] " QUANT_BE_STAT_PRINT_HELPER("bi", i); + i++; + } + */ + + // Return the Risk Reduction Ratio for B1, B2, B3, B4, B5 + params_double.insert( + std::pair>>("fv", be_stats["fv"])); + params_double.insert( + std::pair>>("rrr", be_stats["rrr"])); + params_double.insert( + std::pair>>("rir", be_stats["rir"])); + + // Return the Risk Reduction Difference for B1, B2, B3, B4, B5 + params_double.insert( + std::pair>>("rri", be_stats["rri"])); + params_double.insert( + std::pair>>("rii", be_stats["rii"])); + params_double.insert( + std::pair>>("bi", be_stats["bi"])); +} + +/* + * Function for interpolating the hazard curve based on range of intensity + * measures + */ +/*!private*/ +std::vector +FTAUtils::Quantification::hazInterp(std::vector> hazard, + std::vector im_bins) +/*!endprivate*/ +{ + std::vector data; + // Convert all hazard values to log10 + for (int index = 0; index < hazard.size(); index++) + { + ASSERT(hazard[index].size() == 2, "hazard dimension error"); + hazard[index][0] = log10(hazard[index][0]); + hazard[index][1] = log10(hazard[index][1]); + } + + // Interpolate + for (int index = 0; index < im_bins.size(); index++) + { + data.push_back(pow(10, interpolate(hazard, log10(im_bins[index]), true))); + } + return data; +} + +/* + * Function for getting BIN mean values of intensity measure + */ +/*!private*/ +std::vector +FTAUtils::Quantification::getBinMeans(double im_lower, double im_upper, int n_bins) +/*!endprivate*/ +{ + std::vector bins; + double delta = (im_upper - im_lower) / (n_bins); + for (int index = 0; index < n_bins; index++) + { + double bin = im_lower + (delta * (index + 0.5)); + bins.push_back(bin); + } + return bins; +} + +/* + * Translator function for string -> double + */ +/*!private*/ +std::vector> +FTAUtils::Quantification::linesToDouble(std::vector> lines) +/*!endprivate*/ +{ + std::vector> lines_double; + for (int index = 0; index < lines.size(); index++) + { + std::vector double_vector(lines[index].size()); + transform(lines[index].begin(), + lines[index].end(), + double_vector.begin(), + [](const std::string & val) { return stod(val); }); + lines_double.push_back(double_vector); + } + return lines_double; +} + +/* + * Generates probability vector for a specified distribution + * Nomenclature of a and b changes with distribution. + * eg., a => mean, b => std for NORM + */ +/*!private*/ +std::vector +FTAUtils::Quantification::getProbVector(_dist_t dist, + double a, + double b, + int n, + int seed, + std::vector im, + _analysis_t analysis, + bool uc) +/*!endprivate*/ +{ + std::vector rv; + + std::default_random_engine gen(seed); + + if (analysis == FRAGILITY) + { + ASSERT(dist == LNORM, "unsupported distribution for fragility analysis"); + for (int index = 0; index < im.size(); index++) + { + rv.push_back(FTAUtils::normalCDF(log(im[index] / a) / b)); + } + } + else + { + if (!uc) + { + switch (dist) + { + case PE: + { + GEN_QUALIFICATION_R_VEC(a, double, 2, rv); + } + break; + case NORM: + { + GEN_QUALIFICATION_R_VEC(a, double, 2, rv); + } + break; + default: + ASSERT(false, "Un-supported dist found"); + } + } + else + { + switch (dist) + { + case PE: + { + GEN_QUALIFICATION_R_VEC(a, double, n, rv); + } + break; + case NORM: + { + rv.push_back(CLIP(a, 0, 1)); + std::normal_distribution d(a, b); + GEN_QUALIFICATION_R_VEC(gen, d, n, rv); + } + break; + default: + ASSERT(false, "Un-supported dist found"); + } + } + } + return rv; +} + +/* + * Parses and floods probabilties of basic elements + */ +/*!private*/ +std::vector> +FTAUtils::Quantification::beProb(FTAUtils::Parser parser, + int n_sample, + int seed, + _analysis_t analysis, + std::vector intmes, + bool uncert) +/*!endprivate*/ +{ + std::vector line; + std::vector> lines; + while (true) + { + line = parser.yieldLine(); + + // Stop if no new line to process + if (line.size() != 0) + lines.push_back(line); + else + break; + + // Stash name, probability vector + double b = line.size() > 3 ? stod(line[3]) : 0; + _b_nodes[line[0]] = getProbVector( + _str2dist[line[1]], stod(line[2]), b, n_sample, seed, intmes, analysis, uncert); + } + + return lines; +} + +/* + * Computes accumulated probability for the entire cut set + * 3 Digests are computed: + * 0. Min Max + * 1. Upper Bound + * 2. Top Rare + * NOTE: 1. Its better to compute these 2 together as they have common loops + * which can save time + * 2. Sets the vector cutSetProb + */ +/*!private*/ +std::vector * +FTAUtils::Quantification::cutSetProbWDigest(std::set> cut_sets, + int n, + bool only_min_max) +/*!endprivate*/ +{ + std::vector * digest = new std::vector[3]; + + // Digest 2: Generate power set to generate all possible combinations + std::vector> ps_p; + int pow_set_size = pow(2, cut_sets.size()); + // Start from 1 as NULL set is not needed + for (int index = 1; index < pow_set_size; index++) + { + std::set power_set_row; + for (int j = 0; j < cut_sets.size(); j++) + { + // Check if jth bit in the index is set + // If set then form row of power set + if (index & (1 << j)) + { + std::set j_elem = *next(cut_sets.begin(), j); + for (std::set::iterator it = j_elem.begin(); it != j_elem.end(); ++it) + { + power_set_row.insert(*it); + } + } + } + // For elements taken n at a time, if n is odd, positive else negative + bool is_positive = (__builtin_popcount(index) % 2); + std::vector cut_prob = cutSetRowProb(power_set_row, n, is_positive); + ps_p.push_back(cut_prob); + } + + // Digest 0, 1, 2 + // Avoiding an extra function call to interm event like in the python + // counterpart for performance reasons (loop merging) + _cut_set_prob = computeCutSetProb(cut_sets, n); + for (int index = 0; index < n; index++) + { + double min_max = 0; + double upper_bound = 1; + double top_rare = 0; + if (!only_min_max) + { + for (std::vector>::iterator it = _cut_set_prob.begin(); + it != _cut_set_prob.end(); + ++it) + { + top_rare += (*it)[index]; + upper_bound *= 1 - (*it)[index]; + } + digest[1].push_back(1 - upper_bound); + digest[2].push_back(top_rare); + } + for (std::vector>::iterator it = ps_p.begin(); it != ps_p.end(); ++it) + { + min_max += (*it)[index]; + } + digest[0].push_back(min_max); + } + + return digest; +} + +/* + * Computes probability for a given cut set on a per set basis based on + * pre-flooded basic elem probs + */ +/*!private*/ +std::vector> +FTAUtils::Quantification::computeCutSetProb(std::set> cut_sets, + int n, + bool bypass, + std::string bypass_key, + double bypass_value) +/*!endprivate*/ +{ + // For each row in cut set, compute: + // 1. For each sample in the generated vector, compute AND gate + // probability analysis for each basic element + std::vector> quant; + for (std::set>::iterator row = cut_sets.begin(); row != cut_sets.end(); + ++row) + { + quant.push_back(cutSetRowProb(*row, n, true, bypass, bypass_key, bypass_value)); + } + return quant; +} + +/*!private*/ +std::vector +FTAUtils::Quantification::cutSetRowProb(std::set row, + int n, + bool sign_positive, + bool bypass, + std::string bypass_key, + double bypass_value) +/*!endprivate*/ +{ + std::vector prob; + double sign_m = sign_positive ? 1 : -1; + for (int index = 0; index < n; index++) + { + double value = 1; + for (std::set::iterator col = row.begin(); col != row.end(); ++col) + { + ASSERT( + _b_nodes.find(*col) != _b_nodes.end(), "'%s' key not found in _b_nodes", (*col).c_str()); + + double bv = (bypass && (bypass_key == *col)) ? bypass_value : _b_nodes[*col][index]; + value = getGateProb(value, bv, true); + } + prob.push_back(value * sign_m); + } + return prob; +} + +/* + * Wrapper function for gate probability arith + * AND => a * b + * OR => (1 - a) * (1 - b) + */ +/*!private*/ +double +FTAUtils::Quantification::getGateProb(double a, double b, bool is_and) +/*!endprivate*/ { return is_and ? a * b : (1.0 - a) * (1.0 - b); } + +/* + * Computes cut-set Fussel-Vesely Important measures + */ +/*!private*/ +std::vector> +FTAUtils::Quantification::minCutIM(std::vector upper_bound) +/*!endprivate*/ +{ + std::vector> mc_i_m; + for (int row = 0; row < _cut_set_prob.size(); row++) + { + std::vector mc_i_m_row; + for (int col = 0; col < upper_bound.size(); col++) + { + mc_i_m_row.push_back((100.0 * _cut_set_prob[row][col]) / upper_bound[col]); + } + mc_i_m.push_back(mc_i_m_row); + } + return mc_i_m; +} + +/* + * Function for calculating risk by convoluting hazard and fragility + * Assign risk to basic event probabilites + * WARNING: This consumes _b_nodes and then overwrites it + */ +/*!private*/ +void +FTAUtils::Quantification::computeRisk(int n, std::vector hazard) +/*!endprivate*/ +{ + for (std::map>::iterator bn_it = _b_nodes.begin(); + bn_it != _b_nodes.end(); + ++bn_it) + { + double risk = 0; + for (int index = 0; index < n; index++) + { + risk += bn_it->second[index] * hazard[index]; + } + _b_nodes[bn_it->first] = getProbVector(PE, risk, 0, 1, 0, std::vector(), RISK, false); + } +} + +/* + * Function for top event fragility + */ +/*!private*/ +std::vector +FTAUtils::Quantification::fragility(std::set> cut_sets, + int n, + std::vector im_bins, + double & mu, + double & sigma) +/*!endprivate*/ +{ + // 1. Calculate TOP event fragility using min-max approach (exact) {digest[0]} + std::vector * digest = cutSetProbWDigest(cut_sets, n, true); + + // 2. lognormal parameters of TOP Event fragility + // solveLnParams(im_bins, digest[0], mu, sigma); + + // Use the maximizelikehood function instead of solver function + // Use SGD for optimization + std::vector loc_space = {0.01, 2}; + std::vector sca_space = {0.01, 1}; + std::vector max_values = MastodonUtils::maximizeLogLikelihood( + im_bins, digest[0], loc_space, sca_space, 1000, false, 1e-05, 0.0001, 1000, 1028); + mu = max_values[0]; + sigma = max_values[1]; + + // TODO: 3. TOP event Fragility plot + return digest[0]; +} + +/* For each basic element: + * 1. Calculate #occurrence in minimal cut set (beCount) + * 2. Store row index of all occurrence (beIndex) + */ +/*!private*/ +std::map>> +FTAUtils::Quantification::beIM(std::set> cut_sets, + int n, + std::vector upper_bound, + std::vector & count_v) +/*!endprivate*/ +{ + std::map>> stats; + for (std::map>::iterator bn_it = _b_nodes.begin(); + bn_it != _b_nodes.end(); + ++bn_it) + { + // Generate available vector to save computation on per index loop + std::vector available; + int count = 0; + for (std::set>::iterator cs_it = cut_sets.begin(); + cs_it != cut_sets.end(); + ++cs_it) + { + bool is_a = cs_it->find(bn_it->first) != cs_it->end(); + count += is_a; + available.push_back(is_a); + } + count_v.push_back(count); + if (count != 0) + { + std::vector fv, rrr, rir, rri, rii, bi; + + std::vector> mc_p1 = + computeCutSetProb(cut_sets, n, true, bn_it->first, 1); + std::vector> mc_p0 = + computeCutSetProb(cut_sets, n, true, bn_it->first, 0); + + // OR it with loop merging + for (int index = 0; index < n; index++) + { + double f0_val = 1; + double f1_val = 1; + double fi_val = 1; + for (int row = 0; row < mc_p0.size(); row++) + { + f0_val *= 1 - mc_p0[row][index]; + f1_val *= 1 - mc_p1[row][index]; + if (available[row]) + fi_val *= 1 - _cut_set_prob[row][index]; + } + f0_val = 1 - f0_val; + f1_val = 1 - f1_val; + fi_val = 1 - fi_val; + + fv.push_back(fi_val / upper_bound[index]); + rrr.push_back(upper_bound[index] / f0_val); + rri.push_back(upper_bound[index] - f0_val); + rir.push_back(f1_val / upper_bound[index]); + rii.push_back(f1_val - upper_bound[index]); + bi.push_back(f1_val - f0_val); + } + // Calculate stats for current row and stash + Stats fv_stats(fv); + Stats rrr_stats(rrr); + Stats rri_stats(rri); + Stats rir_stats(rir); + Stats rii_stats(rii); + Stats bi_stats(bi); + stats["fv"].push_back({fv_stats._pe, + fv_stats._mean, + fv_stats._median, + fv_stats._sd, + fv_stats._p5, + fv_stats._p95}); + stats["rrr"].push_back({rrr_stats._pe, + rrr_stats._mean, + rrr_stats._median, + rrr_stats._sd, + rrr_stats._p5, + rrr_stats._p95}); + stats["rri"].push_back({rri_stats._pe, + rri_stats._mean, + rri_stats._median, + rri_stats._sd, + rri_stats._p5, + rri_stats._p95}); + stats["rir"].push_back({rir_stats._pe, + rir_stats._mean, + rir_stats._median, + rir_stats._sd, + rir_stats._p5, + rir_stats._p95}); + stats["rii"].push_back({rii_stats._pe, + rii_stats._mean, + rii_stats._median, + rii_stats._sd, + rii_stats._p5, + rii_stats._p95}); + stats["bi"].push_back({bi_stats._pe, + bi_stats._mean, + bi_stats._median, + bi_stats._sd, + bi_stats._p5, + bi_stats._p95}); + } + } + return stats; +} diff --git a/unit/hazard.txt b/unit/hazard.txt new file mode 100644 index 0000000000..63941fbd1e --- /dev/null +++ b/unit/hazard.txt @@ -0,0 +1,6 @@ +0.0608,1.00E-02 +0.2124,1.00E-03 +0.4,1.00E-04 +0.629,1.00E-05 +0.9344,1.00E-06 +1.3055,1.00E-07 \ No newline at end of file diff --git a/unit/logic1.txt b/unit/logic1.txt new file mode 100644 index 0000000000..a9a55915ae --- /dev/null +++ b/unit/logic1.txt @@ -0,0 +1,5 @@ +TE,OR,IE3,IE4 +IE4,OR,C4 +IE3,OR,C3,IE2 +IE2,AND,C2,IE1 +IE1,OR,C1 \ No newline at end of file diff --git a/unit/logic1_bas_events_LNORM.txt b/unit/logic1_bas_events_LNORM.txt new file mode 100644 index 0000000000..125ca2897a --- /dev/null +++ b/unit/logic1_bas_events_LNORM.txt @@ -0,0 +1,4 @@ +C1,LNORM,1.88,0.5 +C2,LNORM,3.78,0.79 +C3,LNORM,2.33,0.76 +C4,LNORM,3.66,0.45 \ No newline at end of file diff --git a/unit/logic2.txt b/unit/logic2.txt new file mode 100644 index 0000000000..ff88b3a821 --- /dev/null +++ b/unit/logic2.txt @@ -0,0 +1,7 @@ +TOP, AND, GATE1, GATE2 +GATE1, OR, FT-N/m-1, FT-N/m-2, FT-N/m-3 +GATE2, OR, B1, B3, B4 +GATE3, OR, B2, B4 +FT-N/m-1, AND, GATE3, B3, B5 +FT-N/m-2, AND, GATE3, B1 +FT-N/m-3, AND, B3, B5, B1 \ No newline at end of file diff --git a/unit/logic2_bas_events_PE.txt b/unit/logic2_bas_events_PE.txt new file mode 100644 index 0000000000..87c9c51b2c --- /dev/null +++ b/unit/logic2_bas_events_PE.txt @@ -0,0 +1,5 @@ +B1,PE,0.01 +B2,PE,0.02 +B3,PE,0.03 +B4,PE,0.04 +B5,PE,0.05 \ No newline at end of file diff --git a/unit/src/TestFaultTreeUtils.C b/unit/src/TestFaultTreeUtils.C new file mode 100644 index 0000000000..2c909cc95e --- /dev/null +++ b/unit/src/TestFaultTreeUtils.C @@ -0,0 +1,176 @@ +// MOOSE includes +#include "gtest/gtest.h" +#include "MooseUtils.h" +#include "Conversion.h" + +// Boost distribution includes +#include "BoostDistribution.h" + +// MASTODON includes +#include "MastodonUtils.h" + +// QUANTIFICATIONUTILS includes +#include "FaultTreeUtils.h" + +// Helper to make sure the test tree is correct. +void +assertTree(std::map tree_node) +{ + // TOP + EXPECT_EQ(tree_node["TOP"]->_name, "TOP"); + EXPECT_EQ(tree_node["TOP"]->_op, FTAUtils::FaultTree::AND); + EXPECT_EQ(tree_node["TOP"]->_child.at(0), "GATE1"); + EXPECT_EQ(tree_node["TOP"]->_child.at(1), "GATE2"); + + // GATE1 + EXPECT_EQ(tree_node["GATE1"]->_name, "GATE1"); + EXPECT_EQ(tree_node["GATE1"]->_op, FTAUtils::FaultTree::OR); + EXPECT_EQ(tree_node["GATE1"]->_child.at(0), "FT-N/m-1"); + EXPECT_EQ(tree_node["GATE1"]->_child.at(1), "FT-N/m-2"); + EXPECT_EQ(tree_node["GATE1"]->_child.at(2), "FT-N/m-3"); + + // GATE2 + EXPECT_EQ(tree_node["GATE2"]->_name, "GATE2"); + EXPECT_EQ(tree_node["GATE2"]->_op, FTAUtils::FaultTree::OR); + EXPECT_EQ(tree_node["GATE2"]->_child.at(0), "B1"); + EXPECT_EQ(tree_node["GATE2"]->_child.at(1), "B3"); + EXPECT_EQ(tree_node["GATE2"]->_child.at(2), "B4"); + + // FT-N/m-1 + EXPECT_EQ(tree_node["FT-N/m-1"]->_name, "FT-N/m-1"); + EXPECT_EQ(tree_node["FT-N/m-1"]->_op, FTAUtils::FaultTree::AND); + EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(0), "GATE3"); + EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(1), "B3"); + EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(2), "B5"); + + // FT-N/m-2 + EXPECT_EQ(tree_node["FT-N/m-2"]->_name, "FT-N/m-2"); + EXPECT_EQ(tree_node["FT-N/m-2"]->_op, FTAUtils::FaultTree::AND); + EXPECT_EQ(tree_node["FT-N/m-2"]->_child.at(0), "GATE3"); + EXPECT_EQ(tree_node["FT-N/m-2"]->_child.at(1), "B1"); + + // FT-N/m-3 + EXPECT_EQ(tree_node["FT-N/m-3"]->_name, "FT-N/m-3"); + EXPECT_EQ(tree_node["FT-N/m-3"]->_op, FTAUtils::FaultTree::AND); + EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(0), "B3"); + EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(1), "B5"); + EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(2), "B1"); + + // GATE3 + EXPECT_EQ(tree_node["GATE3"]->_name, "GATE3"); + EXPECT_EQ(tree_node["GATE3"]->_op, FTAUtils::FaultTree::OR); + EXPECT_EQ(tree_node["GATE3"]->_child.at(0), "B2"); + EXPECT_EQ(tree_node["GATE3"]->_child.at(1), "B4"); +} + +// Helper to make sure the mocus is correct. +void +assertMocus(std::set> tree_mocus) +{ + + std::set>::iterator it_tree_mocus; + std::set tree_mocus_children, set_mocus, set_gold; + std::set::iterator it_tree_mocus_children; + + std::set> matrix_gold = { + {"B2", "B1"}, {"B4", "B1"}, {"B1", "B5", "B3"}, {"B5", "B2", "B3"}, {"B4", "B5", "B3"}}; + + std::set>::iterator it_matrix_gold = matrix_gold.begin(); + + EXPECT_EQ(tree_mocus.size(), matrix_gold.size()); + + it_tree_mocus = tree_mocus.begin(); + int i = 0; + while (it_tree_mocus != tree_mocus.end()) + { + set_mocus = *it_tree_mocus++; + set_gold = *it_matrix_gold++; + EXPECT_EQ(set_mocus, set_gold); + i++; + } +} + +// Create a test event list. +std::map +eventList() +{ + + std::map _node_d_b; + + std::vector> line = {{"TOP", "AND", "GATE1", "GATE2"}, + {"GATE1", "OR", "FT-N/m-1", "FT-N/m-2", "FT-N/m-3"}, + {"GATE2", "OR", "B1", "B3", "B4"}, + {"GATE3", "OR", "B2", "B4"}, + {"FT-N/m-1", "AND", "GATE3", "B3", "B5"}, + {"FT-N/m-2", "AND", "GATE3", "B1"}, + {"FT-N/m-3", "AND", "B3", "B5", "B1"}}; + + std::vector op = {FTAUtils::FaultTree::AND, + FTAUtils::FaultTree::OR, + FTAUtils::FaultTree::OR, + FTAUtils::FaultTree::OR, + FTAUtils::FaultTree::AND, + FTAUtils::FaultTree::AND, + FTAUtils::FaultTree::AND}; + + for (int i = 0; i < line.size(); i++) + { + FTAUtils::FaultTree::_node * node = new FTAUtils::FaultTree::_node(line[i][0], op[i]); + // Add children + for (int j = 2; j < line[i].size(); j++) + node->_child.push_back(line[i][j]); + // Add the newly created node to node lookup hashmap + _node_d_b[line[i][0]] = node; + } + + return _node_d_b; +} + +// For input +std::vector file_lists_fault_tree_1, file_lists_fault_tree_2, file_lists_fault_tree_3; + +TEST(FTAUtils, FaultTree) +{ + // ==================== TestCase for FaultTree Object =================== + + // ---------------- Test error message for events input ---------------- + + // ###### File Inputs ###### + + try + { + file_lists_fault_tree_1 = {"not_a_valid_filename.txt", ""}; + + FTAUtils::FaultTree(file_lists_fault_tree_1[0], file_lists_fault_tree_1[1]); + } + catch (FTAUtils::CException e) + { + /* + * ------- IO Error ------- + * does not exist + */ + EXPECT_EQ(e.msg, "Unable to open file."); + } + + // ---------------- Test creating tree from list ---------------- + + std::map tree_node_from_list = eventList(); + std::set> tree_mocus_from_list; + FTAUtils::FaultTree ft_from_list = FTAUtils::FaultTree(tree_mocus_from_list, tree_node_from_list); + assertTree(tree_node_from_list); + assertMocus(tree_mocus_from_list); + + // ---------------- Test creating tree from file ---------------- + + file_lists_fault_tree_2 = {"logic2.txt", ""}; + + FTAUtils::FaultTree ft_from_file_2 = + FTAUtils::FaultTree(file_lists_fault_tree_2[0], file_lists_fault_tree_2[1]); + FTAUtils::Parser parser_events_2 = + FTAUtils::Parser(file_lists_fault_tree_2[0], FTAUtils::Parser::FORMAT_CSV); + std::map tree_node_from_file = + ft_from_file_2.buildTree(parser_events_2); + std::set> tree_mocus_from_file = ft_from_file_2.computeMinimumCutSets(); + assertTree(tree_node_from_file); + assertMocus(tree_mocus_from_file); +} diff --git a/unit/src/TestQuantificationUtils.C b/unit/src/TestQuantificationUtils.C new file mode 100644 index 0000000000..37037e9c0a --- /dev/null +++ b/unit/src/TestQuantificationUtils.C @@ -0,0 +1,456 @@ +// MOOSE includes +#include "gtest/gtest.h" +#include "MooseUtils.h" +#include "Conversion.h" + +// Boost distribution includes +#include "BoostDistribution.h" + +// MASTODON includes +#include "MastodonUtils.h" + +// QUANTIFICATIONUTILS includes +#include "QuantificationUtils.h" + +// For input +std::vector file_lists_quantification_utils; +std::vector IM; +int nbins; +bool uncertainty; +int nsamp; +int seed; + +// For output +std::map>> params_double_1, params_double_2, + params_double_3, params_double_4, params_double_5, params_double_6, params_double_7, + params_double_8, params_double_9; + +std::map>> params_string_1, params_string_2, + params_string_3, params_string_4, params_string_5, params_string_6, params_string_7, + params_string_8, params_string_9; + +std::map params_bool_1, params_bool_2, params_bool_3, params_bool_4, + params_bool_5, params_bool_6, params_bool_7, params_bool_8, params_bool_9; + +std::map params_analysis_t_1, + params_analysis_t_2, params_analysis_t_3, params_analysis_t_4, params_analysis_t_5, + params_analysis_t_6, params_analysis_t_7, params_analysis_t_8, params_analysis_t_9; + +std::map params_int_1, params_int_2, params_int_3, params_int_4, params_int_5, + params_int_6, params_int_7, params_int_8, params_int_9; + +// TestCase for Quantification Object +TEST(FTAUtils, Quantification) +{ + + // =========================================================================== + // =============================== Test Inputs =============================== + // =========================================================================== + + // +++++++++++++++++++++++ Input values +++++++++++++++++++++++ + + // ############### File Inputs ############### + + try + { + file_lists_quantification_utils = { + "not_a_valid_filename.txt", "not_a_valid_filename.txt", "not_a_valid_filename.txt"}; + IM = {0.1, 4}; + nbins = 15; + + FTAUtils::Quantification(params_double_1, + params_string_1, + params_int_1, + params_bool_1, + params_analysis_t_1, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins); + } + catch (FTAUtils::CException e) + { + /* + * ------- IO Error ------- + * does not exist + */ + EXPECT_EQ(e.msg, "Unable to open file."); + } + + try + { + file_lists_quantification_utils = {"", "", ""}; + IM = {0.1, 4}; + nbins = 15; + + FTAUtils::Quantification(params_double_2, + params_string_2, + params_int_2, + params_bool_2, + params_analysis_t_2, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins); + } + catch (FTAUtils::CException e) + { + /* + * ------- Type Error ------- + * must be a filename or a list + */ + EXPECT_EQ(e.msg, "Unable to open file."); + } + + // ############## Inputs for Fragility ############## + + file_lists_quantification_utils = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; + IM = {0.1, 4}; + nbins = 15; + FTAUtils::Quantification(params_double_3, + params_string_3, + params_int_3, + params_bool_3, + params_analysis_t_3, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::FRAGILITY, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins); + + // >>>>>>>> logic Value Check <<<<<<<< + + std::vector> event_files_matrix3{{"TE", "OR", "IE3", "IE4"}, + {"IE4", "OR", "C4"}, + {"IE3", "OR", "C3", "IE2"}, + {"IE2", "AND", "C2", "IE1"}, + {"IE1", "OR", "C1"}}; + EXPECT_EQ(params_string_3["events_files"], event_files_matrix3); + + // >>>>>>>> Basic Events Value Check <<<<<<<< + + std::vector> basic_events_matrix3{{"C1", "LNORM", "1.88", "0.5"}, + {"C2", "LNORM", "3.78", "0.79"}, + {"C3", "LNORM", "2.33", "0.76"}, + {"C4", "LNORM", "3.66", "0.45"}}; + std::vector> basic_events_3 = params_string_3["basic_events"]; + EXPECT_EQ(basic_events_3, basic_events_matrix3); + + // >>>>>>>> antype Value Check <<<<<<<< + + EXPECT_EQ(params_analysis_t_3["analysis"], FTAUtils::Quantification::FRAGILITY); + + // >>>>>>>> hazard Value Check <<<<<<<< + + std::vector> matrix_hazard{{0.0608, 0.01}, + {0.2124, 0.001}, + {0.4, 0.0001}, + {0.629, 1e-05}, + {0.9344, 1e-06}, + {1.3055, 1e-07}}; + std::vector> hazard = params_double_3["hazard"]; + EXPECT_EQ(hazard, matrix_hazard); + + // >>>>>>>> imrange Value Check <<<<<<<< + + EXPECT_EQ(params_double_3["IM"][0][0], 0.1); + EXPECT_EQ(params_double_3["IM"][0][1], 4); + + // >>>>>>>> nbins Value Check <<<<<<<< + + EXPECT_EQ(params_int_3["n_bins"], 15); + + // >>>>>>>> nbins Value Check <<<<<<<< + + try + { + file_lists_quantification_utils = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; + IM = {0.1, 4}; + + nbins = -15; + if (!(nbins > 0)) + throw FTAUtils::CException("ValueError"); + + FTAUtils::Quantification(params_double_4, + params_string_4, + params_int_4, + params_bool_4, + params_analysis_t_4, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::FRAGILITY, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins); + } + catch (FTAUtils::CException e) + { + /* + * ------- Value Error ------- + * The supplied value of nbins must be a +ve integer + */ + EXPECT_EQ(e.msg, "ValueError"); + } + + // ############## Inputs for Risk analysis (not fragility) ############## + + file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt"}; + FTAUtils::Quantification(params_double_5, + params_string_5, + params_int_5, + params_bool_5, + params_analysis_t_5, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK); + + // >>>>>>>> logic Value Check <<<<<<<< + + std::vector> event_files_matrix5{ + {"TOP", "AND", "GATE1", "GATE2"}, + {"GATE1", "OR", "FT-N/m-1", "FT-N/m-2", "FT-N/m-3"}, + {"GATE2", "OR", "B1", "B3", "B4"}, + {"GATE3", "OR", "B2", "B4"}, + {"FT-N/m-1", "AND", "GATE3", "B3", "B5"}, + {"FT-N/m-2", "AND", "GATE3", "B1"}, + {"FT-N/m-3", "AND", "B3", "B5", "B1"}}; + + EXPECT_EQ(params_string_5["events_files"], event_files_matrix5); + + // >>>>>>>> Basic Events Value Check <<<<<<<< + + std::vector> basic_events_matrix5{{"B1", "PE", "0.01"}, + {"B2", "PE", "0.02"}, + {"B3", "PE", "0.03"}, + {"B4", "PE", "0.04"}, + {"B5", "PE", "0.05"}}; + std::vector> basic_events_5 = params_string_5["basic_events"]; + EXPECT_EQ(basic_events_5, basic_events_matrix5); + + // >>>>>>>> antype Value Check <<<<<<<< + + EXPECT_NE(params_analysis_t_5["analysis"], FTAUtils::Quantification::FRAGILITY); + + // >>>>>>>> uncertainty Value Check <<<<<<<< + + EXPECT_EQ(params_bool_5["uncertainty"], false); + + // >>>>>>>> nsamp Value Check <<<<<<<< + + EXPECT_EQ(params_int_5["nsamp"], 1); + + // >>>>>>>> seed Value Check <<<<<<<< + + EXPECT_EQ(params_int_5["seed"], 0); + + // ###### Testing for input errors making sure parameters are input correctly ###### + + file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + IM = {0.1, 4}; + nbins = 15; + uncertainty = true; + nsamp = 1000; + seed = 436546754; + + FTAUtils::Quantification(params_double_6, + params_string_6, + params_int_6, + params_bool_6, + params_analysis_t_6, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins, + uncertainty, + "", + nsamp, + seed); + + // >>>>>>>> uncertainty Value Check <<<<<<<< + + EXPECT_EQ(params_bool_6["uncertainty"], true); + + // >>>>>>>> nsamp Value Check <<<<<<<< + + EXPECT_EQ(params_int_6["nsamp"], 1000); + + // >>>>>>>> seed Value Check <<<<<<<< + + EXPECT_EQ(params_int_6["seed"], 436546754); + + // ############## Testing for input type errors and value errors ############## + + // >>>>>>>> nsamp Value Check <<<<<<<< + + try + { + file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + IM = {0.1, 4}; + nbins = 15; + uncertainty = true; + + nsamp = -10; + if (!(nsamp > 0)) + throw FTAUtils::CException("ValueError"); + + seed = 42.0; + + FTAUtils::Quantification(params_double_7, + params_string_7, + params_int_7, + params_bool_7, + params_analysis_t_7, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins, + uncertainty, + "", + nsamp, + seed); + } + catch (FTAUtils::CException e) + { + /* + * ------- Value Error ------- + * The supplied value of nsamp must be a +ve integer. + */ + EXPECT_EQ(e.msg, "ValueError"); + } + + // >>>>>>>> seed Value Check <<<<<<<< + + try + { + file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + IM = {0.1, 4}; + nbins = 15; + uncertainty = true; + nsamp = 10; + + seed = -42; + if (!(seed > 0)) + throw FTAUtils::CException("ValueError"); + + FTAUtils::Quantification(params_double_8, + params_string_8, + params_int_8, + params_bool_8, + params_analysis_t_8, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins, + uncertainty, + "", + nsamp, + seed); + } + catch (FTAUtils::CException e) + { + /* + * ------- Value Error ------- + * The supplied value of seed must be a +ve integer. + */ + EXPECT_EQ(e.msg, "ValueError"); + } + + // =========================================================================== + // ============================== Test TOP Risk ============================== + // =========================================================================== + + // ############ Function for asserting FTA top event risk ############ + + // ++++++++++ Check FTA top event risk ++++++++++ + + file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + IM = {0.1, 4}; + nbins = 15; + FTAUtils::Quantification(params_double_9, + params_string_9, + params_int_9, + params_bool_9, + params_analysis_t_9, + file_lists_quantification_utils[0], + file_lists_quantification_utils[1], + FTAUtils::Quantification::RISK, + file_lists_quantification_utils[2], + IM[0], + IM[1], + nbins); + + // >>>>>>>> Risk Value Check <<<<<<<< + + // min-max + EXPECT_EQ(std::to_string(params_double_9["fta"][0][0]), "0.000694"); + + // upper bound + EXPECT_EQ(std::to_string(params_double_9["fta"][0][1]), "0.000705"); + + // rare event + EXPECT_EQ(std::to_string(params_double_9["fta"][0][2]), "0.000705"); + + // >>>>>>>> IMratio Value Check <<<<<<<< + + // Fussell-Vesely for B1, B2, B3, B4, B5 + EXPECT_EQ(std::to_string(params_double_9["fv"][0][0]), "0.872395"); + EXPECT_EQ(std::to_string(params_double_9["fv"][1][0]), "0.326300"); + EXPECT_EQ(std::to_string(params_double_9["fv"][2][0]), "0.148963"); + EXPECT_EQ(std::to_string(params_double_9["fv"][3][0]), "0.652584"); + EXPECT_EQ(std::to_string(params_double_9["fv"][4][0]), "0.148963"); + + // Risk Reduction Ratio for B1, B2, B3, B4, B5 + EXPECT_EQ(std::to_string(params_double_9["rrr"][0][0]), "7.831866"); + EXPECT_EQ(std::to_string(params_double_9["rrr"][1][0]), "1.483999"); + EXPECT_EQ(std::to_string(params_double_9["rrr"][2][0]), "1.174913"); + EXPECT_EQ(std::to_string(params_double_9["rrr"][3][0]), "2.877066"); + EXPECT_EQ(std::to_string(params_double_9["rrr"][4][0]), "1.174913"); + + // Risk Increase Ratio for B1, B2, B3, B4, B5 + EXPECT_EQ(std::to_string(params_double_9["rir"][0][0]), "86.111103"); + EXPECT_EQ(std::to_string(params_double_9["rir"][1][0]), "16.960273"); + EXPECT_EQ(std::to_string(params_double_9["rir"][2][0]), "5.808755"); + EXPECT_EQ(std::to_string(params_double_9["rir"][3][0]), "16.637742"); + EXPECT_EQ(std::to_string(params_double_9["rir"][4][0]), "3.826894"); + + // >>>>>>>> IMdiff Value Check <<<<<<<< + + // Risk Reduction Difference for B1, B2, B3, B4, B5 + EXPECT_EQ(std::to_string(params_double_9["rri"][0][0]), "0.000615"); + EXPECT_EQ(std::to_string(params_double_9["rri"][1][0]), "0.000230"); + EXPECT_EQ(std::to_string(params_double_9["rri"][2][0]), "0.000105"); + EXPECT_EQ(std::to_string(params_double_9["rri"][3][0]), "0.000460"); + EXPECT_EQ(std::to_string(params_double_9["rri"][4][0]), "0.000105"); + + // Risk Increase Difference for B1, B2, B3, B4, B5 + EXPECT_EQ(std::to_string(params_double_9["rii"][0][0]), "0.059991"); + EXPECT_EQ(std::to_string(params_double_9["rii"][1][0]), "0.011250"); + EXPECT_EQ(std::to_string(params_double_9["rii"][2][0]), "0.003389"); + EXPECT_EQ(std::to_string(params_double_9["rii"][3][0]), "0.011022"); + EXPECT_EQ(std::to_string(params_double_9["rii"][4][0]), "0.001993"); + + // Birnbaum Difference for B1, B2, B3, B4, B5 + EXPECT_EQ(std::to_string(params_double_9["bi"][0][0]), "0.060606"); + EXPECT_EQ(std::to_string(params_double_9["bi"][1][0]), "0.011480"); + EXPECT_EQ(std::to_string(params_double_9["bi"][2][0]), "0.003494"); + EXPECT_EQ(std::to_string(params_double_9["bi"][3][0]), "0.011482"); + EXPECT_EQ(std::to_string(params_double_9["bi"][4][0]), "0.002097"); +} From e525933c5361ce9f0c8d6aeded3556b313a789ef Mon Sep 17 00:00:00 2001 From: HugoNip Date: Fri, 27 Mar 2020 10:17:05 -0400 Subject: [PATCH 02/10] remove the banned keywords --- .vscode/settings.json | 8 ++++++++ include/utils/FaultTreeUtils.h | 6 +++--- include/utils/QuantificationUtils.h | 2 ++ src/utils/FaultTreeUtils.C | 18 +++++++++++------- 4 files changed, 24 insertions(+), 10 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000000..abcd643121 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,8 @@ +{ + "files.associations": { + "QuantificationUtils.C": "cpp", + "FaultTreeUtils.C": "cpp", + "TestMastodonUtils.C": "cpp", + "TestISoilUtils.C": "cpp" + } +} \ No newline at end of file diff --git a/include/utils/FaultTreeUtils.h b/include/utils/FaultTreeUtils.h index 878b9591c0..df8f12a62a 100644 --- a/include/utils/FaultTreeUtils.h +++ b/include/utils/FaultTreeUtils.h @@ -49,9 +49,9 @@ class FTAUtils::FaultTree FaultTree(std::set> & sets_link, std::map _node_base); ~FaultTree(); std::string getRoot(); - void printSets(); - static void printSets(std::set> sets); - static void printRow(std::set row); + // void printSets(); + // static void printSets(std::set> sets); + // static void printRow(std::set row); std::set> getCutSets(); std::map buildTree(Parser parser); std::set> computeMinimumCutSets(); diff --git a/include/utils/QuantificationUtils.h b/include/utils/QuantificationUtils.h index 7a004d36b7..cad5b8f0a4 100644 --- a/include/utils/QuantificationUtils.h +++ b/include/utils/QuantificationUtils.h @@ -114,11 +114,13 @@ class FTAUtils::Quantification p_max = (p_max > (vec.size() - 1)) ? vec.size() : p_max; return vec[p_min - 1] + ((p_i - p_min) * (vec[p_max - 1] - vec[p_min - 1])); } + /* void printStats() { std::cout << "Mean: " << _mean << ", Median: " << _median << ", SD: " << _sd << ", 5th: " << _p5 << ", 95th: " << _p95 << "\n\n"; } + */ }; // TODO: 3 of the following are not supported as they need 1 arg rather than 2 diff --git a/src/utils/FaultTreeUtils.C b/src/utils/FaultTreeUtils.C index 171b48850e..990d23623d 100644 --- a/src/utils/FaultTreeUtils.C +++ b/src/utils/FaultTreeUtils.C @@ -305,9 +305,9 @@ FTAUtils::FaultTree::getRoot() /*!endpublic*/ { return _root; } /*!public*/ -void -FTAUtils::FaultTree::printSets(std::set> sets) +// void FTAUtils::FaultTree::printSets(std::set> sets) /*!endpublic*/ +/* { std::cout << "------------- SETS BEGIN ----------------- " << std::endl; for (std::set>::iterator row = sets.begin(); row != sets.end(); ++row) @@ -317,22 +317,26 @@ FTAUtils::FaultTree::printSets(std::set> sets) } std::cout << "-------------- SETS END ------------------ " << std::endl; } +*/ /*!public*/ -void -FTAUtils::FaultTree::printSets() -/*!endpublic*/ { printSets(_sets); } +// void +// FTAUtils::FaultTree::printSets() +/*!endpublic*/ +// { printSets(_sets); } /*!public*/ -void -FTAUtils::FaultTree::printRow(std::set row) +// void +// FTAUtils::FaultTree::printRow(std::set row) /*!endpublic*/ +/* { for (std::set::iterator col = row.begin(); col != row.end(); ++col) { std::cout << *col << ", "; } } +*/ /* * Function returns cut sets at the given point From 48849343bebd37f151e9970715c4e767b75e7ac0 Mon Sep 17 00:00:00 2001 From: HugoNip Date: Sat, 28 Mar 2020 22:36:54 -0400 Subject: [PATCH 03/10] 1. remove json file 2. remove Utils.h 3. make macro #define a method 4. remove all commented sections 5. remove banned keywords --- .vscode/settings.json | 8 - include/utils/FaultTreeUtils.h | 114 ++++++++- include/utils/QuantificationUtils.h | 148 +++++++---- include/utils/Utils.h | 19 -- src/utils/FaultTreeUtils.C | 275 ++++++++------------ src/utils/QuantificationUtils.C | 381 ++++++++++------------------ unit/src/TestFaultTreeUtils.C | 40 +-- unit/src/TestQuantificationUtils.C | 75 +++--- 8 files changed, 508 insertions(+), 552 deletions(-) delete mode 100644 .vscode/settings.json delete mode 100644 include/utils/Utils.h diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index abcd643121..0000000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,8 +0,0 @@ -{ - "files.associations": { - "QuantificationUtils.C": "cpp", - "FaultTreeUtils.C": "cpp", - "TestMastodonUtils.C": "cpp", - "TestISoilUtils.C": "cpp" - } -} \ No newline at end of file diff --git a/include/utils/FaultTreeUtils.h b/include/utils/FaultTreeUtils.h index df8f12a62a..57df9b3a81 100644 --- a/include/utils/FaultTreeUtils.h +++ b/include/utils/FaultTreeUtils.h @@ -6,7 +6,10 @@ #include #include -#include "Utils.h" +typedef std::vector> vector_double; +typedef std::vector> vector_string; +typedef std::set> set_string; +typedef std::pair pair_string_double; namespace FTAUtils { @@ -14,15 +17,37 @@ class Parser; class Quantification; class FaultTree; class CException; + +/** + * String trim method which removes all leading and lagging whitespace in a string + */ std::string trim(const std::string & str); + +/** + * Converts ASCII string to upper case + */ std::string str2Upper(const std::string & str_in, bool trim_input = false); -double interpolate(std::vector> data, double x, bool extrapolate); -double normalCDF(double x); // Phi(-∞, x) aka N(x) + +/** + * Returns interpolated value at x from parallel arrays ( xData, yData ) + * Assumes that xData has at least two elements, is sorted and is strictly monotonic increasing + * boolean argument extrapolate determines behaviour beyond ends of array (if needed) + */ +double interpolate(vector_double data, double x, bool extrapolate); + +/** + * Phi(-∞, x) aka N(x) + */ +double normalCDF(double x); double normalCDF(double x, double mu, double sigma); -double lnCDF(double x, double mu, double sigma); + +double Clip(double a, double min, double max); + +std::vector genQuantificationRVec(double dataPoint, int n, std::vector rv); } // namespace FTAUtils -// Fault Tree class +/********************* Fault Tree Definition *********************/ + class FTAUtils::FaultTree { public: @@ -45,16 +70,35 @@ class FTAUtils::FaultTree } }; + /** + * Constructor for fault tree class + * Default value of root is the first node in file + */ FaultTree(std::string file_name, std::string root = ""); - FaultTree(std::set> & sets_link, std::map _node_base); + FaultTree(set_string & sets_link, std::map _node_base); + + /** + * Destructor + */ ~FaultTree(); std::string getRoot(); - // void printSets(); - // static void printSets(std::set> sets); - // static void printRow(std::set row); - std::set> getCutSets(); + + /** + * Function returns cut sets at the given point + * NOTE + * If MOCUS ran before this function call, min cut sets will be returned + */ + set_string getCutSets(); + + /** + * Builds m-ary fault tree + */ std::map buildTree(Parser parser); - std::set> computeMinimumCutSets(); + + /** + * Computes minimum cut sets based on MOCUS Algorithm + */ + set_string computeMinimumCutSets(); std::vector event(std::map); private: @@ -66,16 +110,35 @@ class FTAUtils::FaultTree std::string _root; // Cut sets container for in-place computations - std::set> _sets; + set_string _sets; - // Member functions + /** + * Translates string to opeartor + */ _operator_t str2Operator(std::string op); void rmSets(); _node * getNode(std::string name); + + /** + * Recursive call function to flood fill sets by expanding them based on + * operation ALGO + * 1. Iterate through entire list and match for own node's name + * 2. Replace self with children based on operation + * (i) . Replace self with children in same row if AND + * (ii). Replace self with child one per row if OR + * 3. Recurse on each non leaf child + * + * NOTE + * 1. Updates "sets" variable + * 2. Uses std::set 2d array, hence absorption and idempotence properties + * are implicit + */ void cutSetsExpand(_node * node); void removeSubsets(); }; +/************************** Parser Definition **************************/ + class FTAUtils::Parser { public: @@ -86,10 +149,33 @@ class FTAUtils::Parser FORMAT_UNDEF }; + /** + * Constructor for parser class + */ Parser(std::string fileName, parseFormatT format); + + /** + * Destructor for parser class + */ ~Parser(); + + /** + * Yields all records, populates the standard structure and returns + * This function acts as an abstract layer to hide different formats + * that might be supported in future + * + * Returns: Array of strings + */ + vector_string yieldLines(); + + /** + * Yields a single record, populates the standard structure and returns + * This function acts as an abstract layer to hide different formats + * that might be supported in future + * + * Returns: Array of strings + */ std::vector yieldLine(); - std::vector> yieldLines(); private: // Handle to file diff --git a/include/utils/QuantificationUtils.h b/include/utils/QuantificationUtils.h index cad5b8f0a4..ad72c5a7bd 100644 --- a/include/utils/QuantificationUtils.h +++ b/include/utils/QuantificationUtils.h @@ -5,16 +5,8 @@ #include "FaultTreeUtils.h" #include "MastodonUtils.h" -#define CLIP(A, MIN, MAX) (((A) < MIN) ? MIN : ((A) > MAX ? MAX : (A))) +/********************* Quantification Definition *********************/ -#define GEN_QUALIFICATION_R_VEC(gen, d, n, rv) \ - for (int index = 0; index < n; index++) \ - { \ - double dataPoint = d(gen); \ - rv.push_back(CLIP(dataPoint, 0.0, 1.0)); \ - } - -// Quantification class class FTAUtils::Quantification { public: @@ -39,8 +31,8 @@ class FTAUtils::Quantification RISK // = 1 }; - Quantification(std::map>> & params_double, - std::map>> & params_string, + Quantification(std::map & params_double, + std::map & params_string, std::map & params_int, std::map & params_bool, std::map & params_analysis_t, @@ -82,7 +74,18 @@ class FTAUtils::Quantification } int size = vec.size(); - ASSERT(size > 0, "Size 0 vector passed for stats computation"); + + // ASSERT + if (!(size > 0)) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => Root node " + "requested (%s) not found in heirarchy. " + "Size 0 vector passed for stats computation.\n", + __FILE__, + __LINE__); + abort(); + } _pe = vec[0]; @@ -114,33 +117,29 @@ class FTAUtils::Quantification p_max = (p_max > (vec.size() - 1)) ? vec.size() : p_max; return vec[p_min - 1] + ((p_i - p_min) * (vec[p_max - 1] - vec[p_min - 1])); } - /* - void printStats() - { - std::cout << "Mean: " << _mean << ", Median: " << _median << ", SD: " << _sd - << ", 5th: " << _p5 << ", 95th: " << _p95 << "\n\n"; - } - */ }; // TODO: 3 of the following are not supported as they need 1 arg rather than 2 // for their distribution. Remove these 3 if not needed else add support std::map _str2dist = { - /*{"EXP" , EXP },*/ {"GAMMA", GAMMA}, {"WEIBULL", WEIBULL}, {"EXTREME", EXTREME}, {"NORM", NORM}, {"LNORM", LNORM}, {"PE", PE}, - /*{"CHI_SQ" , CHI_SQ },*/ {"CAUCHY", CAUCHY}, {"FISHER_F", FISHER_F}, - /*{"STUDENT_T", STUDENT_T}*/ }; - std::vector> _cut_set_prob; + vector_double _cut_set_prob; std::map> _b_nodes; + + /** + * Generates probability vector for a specified distribution + * Nomenclature of a and b changes with distribution. + * eg., a => mean, b => std for NORM + */ std::vector getProbVector(_dist_t dist, double a, double b, @@ -149,42 +148,103 @@ class FTAUtils::Quantification std::vector im, _analysis_t analysis, bool uc); - std::vector> beProb(FTAUtils::Parser parser, - int n_sample, - int seed, - _analysis_t analysis, - std::vector intmes, - bool uncert); - - std::vector> computeCutSetProb(std::set> cut_sets, - int n, - bool bypass = false, - std::string bypass_key = "", - double bypass_value = 0); + + /** + * Parses and floods probabilties of basic elements + */ + vector_string beProb(FTAUtils::Parser parser, + int n_sample, + int seed, + _analysis_t analysis, + std::vector intmes, + bool uncert); + + /** + * Computes probability for a given cut set on a per set basis based on + * pre-flooded basic elem probs + */ + vector_double computeCutSetProb(std::set> cut_sets, + int n, + bool bypass = false, + std::string bypass_key = "", + double bypass_value = 0); + std::vector cutSetRowProb(std::set row, int n, bool sign_positive = true, bool bypass = false, std::string bypass_key = "", double bypass_value = 0); + + /** + * Computes accumulated probability for the entire cut set + * 3 Digests are computed: + * 0. Min Max + * 1. Upper Bound + * 2. Top Rare + * + * NOTE + * 1. Its better to compute these 2 together as they have common loops + * which can save time + * 2. Sets the vector cutSetProb + */ std::vector * cutSetProbWDigest(std::set> cut_sets, int n, bool only_min_max = false); + + /** + * Wrapper function for gate probability arith + * AND => a * b + * OR => (1 - a) * (1 - b) + */ double getGateProb(double a, double b, bool is_and); - std::vector> minCutIM(std::vector upperBound); - std::map>> - beIM(std::set> cut_sets, - int n, - std::vector upper_bound, - std::vector & count_v); - std::vector> linesToDouble(std::vector> lines); + + /** + * Computes cut-set Fussel-Vesely Important measures + */ + vector_double minCutIM(std::vector upperBound); + + /** + * For each basic element: + * 1. Calculate #occurrence in minimal cut set (beCount) + * 2. Store row index of all occurrence (beIndex) + */ + std::map beIM(std::set> cut_sets, + int n, + std::vector upper_bound, + std::vector & count_v); + + /** + * Translator function for string -> double + */ + vector_double linesToDouble(vector_string lines); + + /** + * Function for getting BIN mean values of intensity measure + */ std::vector getBinMeans(double im_lower, double im_upper, int n_bins); - std::vector hazInterp(std::vector> hazard, - std::vector im_bins); + + /** + * Function for interpolating the hazard curve based on range of intensity + * measures + */ + std::vector hazInterp(vector_double hazard, std::vector im_bins); + + /** + * Function for top event fragility + */ std::vector fragility(std::set> cut_sets, int n, std::vector im_bins, double & mu, double & sigma); + + /** + * Function for calculating risk by convoluting hazard and fragility + * Assign risk to basic event probabilites + * + * WARNING + * This consumes _b_nodes and then overwrites it + */ void computeRisk(int n, std::vector hazard); }; diff --git a/include/utils/Utils.h b/include/utils/Utils.h deleted file mode 100644 index e58a3dde50..0000000000 --- a/include/utils/Utils.h +++ /dev/null @@ -1,19 +0,0 @@ -#ifndef _UTILS_H -#define _UTILS_H - -#include -#include - -//-------------------------- MACRO UTILITIES BEGIN --------------------- -// A custom assert block for assertions -#define ASSERT(condition, statement, ...) \ - if (!(condition)) \ - { \ - fprintf(stderr, \ - "[ASSERT] In File: %s, Line: %d => " #statement "\n", \ - __FILE__, \ - __LINE__, \ - ##__VA_ARGS__); \ - abort(); \ - } -#endif // _UTILS_H diff --git a/src/utils/FaultTreeUtils.C b/src/utils/FaultTreeUtils.C index 990d23623d..2f79347af4 100644 --- a/src/utils/FaultTreeUtils.C +++ b/src/utils/FaultTreeUtils.C @@ -1,14 +1,9 @@ #include "FaultTreeUtils.h" #include "MastodonUtils.h" -// #include "DelimitedFileReader.h" -/* - * Constructor for fault tree class - * NOTE: Default value of root is the first node in file - */ -/*!public*/ +typedef FTAUtils::FaultTree::_node ft_node; + FTAUtils::FaultTree::FaultTree(std::string file_name, std::string root) -/*!endpublic*/ { FTAUtils::Parser parser = FTAUtils::Parser(file_name, FTAUtils::Parser::FORMAT_CSV); buildTree(parser); @@ -16,21 +11,27 @@ FTAUtils::FaultTree::FaultTree(std::string file_name, std::string root) // Override root node if not default if (root != "") { - ASSERT(getNode(root), - "Root node requested (%s) not found in heirarchy. Please the check " - "file %s", - root.c_str(), - file_name.c_str()); + // ASSERT + if (!(getNode(root))) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => Root node " + "requested (%s) not found in heirarchy. " + "Please the check file %s.\n", + __FILE__, + __LINE__, + root.c_str(), + file_name.c_str()); + abort(); + } + this->_root = root; } computeMinimumCutSets(); } -/*!public*/ -FTAUtils::FaultTree::FaultTree(std::set> & sets_link, - std::map _node_base) -/*!endpublic*/ +FTAUtils::FaultTree::FaultTree(set_string & sets_link, std::map _node_base) { std::string root = ""; _node_d_b = _node_base; @@ -60,20 +61,10 @@ FTAUtils::FaultTree::FaultTree(std::set> & sets_link, sets_link = _sets; } -/* - * Destructor - */ -/*!public*/ -FTAUtils::FaultTree::~FaultTree() -/*!endpublic*/ {} - -/* - * Builds m-ary fault tree - */ -/*!public*/ -std::map +FTAUtils::FaultTree::~FaultTree() {} + +std::map FTAUtils::FaultTree::buildTree(FTAUtils::Parser parser) -/*!endpublic*/ { std::vector line; while (true) @@ -102,27 +93,26 @@ FTAUtils::FaultTree::buildTree(FTAUtils::Parser parser) return _node_d_b; } -/* - * Translates string to opeartor - */ -/*!private*/ FTAUtils::FaultTree::_operator_t FTAUtils::FaultTree::str2Operator(std::string op) -/*!endprivate*/ { std::string op_s = FTAUtils::str2Upper(op, true); - ASSERT(_opDict.count(op_s) != 0, "Illegal Operator found: %s", op.c_str()); - + // ASSERT + if (_opDict.count(op_s) == 0) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Illegal Operator found: %s.\n", + __FILE__, + __LINE__, + op.c_str()); + abort(); + } return _opDict[op_s]; } -/* - * Computes minimum cut sets based on MOCUS Algorithm - */ -/*!public*/ -std::set> +set_string FTAUtils::FaultTree::computeMinimumCutSets() -/*!endpublic*/ { // Clearing up sets vector for safety rmSets(); @@ -149,21 +139,16 @@ FTAUtils::FaultTree::computeMinimumCutSets() return _sets; } -/*!private*/ void FTAUtils::FaultTree::removeSubsets() -/*!endprivate*/ { std::set rm_its; uint64_t index1 = 0; - for (std::set>::iterator row1_it = _sets.begin(); row1_it != _sets.end(); - ++row1_it) + for (set_string::iterator row1_it = _sets.begin(); row1_it != _sets.end(); ++row1_it) { uint64_t index2 = index1 + 1; - for (std::set>::iterator row2_it = next(row1_it, 1); - row2_it != _sets.end(); - ++row2_it) + for (set_string::iterator row2_it = next(row1_it, 1); row2_it != _sets.end(); ++row2_it) { bool row1_is_subset = includes(row2_it->begin(), row2_it->end(), row1_it->begin(), row1_it->end()); @@ -189,34 +174,38 @@ FTAUtils::FaultTree::removeSubsets() } } -/* - * Recursive call function to flood fill sets by expanding them based on - * operation ALGO: - * 1. Iterate through entire list and match for own node's name - * 2. Replace self with children based on operation - * (i) . Replace self with children in same row if AND - * (ii). Replace self with child one per row if OR - * 3. Recurse on each non leaf child - * - * NOTE: 1. Updates "sets" variable - * 2. Uses std::set 2d array, hence absorption and idempotence properties - * are implicit - */ -/*!private*/ void FTAUtils::FaultTree::cutSetsExpand(_node * node) -/*!endprivate*/ { - ASSERT(node, "NULL node encountered"); - ASSERT(node->_child.size() != 0, "Empty child list found for node '%s'", node->_name.c_str()); + // ASSERT + if (!node) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "NULL node encountered.\n", + __FILE__, + __LINE__); + abort(); + } + + // ASSERT + if (node->_child.size() == 0) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Empty child list found for node '%s'.\n", + __FILE__, + __LINE__, + node->_name.c_str()); + abort(); + } // New rows which have to be appended at end. Adding them might // result in data hazards - std::set> mk_rows; + set_string mk_rows; // 1. Iterate through entire list and match for own node's name - for (std::set>::iterator row_it = _sets.begin(); row_it != _sets.end(); - ++row_it) + for (set_string::iterator row_it = _sets.begin(); row_it != _sets.end(); ++row_it) { std::set row(*row_it); // Search for the element in row @@ -253,8 +242,16 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) break; default: - ASSERT(false, "Unknown Operator found!"); - break; + { + // ASSERT + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Unknown Operator found!.\n", + __FILE__, + __LINE__); + abort(); + } + break; } } @@ -267,7 +264,7 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) } // Add newly created rows - for (std::set>::iterator it = mk_rows.begin(); it != mk_rows.end(); ++it) + for (set_string::iterator it = mk_rows.begin(); it != mk_rows.end(); ++it) { std::set row(*it); _sets.insert(*it); @@ -282,109 +279,59 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) } } -/*!private*/ void FTAUtils::FaultTree::rmSets() -/*!endprivate*/ { - for (std::set>::iterator row_it = _sets.begin(); row_it != _sets.end(); - ++row_it) + for (set_string::iterator row_it = _sets.begin(); row_it != _sets.end(); ++row_it) { _sets.erase(row_it); } } -/*!private*/ -FTAUtils::FaultTree::_node * +ft_node * FTAUtils::FaultTree::getNode(std::string name) -/*!endprivate*/ { return (_node_d_b.count(name) != 0) ? _node_d_b[name] : NULL; } +{ + return (_node_d_b.count(name) != 0) ? _node_d_b[name] : NULL; +} -/*!public*/ std::string FTAUtils::FaultTree::getRoot() -/*!endpublic*/ { return _root; } - -/*!public*/ -// void FTAUtils::FaultTree::printSets(std::set> sets) -/*!endpublic*/ -/* { - std::cout << "------------- SETS BEGIN ----------------- " << std::endl; - for (std::set>::iterator row = sets.begin(); row != sets.end(); ++row) - { - printRow(*row); - std::cout << std::endl; - } - std::cout << "-------------- SETS END ------------------ " << std::endl; + return _root; } -*/ - -/*!public*/ -// void -// FTAUtils::FaultTree::printSets() -/*!endpublic*/ -// { printSets(_sets); } - -/*!public*/ -// void -// FTAUtils::FaultTree::printRow(std::set row) -/*!endpublic*/ -/* + +set_string +FTAUtils::FaultTree::getCutSets() { - for (std::set::iterator col = row.begin(); col != row.end(); ++col) - { - std::cout << *col << ", "; - } + return _sets; } -*/ - -/* - * Function returns cut sets at the given point - * NOTE: If MOCUS ran before this function call, min cut sets will be returned - */ -/*!public*/ -std::set> -FTAUtils::FaultTree::getCutSets() -/*!endpublic*/ { return _sets; } -// **************************Parser Definition************************** -/* - * Constructor for parser class - */ -/*!public*/ FTAUtils::Parser::Parser(std::string fileName, FTAUtils::Parser::parseFormatT format) -/*!endpublic*/ { // Assertion on supported parsing formats - ASSERT(format == FORMAT_CSV, "Unsupported parse format"); + if (format != FORMAT_CSV) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Unsupported parse format.\n", + __FILE__, + __LINE__); + abort(); + } + fileP = new std::ifstream; fileP->open(fileName, std::ifstream::in); - // ASSERT( fileP->is_open(), "Unable to open file: %s", fileName.c_str()); if (!fileP->is_open()) throw FTAUtils::CException("Unable to open file."); } -/* - * Destructor for parser class - */ -/*!public*/ -FTAUtils::Parser::~Parser() -/*!endpublic*/ {} - -/* - * Yields all records, populates the standard structure and returns - * This function acts as an abstract layer to hide different formats - * that might be supported in future - * - * Returns: Array of strings - */ -/*!public*/ -std::vector> +FTAUtils::Parser::~Parser() {} + +vector_string FTAUtils::Parser::yieldLines() -/*!endpublic*/ { - std::vector> lines; + vector_string lines; std::vector line; while (true) { @@ -398,19 +345,19 @@ FTAUtils::Parser::yieldLines() return lines; } -/* - * Yields a single record, populates the standard structure and returns - * This function acts as an abstract layer to hide different formats - * that might be supported in future - * - * Returns: Array of strings - */ -/*!public*/ std::vector FTAUtils::Parser::yieldLine() -/*!endpublic*/ { - ASSERT(fileP->is_open(), "Illegal call to yieldLine. File not open yet!"); + // ASSERT + if (!(fileP->is_open())) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Unsupported parse format.\n", + __FILE__, + __LINE__); + abort(); + } std::string buffer; std::vector line; @@ -428,9 +375,6 @@ FTAUtils::Parser::yieldLine() return line; } -/* - * String trim method which removes all leading and lagging whitespace in a string - */ std::string FTAUtils::trim(const std::string & str) { @@ -445,11 +389,6 @@ FTAUtils::trim(const std::string & str) return str.substr(first, (last - first + 1)); } -// defined in Utils.cpp in FTA code - -/* - * Converts ASCII string to upper case - */ std::string FTAUtils::str2Upper(const std::string & str_in, bool trim_input) { @@ -462,13 +401,8 @@ FTAUtils::str2Upper(const std::string & str_in, bool trim_input) return str; } -/* - * Returns interpolated value at x from parallel arrays ( xData, yData ) - * Assumes that xData has at least two elements, is sorted and is strictly monotonic increasing - * boolean argument extrapolate determines behaviour beyond ends of array (if needed) - */ double -FTAUtils::interpolate(std::vector> data, double x, bool extrapolate) +FTAUtils::interpolate(vector_double data, double x, bool extrapolate) { int size = data.size(); @@ -505,7 +439,6 @@ FTAUtils::interpolate(std::vector> data, double x, bool extr return yL + dydx * (x - xL); } -// Phi(-∞, x) aka N(x) double FTAUtils::normalCDF(double x) { diff --git a/src/utils/QuantificationUtils.C b/src/utils/QuantificationUtils.C index 5e9e29603a..e740a5d4fe 100644 --- a/src/utils/QuantificationUtils.C +++ b/src/utils/QuantificationUtils.C @@ -1,45 +1,30 @@ #include "QuantificationUtils.h" -// #include "DelimitedFileReader.h" - -// A macro to reduce typing while printing stats -#define QUANT_BE_STAT_PRINT_HELPER(index, i) \ - "PE: " << be_stats[index][i][0] << ", Mean: " << be_stats[index][i][1] \ - << ", Median: " << be_stats[index][i][2] << ", SD: " << be_stats[index][i][3] \ - << ", 5th: " << be_stats[index][i][4] << ", 95th: " << be_stats[index][i][5] << "\n" - -/* - * Constructor for qualifications class - */ -/*!public*/ -FTAUtils::Quantification::Quantification( - std::map>> & params_double, - std::map>> & params_string, - std::map & params_int, - std::map & params_bool, - std::map & params_analysis_t, - std::string events_file, - std::string events_prob_file, - _analysis_t analysis, - std::string hazard_file, - double im_lower, - double im_upper, - int n_bins, - bool uncertainty, - std::string root, - int n_sample, - int seed) -/*!endpublic*/ + +FTAUtils::Quantification::Quantification(std::map & params_double, + std::map & params_string, + std::map & params_int, + std::map & params_bool, + std::map & params_analysis_t, + std::string events_file, + std::string events_prob_file, + _analysis_t analysis, + std::string hazard_file, + double im_lower, + double im_upper, + int n_bins, + bool uncertainty, + std::string root, + int n_sample, + int seed) { //---------------- SAVE PARAMETERS ------------------------------------------ std::vector params_double1; params_double1.push_back(im_lower); params_double1.push_back(im_upper); - std::vector> params_double2; + vector_double params_double2; params_double2.push_back(params_double1); - params_double.insert( - std::pair>>("IM", params_double2)); - + params_double.insert(pair_string_double("IM", params_double2)); params_int.insert(std::pair("n_bins", n_bins)); params_int.insert(std::pair("nsamp", n_sample)); params_int.insert(std::pair("seed", seed)); @@ -50,25 +35,22 @@ FTAUtils::Quantification::Quantification( // Construct the fault tree with MOCUS algo for minimal cut sets FTAUtils::FaultTree ft = FTAUtils::FaultTree(events_file, root); Parser parser_events = Parser(events_file, Parser::FORMAT_CSV); - std::vector> lines_events = parser_events.yieldLines(); - std::set> cut_sets = ft.getCutSets(); - // ft.printSets(); - params_string.insert( - std::pair>>("events_files", lines_events)); + vector_string lines_events = parser_events.yieldLines(); + set_string cut_sets = ft.getCutSets(); + params_string.insert(std::pair("events_files", lines_events)); // Read basic events file FTAUtils::Parser parser_event_prob = FTAUtils::Parser(events_prob_file, FTAUtils::Parser::FORMAT_CSV); - std::vector> lines_events_prob; + vector_string lines_events_prob; if (analysis == FRAGILITY) { // Read and interpolate hazard curve FTAUtils::Parser parser_hazard = FTAUtils::Parser(hazard_file, FTAUtils::Parser::FORMAT_CSV); - std::vector> lines_hazard = parser_hazard.yieldLines(); - std::vector> hazard = FTAUtils::Quantification::linesToDouble(lines_hazard); - params_double.insert( - std::pair>>("hazard", hazard)); + vector_string lines_hazard = parser_hazard.yieldLines(); + vector_double hazard = FTAUtils::Quantification::linesToDouble(lines_hazard); + params_double.insert(pair_string_double("hazard", hazard)); // List of Intensity Measure bin values std::vector im_bins = getBinMeans(im_lower, im_upper, n_bins); @@ -82,12 +64,6 @@ FTAUtils::Quantification::Quantification( // Top event fragility (lognormal parameters) double mu, sigma; fragility(cut_sets, n_bins, im_bins, mu, sigma); - /* - std::cout << "----------- LN PARAMS BEGIN --------------" << std::endl; - std::cout << "mu: " << mu << std::endl; - std::cout << "sigma: " << sigma << std::endl; - std::cout << "------------ LN PARAMS END ---------------" << std::endl; - */ // Dictionary of basic events risk (convoluting fragility and hazard) computeRisk(n_bins, hazard_freq); @@ -97,105 +73,59 @@ FTAUtils::Quantification::Quantification( // Dictionary of basic events risk (risk inputs) lines_events_prob = beProb(parser_event_prob, n_sample, seed, analysis, std::vector(), uncertainty); - /* - std::cout << "------------ LN PARAMS END ---------------" << std::endl; - */ } - params_string.insert(std::pair>>( - "basic_events", lines_events_prob)); + + params_string.insert(std::pair("basic_events", lines_events_prob)); // Calculate minimal cut sets probability // Digest cut set probability to unified probability // Adding 1 to accomodate mean as 0th element std::vector * digest = cutSetProbWDigest(cut_sets, n_sample + 1); - std::vector> mc_im = minCutIM(digest[1]); + vector_double mc_im = minCutIM(digest[1]); Stats s0(digest[0]); Stats s1(digest[1]); Stats s2(digest[2]); std::vector count_v; - std::map>> be_stats = - beIM(cut_sets, n_sample + 1, digest[1], count_v); - - //---------------- PRINTING -------------------------------------------------- - /* - std::cout << "---------- PROBABILITY BEGIN -------------" << std::endl; - std::cout << "1. Exact solution: " << s0._pe << std::endl; - s0.printStats(); - std::cout << "2. Upper Bound : " << s1._pe << std::endl; - s1.printStats(); - std::cout << "3. Rare Event : " << s2._pe << std::endl; - s2.printStats(); - */ + std::map be_stats = beIM(cut_sets, n_sample + 1, digest[1], count_v); // return the FTA top event risk std::vector fta_0; - std::vector> fta_1; + vector_double fta_1; fta_0.push_back(s0._pe); fta_0.push_back(s1._pe); fta_0.push_back(s2._pe); fta_1.push_back(fta_0); params_double.insert(std::make_pair("fta", fta_1)); - /* - std::cout << "----------- PROBABILITY END --------------" << std::endl; - std::cout << "-------- CUT SET DETAILS BEGIN -----------" << std::endl; - - int i = 0; - for (std::set>::iterator row = cut_sets.begin(); row != cut_sets.end(); - ++row) { - ft.printRow(*row); - std::cout << "\t(PE: " << _cut_set_prob[i][0] << ", IM: " << mc_im[i][0] << ")" - << std::endl; - i++; - } - - std::cout << "--------- CUT SET DETAILS END ------------" << std::endl; - i = 0; - for (std::map>::iterator bn_it = _b_nodes.begin(); - bn_it != _b_nodes.end(); ++bn_it) { - std::cout << bn_it->first << " (" << count_v[i] << ")" << std::endl; - std::cout << "\t[FV] " QUANT_BE_STAT_PRINT_HELPER("fv", i); - std::cout << "\t[RRR] " QUANT_BE_STAT_PRINT_HELPER("rrr", i); - std::cout << "\t[RIR] " QUANT_BE_STAT_PRINT_HELPER("rir", i); - std::cout << "\t[RRI] " QUANT_BE_STAT_PRINT_HELPER("rri", i); - std::cout << "\t[RII] " QUANT_BE_STAT_PRINT_HELPER("rii", i); - std::cout << "\t[BI] " QUANT_BE_STAT_PRINT_HELPER("bi", i); - i++; - } - */ - // Return the Risk Reduction Ratio for B1, B2, B3, B4, B5 - params_double.insert( - std::pair>>("fv", be_stats["fv"])); - params_double.insert( - std::pair>>("rrr", be_stats["rrr"])); - params_double.insert( - std::pair>>("rir", be_stats["rir"])); + params_double.insert(pair_string_double("fv", be_stats["fv"])); + params_double.insert(pair_string_double("rrr", be_stats["rrr"])); + params_double.insert(pair_string_double("rir", be_stats["rir"])); // Return the Risk Reduction Difference for B1, B2, B3, B4, B5 - params_double.insert( - std::pair>>("rri", be_stats["rri"])); - params_double.insert( - std::pair>>("rii", be_stats["rii"])); - params_double.insert( - std::pair>>("bi", be_stats["bi"])); + params_double.insert(pair_string_double("rri", be_stats["rri"])); + params_double.insert(pair_string_double("rii", be_stats["rii"])); + params_double.insert(pair_string_double("bi", be_stats["bi"])); } -/* - * Function for interpolating the hazard curve based on range of intensity - * measures - */ -/*!private*/ std::vector -FTAUtils::Quantification::hazInterp(std::vector> hazard, - std::vector im_bins) -/*!endprivate*/ +FTAUtils::Quantification::hazInterp(vector_double hazard, std::vector im_bins) { std::vector data; // Convert all hazard values to log10 for (int index = 0; index < hazard.size(); index++) { - ASSERT(hazard[index].size() == 2, "hazard dimension error"); + // ASSERT + if (hazard[index].size() != 2) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "hazard dimension error.\n", + __FILE__, + __LINE__); + abort(); + } + hazard[index][0] = log10(hazard[index][0]); hazard[index][1] = log10(hazard[index][1]); } @@ -208,13 +138,8 @@ FTAUtils::Quantification::hazInterp(std::vector> hazard, return data; } -/* - * Function for getting BIN mean values of intensity measure - */ -/*!private*/ std::vector FTAUtils::Quantification::getBinMeans(double im_lower, double im_upper, int n_bins) -/*!endprivate*/ { std::vector bins; double delta = (im_upper - im_lower) / (n_bins); @@ -226,15 +151,10 @@ FTAUtils::Quantification::getBinMeans(double im_lower, double im_upper, int n_bi return bins; } -/* - * Translator function for string -> double - */ -/*!private*/ -std::vector> -FTAUtils::Quantification::linesToDouble(std::vector> lines) -/*!endprivate*/ +vector_double +FTAUtils::Quantification::linesToDouble(vector_string lines) { - std::vector> lines_double; + vector_double lines_double; for (int index = 0; index < lines.size(); index++) { std::vector double_vector(lines[index].size()); @@ -247,12 +167,23 @@ FTAUtils::Quantification::linesToDouble(std::vector> li return lines_double; } -/* - * Generates probability vector for a specified distribution - * Nomenclature of a and b changes with distribution. - * eg., a => mean, b => std for NORM - */ -/*!private*/ +double +FTAUtils::Clip(double a, double min, double max) +{ + return ((a) < min) ? min : ((a) > max ? max : (a)); +} + +std::vector +FTAUtils::genQuantificationRVec(double dataPoint, int n, std::vector rv) +{ + for (int index = 0; index < n; index++) + { + // double dataPoint = double(gen); + rv.push_back(FTAUtils::Clip(dataPoint, 0.0, 1.0)); + } + return rv; +} + std::vector FTAUtils::Quantification::getProbVector(_dist_t dist, double a, @@ -262,7 +193,6 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, std::vector im, _analysis_t analysis, bool uc) -/*!endprivate*/ { std::vector rv; @@ -270,7 +200,17 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, if (analysis == FRAGILITY) { - ASSERT(dist == LNORM, "unsupported distribution for fragility analysis"); + // ASSERT + if (dist != LNORM) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "unsupported distribution for fragility analysis.\n", + __FILE__, + __LINE__); + abort(); + } + for (int index = 0; index < im.size(); index++) { rv.push_back(FTAUtils::normalCDF(log(im[index] / a) / b)); @@ -284,16 +224,24 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, { case PE: { - GEN_QUALIFICATION_R_VEC(a, double, 2, rv); + rv = FTAUtils::genQuantificationRVec(double(a), 2, rv); } break; case NORM: { - GEN_QUALIFICATION_R_VEC(a, double, 2, rv); + rv = FTAUtils::genQuantificationRVec(double(a), 2, rv); } break; default: - ASSERT(false, "Un-supported dist found"); + { + // ASSERT + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Un-supported dist found.\n", + __FILE__, + __LINE__); + abort(); + } } } else @@ -302,39 +250,45 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, { case PE: { - GEN_QUALIFICATION_R_VEC(a, double, n, rv); + rv = FTAUtils::genQuantificationRVec(double(a), n, rv); } break; case NORM: { - rv.push_back(CLIP(a, 0, 1)); + double rv_norm = FTAUtils::Clip(a, 0, 1); + rv.push_back(rv_norm); + std::normal_distribution d(a, b); - GEN_QUALIFICATION_R_VEC(gen, d, n, rv); + + rv = FTAUtils::genQuantificationRVec(d(gen), n, rv); } break; default: - ASSERT(false, "Un-supported dist found"); + { + // ASSERT + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Un-supported dist found.\n", + __FILE__, + __LINE__); + abort(); + } } } } return rv; } -/* - * Parses and floods probabilties of basic elements - */ -/*!private*/ -std::vector> +vector_string FTAUtils::Quantification::beProb(FTAUtils::Parser parser, int n_sample, int seed, _analysis_t analysis, std::vector intmes, bool uncert) -/*!endprivate*/ { std::vector line; - std::vector> lines; + vector_string lines; while (true) { line = parser.yieldLine(); @@ -354,27 +308,13 @@ FTAUtils::Quantification::beProb(FTAUtils::Parser parser, return lines; } -/* - * Computes accumulated probability for the entire cut set - * 3 Digests are computed: - * 0. Min Max - * 1. Upper Bound - * 2. Top Rare - * NOTE: 1. Its better to compute these 2 together as they have common loops - * which can save time - * 2. Sets the vector cutSetProb - */ -/*!private*/ std::vector * -FTAUtils::Quantification::cutSetProbWDigest(std::set> cut_sets, - int n, - bool only_min_max) -/*!endprivate*/ +FTAUtils::Quantification::cutSetProbWDigest(set_string cut_sets, int n, bool only_min_max) { std::vector * digest = new std::vector[3]; // Digest 2: Generate power set to generate all possible combinations - std::vector> ps_p; + vector_double ps_p; int pow_set_size = pow(2, cut_sets.size()); // Start from 1 as NULL set is not needed for (int index = 1; index < pow_set_size; index++) @@ -410,9 +350,7 @@ FTAUtils::Quantification::cutSetProbWDigest(std::set> cut_ double top_rare = 0; if (!only_min_max) { - for (std::vector>::iterator it = _cut_set_prob.begin(); - it != _cut_set_prob.end(); - ++it) + for (vector_double::iterator it = _cut_set_prob.begin(); it != _cut_set_prob.end(); ++it) { top_rare += (*it)[index]; upper_bound *= 1 - (*it)[index]; @@ -420,7 +358,7 @@ FTAUtils::Quantification::cutSetProbWDigest(std::set> cut_ digest[1].push_back(1 - upper_bound); digest[2].push_back(top_rare); } - for (std::vector>::iterator it = ps_p.begin(); it != ps_p.end(); ++it) + for (vector_double::iterator it = ps_p.begin(); it != ps_p.end(); ++it) { min_max += (*it)[index]; } @@ -430,32 +368,21 @@ FTAUtils::Quantification::cutSetProbWDigest(std::set> cut_ return digest; } -/* - * Computes probability for a given cut set on a per set basis based on - * pre-flooded basic elem probs - */ -/*!private*/ -std::vector> -FTAUtils::Quantification::computeCutSetProb(std::set> cut_sets, - int n, - bool bypass, - std::string bypass_key, - double bypass_value) -/*!endprivate*/ +vector_double +FTAUtils::Quantification::computeCutSetProb( + set_string cut_sets, int n, bool bypass, std::string bypass_key, double bypass_value) { // For each row in cut set, compute: // 1. For each sample in the generated vector, compute AND gate // probability analysis for each basic element - std::vector> quant; - for (std::set>::iterator row = cut_sets.begin(); row != cut_sets.end(); - ++row) + vector_double quant; + for (set_string::iterator row = cut_sets.begin(); row != cut_sets.end(); ++row) { quant.push_back(cutSetRowProb(*row, n, true, bypass, bypass_key, bypass_value)); } return quant; } -/*!private*/ std::vector FTAUtils::Quantification::cutSetRowProb(std::set row, int n, @@ -463,7 +390,6 @@ FTAUtils::Quantification::cutSetRowProb(std::set row, bool bypass, std::string bypass_key, double bypass_value) -/*!endprivate*/ { std::vector prob; double sign_m = sign_positive ? 1 : -1; @@ -472,8 +398,17 @@ FTAUtils::Quantification::cutSetRowProb(std::set row, double value = 1; for (std::set::iterator col = row.begin(); col != row.end(); ++col) { - ASSERT( - _b_nodes.find(*col) != _b_nodes.end(), "'%s' key not found in _b_nodes", (*col).c_str()); + // ASSERT + if (_b_nodes.find(*col) == _b_nodes.end()) + { + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "'%s' key not found in _b_nodes.\n", + __FILE__, + __LINE__, + (*col).c_str()); + abort(); + } double bv = (bypass && (bypass_key == *col)) ? bypass_value : _b_nodes[*col][index]; value = getGateProb(value, bv, true); @@ -483,25 +418,16 @@ FTAUtils::Quantification::cutSetRowProb(std::set row, return prob; } -/* - * Wrapper function for gate probability arith - * AND => a * b - * OR => (1 - a) * (1 - b) - */ -/*!private*/ double FTAUtils::Quantification::getGateProb(double a, double b, bool is_and) -/*!endprivate*/ { return is_and ? a * b : (1.0 - a) * (1.0 - b); } +{ + return is_and ? a * b : (1.0 - a) * (1.0 - b); +} -/* - * Computes cut-set Fussel-Vesely Important measures - */ -/*!private*/ -std::vector> +vector_double FTAUtils::Quantification::minCutIM(std::vector upper_bound) -/*!endprivate*/ { - std::vector> mc_i_m; + vector_double mc_i_m; for (int row = 0; row < _cut_set_prob.size(); row++) { std::vector mc_i_m_row; @@ -514,15 +440,8 @@ FTAUtils::Quantification::minCutIM(std::vector upper_bound) return mc_i_m; } -/* - * Function for calculating risk by convoluting hazard and fragility - * Assign risk to basic event probabilites - * WARNING: This consumes _b_nodes and then overwrites it - */ -/*!private*/ void FTAUtils::Quantification::computeRisk(int n, std::vector hazard) -/*!endprivate*/ { for (std::map>::iterator bn_it = _b_nodes.begin(); bn_it != _b_nodes.end(); @@ -537,17 +456,9 @@ FTAUtils::Quantification::computeRisk(int n, std::vector hazard) } } -/* - * Function for top event fragility - */ -/*!private*/ std::vector -FTAUtils::Quantification::fragility(std::set> cut_sets, - int n, - std::vector im_bins, - double & mu, - double & sigma) -/*!endprivate*/ +FTAUtils::Quantification::fragility( + set_string cut_sets, int n, std::vector im_bins, double & mu, double & sigma) { // 1. Calculate TOP event fragility using min-max approach (exact) {digest[0]} std::vector * digest = cutSetProbWDigest(cut_sets, n, true); @@ -568,19 +479,13 @@ FTAUtils::Quantification::fragility(std::set> cut_sets, return digest[0]; } -/* For each basic element: - * 1. Calculate #occurrence in minimal cut set (beCount) - * 2. Store row index of all occurrence (beIndex) - */ -/*!private*/ -std::map>> -FTAUtils::Quantification::beIM(std::set> cut_sets, +std::map +FTAUtils::Quantification::beIM(set_string cut_sets, int n, std::vector upper_bound, std::vector & count_v) -/*!endprivate*/ { - std::map>> stats; + std::map stats; for (std::map>::iterator bn_it = _b_nodes.begin(); bn_it != _b_nodes.end(); ++bn_it) @@ -588,9 +493,7 @@ FTAUtils::Quantification::beIM(std::set> cut_sets, // Generate available vector to save computation on per index loop std::vector available; int count = 0; - for (std::set>::iterator cs_it = cut_sets.begin(); - cs_it != cut_sets.end(); - ++cs_it) + for (set_string::iterator cs_it = cut_sets.begin(); cs_it != cut_sets.end(); ++cs_it) { bool is_a = cs_it->find(bn_it->first) != cs_it->end(); count += is_a; @@ -601,10 +504,8 @@ FTAUtils::Quantification::beIM(std::set> cut_sets, { std::vector fv, rrr, rir, rri, rii, bi; - std::vector> mc_p1 = - computeCutSetProb(cut_sets, n, true, bn_it->first, 1); - std::vector> mc_p0 = - computeCutSetProb(cut_sets, n, true, bn_it->first, 0); + vector_double mc_p1 = computeCutSetProb(cut_sets, n, true, bn_it->first, 1); + vector_double mc_p0 = computeCutSetProb(cut_sets, n, true, bn_it->first, 0); // OR it with loop merging for (int index = 0; index < n; index++) diff --git a/unit/src/TestFaultTreeUtils.C b/unit/src/TestFaultTreeUtils.C index 2c909cc95e..ccfc21d31a 100644 --- a/unit/src/TestFaultTreeUtils.C +++ b/unit/src/TestFaultTreeUtils.C @@ -12,9 +12,11 @@ // QUANTIFICATIONUTILS includes #include "FaultTreeUtils.h" +typedef FTAUtils::FaultTree::_node ft_node_test; + // Helper to make sure the test tree is correct. void -assertTree(std::map tree_node) +assertTree(std::map tree_node) { // TOP EXPECT_EQ(tree_node["TOP"]->_name, "TOP"); @@ -65,17 +67,17 @@ assertTree(std::map tree_node) // Helper to make sure the mocus is correct. void -assertMocus(std::set> tree_mocus) +assertMocus(set_string tree_mocus) { - std::set>::iterator it_tree_mocus; + set_string::iterator it_tree_mocus; std::set tree_mocus_children, set_mocus, set_gold; std::set::iterator it_tree_mocus_children; - std::set> matrix_gold = { + set_string matrix_gold = { {"B2", "B1"}, {"B4", "B1"}, {"B1", "B5", "B3"}, {"B5", "B2", "B3"}, {"B4", "B5", "B3"}}; - std::set>::iterator it_matrix_gold = matrix_gold.begin(); + set_string::iterator it_matrix_gold = matrix_gold.begin(); EXPECT_EQ(tree_mocus.size(), matrix_gold.size()); @@ -91,11 +93,11 @@ assertMocus(std::set> tree_mocus) } // Create a test event list. -std::map +std::map eventList() { - std::map _node_d_b; + std::map _node_d_b; std::vector> line = {{"TOP", "AND", "GATE1", "GATE2"}, {"GATE1", "OR", "FT-N/m-1", "FT-N/m-2", "FT-N/m-3"}, @@ -115,7 +117,7 @@ eventList() for (int i = 0; i < line.size(); i++) { - FTAUtils::FaultTree::_node * node = new FTAUtils::FaultTree::_node(line[i][0], op[i]); + ft_node_test * node = new ft_node_test(line[i][0], op[i]); // Add children for (int j = 2; j < line[i].size(); j++) node->_child.push_back(line[i][j]); @@ -127,7 +129,7 @@ eventList() } // For input -std::vector file_lists_fault_tree_1, file_lists_fault_tree_2, file_lists_fault_tree_3; +std::vector file_lists; TEST(FTAUtils, FaultTree) { @@ -139,9 +141,9 @@ TEST(FTAUtils, FaultTree) try { - file_lists_fault_tree_1 = {"not_a_valid_filename.txt", ""}; + file_lists = {"not_a_valid_filename.txt", ""}; - FTAUtils::FaultTree(file_lists_fault_tree_1[0], file_lists_fault_tree_1[1]); + FTAUtils::FaultTree(file_lists[0], file_lists[1]); } catch (FTAUtils::CException e) { @@ -154,23 +156,21 @@ TEST(FTAUtils, FaultTree) // ---------------- Test creating tree from list ---------------- - std::map tree_node_from_list = eventList(); - std::set> tree_mocus_from_list; + std::map tree_node_from_list = eventList(); + set_string tree_mocus_from_list; FTAUtils::FaultTree ft_from_list = FTAUtils::FaultTree(tree_mocus_from_list, tree_node_from_list); assertTree(tree_node_from_list); assertMocus(tree_mocus_from_list); // ---------------- Test creating tree from file ---------------- - file_lists_fault_tree_2 = {"logic2.txt", ""}; + file_lists = {"logic2.txt", ""}; - FTAUtils::FaultTree ft_from_file_2 = - FTAUtils::FaultTree(file_lists_fault_tree_2[0], file_lists_fault_tree_2[1]); - FTAUtils::Parser parser_events_2 = - FTAUtils::Parser(file_lists_fault_tree_2[0], FTAUtils::Parser::FORMAT_CSV); - std::map tree_node_from_file = + FTAUtils::FaultTree ft_from_file_2 = FTAUtils::FaultTree(file_lists[0], file_lists[1]); + FTAUtils::Parser parser_events_2 = FTAUtils::Parser(file_lists[0], FTAUtils::Parser::FORMAT_CSV); + std::map tree_node_from_file = ft_from_file_2.buildTree(parser_events_2); - std::set> tree_mocus_from_file = ft_from_file_2.computeMinimumCutSets(); + set_string tree_mocus_from_file = ft_from_file_2.computeMinimumCutSets(); assertTree(tree_node_from_file); assertMocus(tree_mocus_from_file); } diff --git a/unit/src/TestQuantificationUtils.C b/unit/src/TestQuantificationUtils.C index 37037e9c0a..2295e91c84 100644 --- a/unit/src/TestQuantificationUtils.C +++ b/unit/src/TestQuantificationUtils.C @@ -13,7 +13,7 @@ #include "QuantificationUtils.h" // For input -std::vector file_lists_quantification_utils; +std::vector file_lists_q; std::vector IM; int nbins; bool uncertainty; @@ -53,7 +53,7 @@ TEST(FTAUtils, Quantification) try { - file_lists_quantification_utils = { + file_lists_q = { "not_a_valid_filename.txt", "not_a_valid_filename.txt", "not_a_valid_filename.txt"}; IM = {0.1, 4}; nbins = 15; @@ -63,10 +63,10 @@ TEST(FTAUtils, Quantification) params_int_1, params_bool_1, params_analysis_t_1, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins); @@ -82,7 +82,7 @@ TEST(FTAUtils, Quantification) try { - file_lists_quantification_utils = {"", "", ""}; + file_lists_q = {"", "", ""}; IM = {0.1, 4}; nbins = 15; @@ -91,10 +91,10 @@ TEST(FTAUtils, Quantification) params_int_2, params_bool_2, params_analysis_t_2, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins); @@ -110,7 +110,7 @@ TEST(FTAUtils, Quantification) // ############## Inputs for Fragility ############## - file_lists_quantification_utils = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; + file_lists_q = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; IM = {0.1, 4}; nbins = 15; FTAUtils::Quantification(params_double_3, @@ -118,10 +118,10 @@ TEST(FTAUtils, Quantification) params_int_3, params_bool_3, params_analysis_t_3, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::FRAGILITY, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins); @@ -172,7 +172,7 @@ TEST(FTAUtils, Quantification) try { - file_lists_quantification_utils = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; + file_lists_q = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; IM = {0.1, 4}; nbins = -15; @@ -184,10 +184,10 @@ TEST(FTAUtils, Quantification) params_int_4, params_bool_4, params_analysis_t_4, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::FRAGILITY, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins); @@ -203,14 +203,15 @@ TEST(FTAUtils, Quantification) // ############## Inputs for Risk analysis (not fragility) ############## - file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt"}; + file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt"}; + FTAUtils::Quantification(params_double_5, params_string_5, params_int_5, params_bool_5, params_analysis_t_5, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK); // >>>>>>>> logic Value Check <<<<<<<< @@ -233,7 +234,9 @@ TEST(FTAUtils, Quantification) {"B3", "PE", "0.03"}, {"B4", "PE", "0.04"}, {"B5", "PE", "0.05"}}; + std::vector> basic_events_5 = params_string_5["basic_events"]; + EXPECT_EQ(basic_events_5, basic_events_matrix5); // >>>>>>>> antype Value Check <<<<<<<< @@ -254,7 +257,7 @@ TEST(FTAUtils, Quantification) // ###### Testing for input errors making sure parameters are input correctly ###### - file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; nbins = 15; uncertainty = true; @@ -266,10 +269,10 @@ TEST(FTAUtils, Quantification) params_int_6, params_bool_6, params_analysis_t_6, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins, @@ -296,7 +299,7 @@ TEST(FTAUtils, Quantification) try { - file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; nbins = 15; uncertainty = true; @@ -312,10 +315,10 @@ TEST(FTAUtils, Quantification) params_int_7, params_bool_7, params_analysis_t_7, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins, @@ -337,7 +340,7 @@ TEST(FTAUtils, Quantification) try { - file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; nbins = 15; uncertainty = true; @@ -352,10 +355,10 @@ TEST(FTAUtils, Quantification) params_int_8, params_bool_8, params_analysis_t_8, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins, @@ -381,7 +384,7 @@ TEST(FTAUtils, Quantification) // ++++++++++ Check FTA top event risk ++++++++++ - file_lists_quantification_utils = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; + file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; nbins = 15; FTAUtils::Quantification(params_double_9, @@ -389,10 +392,10 @@ TEST(FTAUtils, Quantification) params_int_9, params_bool_9, params_analysis_t_9, - file_lists_quantification_utils[0], - file_lists_quantification_utils[1], + file_lists_q[0], + file_lists_q[1], FTAUtils::Quantification::RISK, - file_lists_quantification_utils[2], + file_lists_q[2], IM[0], IM[1], nbins); From 191a4e516e727059fe51422844f58814bea503d7 Mon Sep 17 00:00:00 2001 From: HugoNip Date: Sat, 28 Mar 2020 23:14:21 -0400 Subject: [PATCH 04/10] remove all comments --- include/utils/FaultTreeUtils.h | 93 +++++------------------------ include/utils/QuantificationUtils.h | 89 ++++++--------------------- src/utils/FaultTreeUtils.C | 60 ++----------------- src/utils/QuantificationUtils.C | 52 ++-------------- unit/src/TestFaultTreeUtils.C | 28 +-------- unit/src/TestQuantificationUtils.C | 80 ++----------------------- 6 files changed, 50 insertions(+), 352 deletions(-) diff --git a/include/utils/FaultTreeUtils.h b/include/utils/FaultTreeUtils.h index 57df9b3a81..302b11edf4 100644 --- a/include/utils/FaultTreeUtils.h +++ b/include/utils/FaultTreeUtils.h @@ -18,35 +18,22 @@ class Quantification; class FaultTree; class CException; -/** - * String trim method which removes all leading and lagging whitespace in a string - */ + std::string trim(const std::string & str); -/** - * Converts ASCII string to upper case - */ + std::string str2Upper(const std::string & str_in, bool trim_input = false); -/** - * Returns interpolated value at x from parallel arrays ( xData, yData ) - * Assumes that xData has at least two elements, is sorted and is strictly monotonic increasing - * boolean argument extrapolate determines behaviour beyond ends of array (if needed) - */ double interpolate(vector_double data, double x, bool extrapolate); -/** - * Phi(-∞, x) aka N(x) - */ + double normalCDF(double x); double normalCDF(double x, double mu, double sigma); double Clip(double a, double min, double max); std::vector genQuantificationRVec(double dataPoint, int n, std::vector rv); -} // namespace FTAUtils - -/********************* Fault Tree Definition *********************/ +} class FTAUtils::FaultTree { @@ -70,115 +57,63 @@ class FTAUtils::FaultTree } }; - /** - * Constructor for fault tree class - * Default value of root is the first node in file - */ + FaultTree(std::string file_name, std::string root = ""); FaultTree(set_string & sets_link, std::map _node_base); - /** - * Destructor - */ + ~FaultTree(); std::string getRoot(); - /** - * Function returns cut sets at the given point - * NOTE - * If MOCUS ran before this function call, min cut sets will be returned - */ + set_string getCutSets(); - /** - * Builds m-ary fault tree - */ + std::map buildTree(Parser parser); - /** - * Computes minimum cut sets based on MOCUS Algorithm - */ set_string computeMinimumCutSets(); std::vector event(std::map); private: - // Hash map for operators std::map _opDict = {{"AND", AND}, {"OR", OR}}; - // Inverse mapping for printing purpose only std::map<_operator_t, std::string> _opDictInv = {{AND, "AND"}, {OR, "OR"}}; std::map _node_d_b; std::string _root; - // Cut sets container for in-place computations set_string _sets; - /** - * Translates string to opeartor - */ _operator_t str2Operator(std::string op); void rmSets(); _node * getNode(std::string name); - /** - * Recursive call function to flood fill sets by expanding them based on - * operation ALGO - * 1. Iterate through entire list and match for own node's name - * 2. Replace self with children based on operation - * (i) . Replace self with children in same row if AND - * (ii). Replace self with child one per row if OR - * 3. Recurse on each non leaf child - * - * NOTE - * 1. Updates "sets" variable - * 2. Uses std::set 2d array, hence absorption and idempotence properties - * are implicit - */ + void cutSetsExpand(_node * node); void removeSubsets(); }; -/************************** Parser Definition **************************/ class FTAUtils::Parser { public: - // Supported formats for parsing enum parseFormatT { FORMAT_CSV, FORMAT_UNDEF }; - /** - * Constructor for parser class - */ + Parser(std::string fileName, parseFormatT format); - /** - * Destructor for parser class - */ ~Parser(); - /** - * Yields all records, populates the standard structure and returns - * This function acts as an abstract layer to hide different formats - * that might be supported in future - * - * Returns: Array of strings - */ + vector_string yieldLines(); - /** - * Yields a single record, populates the standard structure and returns - * This function acts as an abstract layer to hide different formats - * that might be supported in future - * - * Returns: Array of strings - */ + std::vector yieldLine(); private: - // Handle to file + std::ifstream * fileP; }; @@ -189,4 +124,4 @@ class FTAUtils::CException CException(std::string s) : msg(s) {} }; -#endif // _FAULT_TREE_H +#endif diff --git a/include/utils/QuantificationUtils.h b/include/utils/QuantificationUtils.h index ad72c5a7bd..8f26c6fa5e 100644 --- a/include/utils/QuantificationUtils.h +++ b/include/utils/QuantificationUtils.h @@ -5,16 +5,15 @@ #include "FaultTreeUtils.h" #include "MastodonUtils.h" -/********************* Quantification Definition *********************/ class FTAUtils::Quantification { public: enum _dist_t { - PE, // = 0 - EXP, // = 1 - GAMMA, // = 2 + PE, + EXP, + GAMMA, WEIBULL, EXTREME, NORM, @@ -27,8 +26,8 @@ class FTAUtils::Quantification enum _analysis_t { - FRAGILITY, // = 0 - RISK // = 1 + FRAGILITY, + RISK }; Quantification(std::map & params_double, @@ -64,7 +63,6 @@ class FTAUtils::Quantification std::vector vec; vec.clear(); - // Remove all INF from vector for (int index = 0; index < vec_in.size(); index++) { if (!std::isinf(vec_in[index])) @@ -75,7 +73,6 @@ class FTAUtils::Quantification int size = vec.size(); - // ASSERT if (!(size > 0)) { fprintf(stderr, @@ -89,8 +86,6 @@ class FTAUtils::Quantification _pe = vec[0]; - // TODO: Remove 0th element from further calculations - // vec.erase( vec.begin() ); sort(vec.begin(), vec.end()); _mean = accumulate(vec.begin(), vec.end(), 0.0) / size; @@ -119,8 +114,7 @@ class FTAUtils::Quantification } }; - // TODO: 3 of the following are not supported as they need 1 arg rather than 2 - // for their distribution. Remove these 3 if not needed else add support + std::map _str2dist = { {"GAMMA", GAMMA}, {"WEIBULL", WEIBULL}, @@ -135,11 +129,7 @@ class FTAUtils::Quantification vector_double _cut_set_prob; std::map> _b_nodes; - /** - * Generates probability vector for a specified distribution - * Nomenclature of a and b changes with distribution. - * eg., a => mean, b => std for NORM - */ + std::vector getProbVector(_dist_t dist, double a, double b, @@ -149,9 +139,7 @@ class FTAUtils::Quantification _analysis_t analysis, bool uc); - /** - * Parses and floods probabilties of basic elements - */ + vector_string beProb(FTAUtils::Parser parser, int n_sample, int seed, @@ -159,10 +147,7 @@ class FTAUtils::Quantification std::vector intmes, bool uncert); - /** - * Computes probability for a given cut set on a per set basis based on - * pre-flooded basic elem probs - */ + vector_double computeCutSetProb(std::set> cut_sets, int n, bool bypass = false, @@ -176,76 +161,40 @@ class FTAUtils::Quantification std::string bypass_key = "", double bypass_value = 0); - /** - * Computes accumulated probability for the entire cut set - * 3 Digests are computed: - * 0. Min Max - * 1. Upper Bound - * 2. Top Rare - * - * NOTE - * 1. Its better to compute these 2 together as they have common loops - * which can save time - * 2. Sets the vector cutSetProb - */ + std::vector * cutSetProbWDigest(std::set> cut_sets, int n, bool only_min_max = false); - /** - * Wrapper function for gate probability arith - * AND => a * b - * OR => (1 - a) * (1 - b) - */ + double getGateProb(double a, double b, bool is_and); - /** - * Computes cut-set Fussel-Vesely Important measures - */ + vector_double minCutIM(std::vector upperBound); - /** - * For each basic element: - * 1. Calculate #occurrence in minimal cut set (beCount) - * 2. Store row index of all occurrence (beIndex) - */ + std::map beIM(std::set> cut_sets, int n, std::vector upper_bound, std::vector & count_v); - /** - * Translator function for string -> double - */ + vector_double linesToDouble(vector_string lines); - /** - * Function for getting BIN mean values of intensity measure - */ + std::vector getBinMeans(double im_lower, double im_upper, int n_bins); - /** - * Function for interpolating the hazard curve based on range of intensity - * measures - */ + std::vector hazInterp(vector_double hazard, std::vector im_bins); - /** - * Function for top event fragility - */ + std::vector fragility(std::set> cut_sets, int n, std::vector im_bins, double & mu, double & sigma); - /** - * Function for calculating risk by convoluting hazard and fragility - * Assign risk to basic event probabilites - * - * WARNING - * This consumes _b_nodes and then overwrites it - */ + void computeRisk(int n, std::vector hazard); }; -#endif // _QUANTIFICATION_H +#endif diff --git a/src/utils/FaultTreeUtils.C b/src/utils/FaultTreeUtils.C index 2f79347af4..34a8276814 100644 --- a/src/utils/FaultTreeUtils.C +++ b/src/utils/FaultTreeUtils.C @@ -8,10 +8,8 @@ FTAUtils::FaultTree::FaultTree(std::string file_name, std::string root) FTAUtils::Parser parser = FTAUtils::Parser(file_name, FTAUtils::Parser::FORMAT_CSV); buildTree(parser); - // Override root node if not default if (root != "") { - // ASSERT if (!(getNode(root))) { fprintf(stderr, @@ -36,21 +34,16 @@ FTAUtils::FaultTree::FaultTree(set_string & sets_link, std::map root_set, sub_sets; - root_set.insert("TOP"); // "TOP", size: 1 + root_set.insert("TOP"); _sets.insert(root_set); - // Start with root node to expand on heirarchy cutSetsExpand(getNode("TOP")); - // Check if first row is empty - // NOTE: Since its std::set, only first row could - // be empty + if (_sets.begin()->size() == 0) { _sets.erase(_sets.begin()); @@ -71,22 +64,17 @@ FTAUtils::FaultTree::buildTree(FTAUtils::Parser parser) { line = parser.yieldLine(); - // Stop if no new line to process if (line.size() == 0) break; - // Stash name, operator _node * node = new _node(line[0], str2Operator(line[1])); - // Add children for (int i = 2; i < line.size(); i++) node->_child.push_back(line[i]); - // Stash the first entry as ROOT of tree if (_node_d_b.size() == 0) _root = line[0]; - // Add the newly created node to node lookup hashmap _node_d_b[line[0]] = node; } @@ -97,7 +85,6 @@ FTAUtils::FaultTree::_operator_t FTAUtils::FaultTree::str2Operator(std::string op) { std::string op_s = FTAUtils::str2Upper(op, true); - // ASSERT if (_opDict.count(op_s) == 0) { fprintf(stderr, @@ -114,21 +101,16 @@ FTAUtils::FaultTree::str2Operator(std::string op) set_string FTAUtils::FaultTree::computeMinimumCutSets() { - // Clearing up sets vector for safety rmSets(); - // Add root to set std::set root_set, sub_sets; root_set.insert(getRoot()); _sets.insert(root_set); - // Start with root node to expand on heirarchy cutSetsExpand(getNode(getRoot())); - // Check if first row is empty - // NOTE: Since its std::set, only first row could - // be empty + if (_sets.begin()->size() == 0) { _sets.erase(_sets.begin()); @@ -154,10 +136,8 @@ FTAUtils::FaultTree::removeSubsets() includes(row2_it->begin(), row2_it->end(), row1_it->begin(), row1_it->end()); bool row2_is_subset = includes(row1_it->begin(), row1_it->end(), row2_it->begin(), row2_it->end()); - // Remove row1 if row2 is a subset of row1 by marking it if (row2_is_subset) rm_its.insert(index1); - // Remove row2 if row1 is a subset of row1 by marking it else if (row1_is_subset) rm_its.insert(index2); index2++; @@ -166,7 +146,6 @@ FTAUtils::FaultTree::removeSubsets() } uint64_t index = 0; - // Remove rows marked to be removed for (std::set::iterator it = rm_its.begin(); it != rm_its.end(); ++it) { _sets.erase(next(_sets.begin(), *it - index)); @@ -177,7 +156,6 @@ FTAUtils::FaultTree::removeSubsets() void FTAUtils::FaultTree::cutSetsExpand(_node * node) { - // ASSERT if (!node) { fprintf(stderr, @@ -188,7 +166,6 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) abort(); } - // ASSERT if (node->_child.size() == 0) { fprintf(stderr, @@ -200,28 +177,21 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) abort(); } - // New rows which have to be appended at end. Adding them might - // result in data hazards + set_string mk_rows; - // 1. Iterate through entire list and match for own node's name for (set_string::iterator row_it = _sets.begin(); row_it != _sets.end(); ++row_it) { std::set row(*row_it); - // Search for the element in row std::set::iterator it = find(row.begin(), row.end(), node->_name); if (it != row.end()) { - // Erase self to add children row.erase(it); - // 2. Replace self with children based on operation - // (i) . Replace self with children in same row if AND - // (ii). Replace self with child one per row if OR + switch (node->_op) { case AND: - // Append all children in same row for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) { row.insert(node->_child[c_id]); @@ -229,7 +199,6 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) break; case OR: - // Replicate line and append one child per row for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) { std::set row_set(row); @@ -237,13 +206,11 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) mk_rows.insert(row_set); } - // Clear this row to be removed at end (postprocessing) row.clear(); break; default: { - // ASSERT fprintf(stderr, "[ASSERT] In File: %s, Line: %d => " "Unknown Operator found!.\n", @@ -255,22 +222,17 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) } } - // The following 2 lines might result in an empty row to be pushed. - // If we gate insertion with count, we might have a data hazard due - // to wrong increment of iterators. - // A quick fix to this is to check for first row and delete it at end + _sets.erase(row_it); _sets.insert(row); } - // Add newly created rows for (set_string::iterator it = mk_rows.begin(); it != mk_rows.end(); ++it) { std::set row(*it); _sets.insert(*it); } - // 3. Recurse on each non leaf child for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) { _node * child = getNode(node->_child[c_id]); @@ -308,7 +270,6 @@ FTAUtils::FaultTree::getCutSets() FTAUtils::Parser::Parser(std::string fileName, FTAUtils::Parser::parseFormatT format) { - // Assertion on supported parsing formats if (format != FORMAT_CSV) { fprintf(stderr, @@ -337,7 +298,6 @@ FTAUtils::Parser::yieldLines() { line = FTAUtils::Parser::yieldLine(); - // Stop if no new line to process if (line.size() == 0) break; lines.push_back(line); @@ -348,7 +308,6 @@ FTAUtils::Parser::yieldLines() std::vector FTAUtils::Parser::yieldLine() { - // ASSERT if (!(fileP->is_open())) { fprintf(stderr, @@ -362,7 +321,6 @@ FTAUtils::Parser::yieldLine() std::string buffer; std::vector line; - // Get a line that has something except \n if (getline(*fileP, buffer)) { std::string token; @@ -406,10 +364,8 @@ FTAUtils::interpolate(vector_double data, double x, bool extrapolate) { int size = data.size(); - // find left end of interval for interpolation int i = 0; - // special case: beyond right end if (x >= data[size - 2][0]) { i = size - 2; @@ -420,10 +376,8 @@ FTAUtils::interpolate(vector_double data, double x, bool extrapolate) i++; } - // points on either side (unless beyond ends) double xL = data[i][0], yL = data[i][1], xR = data[i + 1][0], yR = data[i + 1][1]; - // if beyond ends of array and not extrapolating if (!extrapolate) { if (x < xL) @@ -432,10 +386,8 @@ FTAUtils::interpolate(vector_double data, double x, bool extrapolate) yL = yR; } - // gradient double dydx = (yR - yL) / (xR - xL); - // linear interpolation return yL + dydx * (x - xL); } diff --git a/src/utils/QuantificationUtils.C b/src/utils/QuantificationUtils.C index e740a5d4fe..4c578b0cb8 100644 --- a/src/utils/QuantificationUtils.C +++ b/src/utils/QuantificationUtils.C @@ -17,7 +17,6 @@ FTAUtils::Quantification::Quantification(std::map & int n_sample, int seed) { - //---------------- SAVE PARAMETERS ------------------------------------------ std::vector params_double1; params_double1.push_back(im_lower); params_double1.push_back(im_upper); @@ -31,55 +30,44 @@ FTAUtils::Quantification::Quantification(std::map & params_bool.insert(std::pair("uncertainty", uncertainty)); params_analysis_t.insert(std::pair("analysis", analysis)); - //---------------- COMPUTE -------------------------------------------------- - // Construct the fault tree with MOCUS algo for minimal cut sets + FTAUtils::FaultTree ft = FTAUtils::FaultTree(events_file, root); Parser parser_events = Parser(events_file, Parser::FORMAT_CSV); vector_string lines_events = parser_events.yieldLines(); set_string cut_sets = ft.getCutSets(); params_string.insert(std::pair("events_files", lines_events)); - // Read basic events file FTAUtils::Parser parser_event_prob = FTAUtils::Parser(events_prob_file, FTAUtils::Parser::FORMAT_CSV); vector_string lines_events_prob; if (analysis == FRAGILITY) { - // Read and interpolate hazard curve FTAUtils::Parser parser_hazard = FTAUtils::Parser(hazard_file, FTAUtils::Parser::FORMAT_CSV); vector_string lines_hazard = parser_hazard.yieldLines(); vector_double hazard = FTAUtils::Quantification::linesToDouble(lines_hazard); params_double.insert(pair_string_double("hazard", hazard)); - // List of Intensity Measure bin values std::vector im_bins = getBinMeans(im_lower, im_upper, n_bins); - // List of Intensity Measure bin values std::vector hazard_freq = hazInterp(hazard, im_bins); - // Dictionary of basic events for fragility input lines_events_prob = beProb(parser_event_prob, n_sample, seed, analysis, im_bins, uncertainty); - // Top event fragility (lognormal parameters) double mu, sigma; fragility(cut_sets, n_bins, im_bins, mu, sigma); - // Dictionary of basic events risk (convoluting fragility and hazard) computeRisk(n_bins, hazard_freq); } else { - // Dictionary of basic events risk (risk inputs) lines_events_prob = beProb(parser_event_prob, n_sample, seed, analysis, std::vector(), uncertainty); } params_string.insert(std::pair("basic_events", lines_events_prob)); - // Calculate minimal cut sets probability - // Digest cut set probability to unified probability - // Adding 1 to accomodate mean as 0th element + std::vector * digest = cutSetProbWDigest(cut_sets, n_sample + 1); vector_double mc_im = minCutIM(digest[1]); Stats s0(digest[0]); @@ -88,7 +76,6 @@ FTAUtils::Quantification::Quantification(std::map & std::vector count_v; std::map be_stats = beIM(cut_sets, n_sample + 1, digest[1], count_v); - // return the FTA top event risk std::vector fta_0; vector_double fta_1; fta_0.push_back(s0._pe); @@ -97,12 +84,10 @@ FTAUtils::Quantification::Quantification(std::map & fta_1.push_back(fta_0); params_double.insert(std::make_pair("fta", fta_1)); - // Return the Risk Reduction Ratio for B1, B2, B3, B4, B5 params_double.insert(pair_string_double("fv", be_stats["fv"])); params_double.insert(pair_string_double("rrr", be_stats["rrr"])); params_double.insert(pair_string_double("rir", be_stats["rir"])); - // Return the Risk Reduction Difference for B1, B2, B3, B4, B5 params_double.insert(pair_string_double("rri", be_stats["rri"])); params_double.insert(pair_string_double("rii", be_stats["rii"])); params_double.insert(pair_string_double("bi", be_stats["bi"])); @@ -112,10 +97,8 @@ std::vector FTAUtils::Quantification::hazInterp(vector_double hazard, std::vector im_bins) { std::vector data; - // Convert all hazard values to log10 for (int index = 0; index < hazard.size(); index++) { - // ASSERT if (hazard[index].size() != 2) { fprintf(stderr, @@ -130,7 +113,6 @@ FTAUtils::Quantification::hazInterp(vector_double hazard, std::vector im hazard[index][1] = log10(hazard[index][1]); } - // Interpolate for (int index = 0; index < im_bins.size(); index++) { data.push_back(pow(10, interpolate(hazard, log10(im_bins[index]), true))); @@ -178,7 +160,6 @@ FTAUtils::genQuantificationRVec(double dataPoint, int n, std::vector rv) { for (int index = 0; index < n; index++) { - // double dataPoint = double(gen); rv.push_back(FTAUtils::Clip(dataPoint, 0.0, 1.0)); } return rv; @@ -200,7 +181,6 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, if (analysis == FRAGILITY) { - // ASSERT if (dist != LNORM) { fprintf(stderr, @@ -234,7 +214,6 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, break; default: { - // ASSERT fprintf(stderr, "[ASSERT] In File: %s, Line: %d => " "Un-supported dist found.\n", @@ -265,7 +244,6 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, break; default: { - // ASSERT fprintf(stderr, "[ASSERT] In File: %s, Line: %d => " "Un-supported dist found.\n", @@ -293,13 +271,11 @@ FTAUtils::Quantification::beProb(FTAUtils::Parser parser, { line = parser.yieldLine(); - // Stop if no new line to process if (line.size() != 0) lines.push_back(line); else break; - // Stash name, probability vector double b = line.size() > 3 ? stod(line[3]) : 0; _b_nodes[line[0]] = getProbVector( _str2dist[line[1]], stod(line[2]), b, n_sample, seed, intmes, analysis, uncert); @@ -313,17 +289,14 @@ FTAUtils::Quantification::cutSetProbWDigest(set_string cut_sets, int n, bool onl { std::vector * digest = new std::vector[3]; - // Digest 2: Generate power set to generate all possible combinations vector_double ps_p; int pow_set_size = pow(2, cut_sets.size()); - // Start from 1 as NULL set is not needed for (int index = 1; index < pow_set_size; index++) { std::set power_set_row; for (int j = 0; j < cut_sets.size(); j++) { - // Check if jth bit in the index is set - // If set then form row of power set + if (index & (1 << j)) { std::set j_elem = *next(cut_sets.begin(), j); @@ -333,15 +306,12 @@ FTAUtils::Quantification::cutSetProbWDigest(set_string cut_sets, int n, bool onl } } } - // For elements taken n at a time, if n is odd, positive else negative bool is_positive = (__builtin_popcount(index) % 2); std::vector cut_prob = cutSetRowProb(power_set_row, n, is_positive); ps_p.push_back(cut_prob); } - // Digest 0, 1, 2 - // Avoiding an extra function call to interm event like in the python - // counterpart for performance reasons (loop merging) + _cut_set_prob = computeCutSetProb(cut_sets, n); for (int index = 0; index < n; index++) { @@ -372,9 +342,7 @@ vector_double FTAUtils::Quantification::computeCutSetProb( set_string cut_sets, int n, bool bypass, std::string bypass_key, double bypass_value) { - // For each row in cut set, compute: - // 1. For each sample in the generated vector, compute AND gate - // probability analysis for each basic element + vector_double quant; for (set_string::iterator row = cut_sets.begin(); row != cut_sets.end(); ++row) { @@ -398,7 +366,6 @@ FTAUtils::Quantification::cutSetRowProb(std::set row, double value = 1; for (std::set::iterator col = row.begin(); col != row.end(); ++col) { - // ASSERT if (_b_nodes.find(*col) == _b_nodes.end()) { fprintf(stderr, @@ -460,14 +427,9 @@ std::vector FTAUtils::Quantification::fragility( set_string cut_sets, int n, std::vector im_bins, double & mu, double & sigma) { - // 1. Calculate TOP event fragility using min-max approach (exact) {digest[0]} std::vector * digest = cutSetProbWDigest(cut_sets, n, true); - // 2. lognormal parameters of TOP Event fragility - // solveLnParams(im_bins, digest[0], mu, sigma); - // Use the maximizelikehood function instead of solver function - // Use SGD for optimization std::vector loc_space = {0.01, 2}; std::vector sca_space = {0.01, 1}; std::vector max_values = MastodonUtils::maximizeLogLikelihood( @@ -475,7 +437,6 @@ FTAUtils::Quantification::fragility( mu = max_values[0]; sigma = max_values[1]; - // TODO: 3. TOP event Fragility plot return digest[0]; } @@ -490,7 +451,6 @@ FTAUtils::Quantification::beIM(set_string cut_sets, bn_it != _b_nodes.end(); ++bn_it) { - // Generate available vector to save computation on per index loop std::vector available; int count = 0; for (set_string::iterator cs_it = cut_sets.begin(); cs_it != cut_sets.end(); ++cs_it) @@ -507,7 +467,6 @@ FTAUtils::Quantification::beIM(set_string cut_sets, vector_double mc_p1 = computeCutSetProb(cut_sets, n, true, bn_it->first, 1); vector_double mc_p0 = computeCutSetProb(cut_sets, n, true, bn_it->first, 0); - // OR it with loop merging for (int index = 0; index < n; index++) { double f0_val = 1; @@ -531,7 +490,6 @@ FTAUtils::Quantification::beIM(set_string cut_sets, rii.push_back(f1_val - upper_bound[index]); bi.push_back(f1_val - f0_val); } - // Calculate stats for current row and stash Stats fv_stats(fv); Stats rrr_stats(rrr); Stats rri_stats(rri); diff --git a/unit/src/TestFaultTreeUtils.C b/unit/src/TestFaultTreeUtils.C index ccfc21d31a..9db7f291e5 100644 --- a/unit/src/TestFaultTreeUtils.C +++ b/unit/src/TestFaultTreeUtils.C @@ -1,71 +1,58 @@ -// MOOSE includes #include "gtest/gtest.h" #include "MooseUtils.h" #include "Conversion.h" -// Boost distribution includes #include "BoostDistribution.h" -// MASTODON includes #include "MastodonUtils.h" -// QUANTIFICATIONUTILS includes #include "FaultTreeUtils.h" typedef FTAUtils::FaultTree::_node ft_node_test; -// Helper to make sure the test tree is correct. void assertTree(std::map tree_node) { - // TOP EXPECT_EQ(tree_node["TOP"]->_name, "TOP"); EXPECT_EQ(tree_node["TOP"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["TOP"]->_child.at(0), "GATE1"); EXPECT_EQ(tree_node["TOP"]->_child.at(1), "GATE2"); - // GATE1 EXPECT_EQ(tree_node["GATE1"]->_name, "GATE1"); EXPECT_EQ(tree_node["GATE1"]->_op, FTAUtils::FaultTree::OR); EXPECT_EQ(tree_node["GATE1"]->_child.at(0), "FT-N/m-1"); EXPECT_EQ(tree_node["GATE1"]->_child.at(1), "FT-N/m-2"); EXPECT_EQ(tree_node["GATE1"]->_child.at(2), "FT-N/m-3"); - // GATE2 EXPECT_EQ(tree_node["GATE2"]->_name, "GATE2"); EXPECT_EQ(tree_node["GATE2"]->_op, FTAUtils::FaultTree::OR); EXPECT_EQ(tree_node["GATE2"]->_child.at(0), "B1"); EXPECT_EQ(tree_node["GATE2"]->_child.at(1), "B3"); EXPECT_EQ(tree_node["GATE2"]->_child.at(2), "B4"); - // FT-N/m-1 EXPECT_EQ(tree_node["FT-N/m-1"]->_name, "FT-N/m-1"); EXPECT_EQ(tree_node["FT-N/m-1"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(0), "GATE3"); EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(1), "B3"); EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(2), "B5"); - // FT-N/m-2 EXPECT_EQ(tree_node["FT-N/m-2"]->_name, "FT-N/m-2"); EXPECT_EQ(tree_node["FT-N/m-2"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["FT-N/m-2"]->_child.at(0), "GATE3"); EXPECT_EQ(tree_node["FT-N/m-2"]->_child.at(1), "B1"); - // FT-N/m-3 EXPECT_EQ(tree_node["FT-N/m-3"]->_name, "FT-N/m-3"); EXPECT_EQ(tree_node["FT-N/m-3"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(0), "B3"); EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(1), "B5"); EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(2), "B1"); - // GATE3 EXPECT_EQ(tree_node["GATE3"]->_name, "GATE3"); EXPECT_EQ(tree_node["GATE3"]->_op, FTAUtils::FaultTree::OR); EXPECT_EQ(tree_node["GATE3"]->_child.at(0), "B2"); EXPECT_EQ(tree_node["GATE3"]->_child.at(1), "B4"); } -// Helper to make sure the mocus is correct. void assertMocus(set_string tree_mocus) { @@ -92,7 +79,6 @@ assertMocus(set_string tree_mocus) } } -// Create a test event list. std::map eventList() { @@ -118,26 +104,19 @@ eventList() for (int i = 0; i < line.size(); i++) { ft_node_test * node = new ft_node_test(line[i][0], op[i]); - // Add children for (int j = 2; j < line[i].size(); j++) node->_child.push_back(line[i][j]); - // Add the newly created node to node lookup hashmap _node_d_b[line[i][0]] = node; } return _node_d_b; } -// For input std::vector file_lists; TEST(FTAUtils, FaultTree) { - // ==================== TestCase for FaultTree Object =================== - // ---------------- Test error message for events input ---------------- - - // ###### File Inputs ###### try { @@ -147,14 +126,10 @@ TEST(FTAUtils, FaultTree) } catch (FTAUtils::CException e) { - /* - * ------- IO Error ------- - * does not exist - */ + EXPECT_EQ(e.msg, "Unable to open file."); } - // ---------------- Test creating tree from list ---------------- std::map tree_node_from_list = eventList(); set_string tree_mocus_from_list; @@ -162,7 +137,6 @@ TEST(FTAUtils, FaultTree) assertTree(tree_node_from_list); assertMocus(tree_mocus_from_list); - // ---------------- Test creating tree from file ---------------- file_lists = {"logic2.txt", ""}; diff --git a/unit/src/TestQuantificationUtils.C b/unit/src/TestQuantificationUtils.C index 2295e91c84..9920d089f1 100644 --- a/unit/src/TestQuantificationUtils.C +++ b/unit/src/TestQuantificationUtils.C @@ -1,18 +1,13 @@ -// MOOSE includes #include "gtest/gtest.h" #include "MooseUtils.h" #include "Conversion.h" -// Boost distribution includes #include "BoostDistribution.h" -// MASTODON includes #include "MastodonUtils.h" -// QUANTIFICATIONUTILS includes #include "QuantificationUtils.h" -// For input std::vector file_lists_q; std::vector IM; int nbins; @@ -20,7 +15,6 @@ bool uncertainty; int nsamp; int seed; -// For output std::map>> params_double_1, params_double_2, params_double_3, params_double_4, params_double_5, params_double_6, params_double_7, params_double_8, params_double_9; @@ -39,18 +33,9 @@ std::map params_analysis_t_1 std::map params_int_1, params_int_2, params_int_3, params_int_4, params_int_5, params_int_6, params_int_7, params_int_8, params_int_9; -// TestCase for Quantification Object TEST(FTAUtils, Quantification) { - // =========================================================================== - // =============================== Test Inputs =============================== - // =========================================================================== - - // +++++++++++++++++++++++ Input values +++++++++++++++++++++++ - - // ############### File Inputs ############### - try { file_lists_q = { @@ -73,10 +58,7 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - /* - * ------- IO Error ------- - * does not exist - */ + EXPECT_EQ(e.msg, "Unable to open file."); } @@ -101,14 +83,10 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - /* - * ------- Type Error ------- - * must be a filename or a list - */ + EXPECT_EQ(e.msg, "Unable to open file."); } - // ############## Inputs for Fragility ############## file_lists_q = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; IM = {0.1, 4}; @@ -126,7 +104,6 @@ TEST(FTAUtils, Quantification) IM[1], nbins); - // >>>>>>>> logic Value Check <<<<<<<< std::vector> event_files_matrix3{{"TE", "OR", "IE3", "IE4"}, {"IE4", "OR", "C4"}, @@ -135,7 +112,6 @@ TEST(FTAUtils, Quantification) {"IE1", "OR", "C1"}}; EXPECT_EQ(params_string_3["events_files"], event_files_matrix3); - // >>>>>>>> Basic Events Value Check <<<<<<<< std::vector> basic_events_matrix3{{"C1", "LNORM", "1.88", "0.5"}, {"C2", "LNORM", "3.78", "0.79"}, @@ -144,11 +120,9 @@ TEST(FTAUtils, Quantification) std::vector> basic_events_3 = params_string_3["basic_events"]; EXPECT_EQ(basic_events_3, basic_events_matrix3); - // >>>>>>>> antype Value Check <<<<<<<< EXPECT_EQ(params_analysis_t_3["analysis"], FTAUtils::Quantification::FRAGILITY); - // >>>>>>>> hazard Value Check <<<<<<<< std::vector> matrix_hazard{{0.0608, 0.01}, {0.2124, 0.001}, @@ -159,16 +133,13 @@ TEST(FTAUtils, Quantification) std::vector> hazard = params_double_3["hazard"]; EXPECT_EQ(hazard, matrix_hazard); - // >>>>>>>> imrange Value Check <<<<<<<< EXPECT_EQ(params_double_3["IM"][0][0], 0.1); EXPECT_EQ(params_double_3["IM"][0][1], 4); - // >>>>>>>> nbins Value Check <<<<<<<< EXPECT_EQ(params_int_3["n_bins"], 15); - // >>>>>>>> nbins Value Check <<<<<<<< try { @@ -194,14 +165,10 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - /* - * ------- Value Error ------- - * The supplied value of nbins must be a +ve integer - */ + EXPECT_EQ(e.msg, "ValueError"); } - // ############## Inputs for Risk analysis (not fragility) ############## file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt"}; @@ -214,7 +181,6 @@ TEST(FTAUtils, Quantification) file_lists_q[1], FTAUtils::Quantification::RISK); - // >>>>>>>> logic Value Check <<<<<<<< std::vector> event_files_matrix5{ {"TOP", "AND", "GATE1", "GATE2"}, @@ -227,7 +193,6 @@ TEST(FTAUtils, Quantification) EXPECT_EQ(params_string_5["events_files"], event_files_matrix5); - // >>>>>>>> Basic Events Value Check <<<<<<<< std::vector> basic_events_matrix5{{"B1", "PE", "0.01"}, {"B2", "PE", "0.02"}, @@ -239,23 +204,18 @@ TEST(FTAUtils, Quantification) EXPECT_EQ(basic_events_5, basic_events_matrix5); - // >>>>>>>> antype Value Check <<<<<<<< EXPECT_NE(params_analysis_t_5["analysis"], FTAUtils::Quantification::FRAGILITY); - // >>>>>>>> uncertainty Value Check <<<<<<<< EXPECT_EQ(params_bool_5["uncertainty"], false); - // >>>>>>>> nsamp Value Check <<<<<<<< EXPECT_EQ(params_int_5["nsamp"], 1); - // >>>>>>>> seed Value Check <<<<<<<< EXPECT_EQ(params_int_5["seed"], 0); - // ###### Testing for input errors making sure parameters are input correctly ###### file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; @@ -281,21 +241,16 @@ TEST(FTAUtils, Quantification) nsamp, seed); - // >>>>>>>> uncertainty Value Check <<<<<<<< EXPECT_EQ(params_bool_6["uncertainty"], true); - // >>>>>>>> nsamp Value Check <<<<<<<< EXPECT_EQ(params_int_6["nsamp"], 1000); - // >>>>>>>> seed Value Check <<<<<<<< EXPECT_EQ(params_int_6["seed"], 436546754); - // ############## Testing for input type errors and value errors ############## - // >>>>>>>> nsamp Value Check <<<<<<<< try { @@ -329,14 +284,10 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - /* - * ------- Value Error ------- - * The supplied value of nsamp must be a +ve integer. - */ + EXPECT_EQ(e.msg, "ValueError"); } - // >>>>>>>> seed Value Check <<<<<<<< try { @@ -369,20 +320,11 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - /* - * ------- Value Error ------- - * The supplied value of seed must be a +ve integer. - */ + EXPECT_EQ(e.msg, "ValueError"); } - // =========================================================================== - // ============================== Test TOP Risk ============================== - // =========================================================================== - - // ############ Function for asserting FTA top event risk ############ - // ++++++++++ Check FTA top event risk ++++++++++ file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; @@ -400,57 +342,45 @@ TEST(FTAUtils, Quantification) IM[1], nbins); - // >>>>>>>> Risk Value Check <<<<<<<< - // min-max EXPECT_EQ(std::to_string(params_double_9["fta"][0][0]), "0.000694"); - // upper bound EXPECT_EQ(std::to_string(params_double_9["fta"][0][1]), "0.000705"); - // rare event EXPECT_EQ(std::to_string(params_double_9["fta"][0][2]), "0.000705"); - // >>>>>>>> IMratio Value Check <<<<<<<< - // Fussell-Vesely for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["fv"][0][0]), "0.872395"); EXPECT_EQ(std::to_string(params_double_9["fv"][1][0]), "0.326300"); EXPECT_EQ(std::to_string(params_double_9["fv"][2][0]), "0.148963"); EXPECT_EQ(std::to_string(params_double_9["fv"][3][0]), "0.652584"); EXPECT_EQ(std::to_string(params_double_9["fv"][4][0]), "0.148963"); - // Risk Reduction Ratio for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rrr"][0][0]), "7.831866"); EXPECT_EQ(std::to_string(params_double_9["rrr"][1][0]), "1.483999"); EXPECT_EQ(std::to_string(params_double_9["rrr"][2][0]), "1.174913"); EXPECT_EQ(std::to_string(params_double_9["rrr"][3][0]), "2.877066"); EXPECT_EQ(std::to_string(params_double_9["rrr"][4][0]), "1.174913"); - // Risk Increase Ratio for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rir"][0][0]), "86.111103"); EXPECT_EQ(std::to_string(params_double_9["rir"][1][0]), "16.960273"); EXPECT_EQ(std::to_string(params_double_9["rir"][2][0]), "5.808755"); EXPECT_EQ(std::to_string(params_double_9["rir"][3][0]), "16.637742"); EXPECT_EQ(std::to_string(params_double_9["rir"][4][0]), "3.826894"); - // >>>>>>>> IMdiff Value Check <<<<<<<< - // Risk Reduction Difference for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rri"][0][0]), "0.000615"); EXPECT_EQ(std::to_string(params_double_9["rri"][1][0]), "0.000230"); EXPECT_EQ(std::to_string(params_double_9["rri"][2][0]), "0.000105"); EXPECT_EQ(std::to_string(params_double_9["rri"][3][0]), "0.000460"); EXPECT_EQ(std::to_string(params_double_9["rri"][4][0]), "0.000105"); - // Risk Increase Difference for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rii"][0][0]), "0.059991"); EXPECT_EQ(std::to_string(params_double_9["rii"][1][0]), "0.011250"); EXPECT_EQ(std::to_string(params_double_9["rii"][2][0]), "0.003389"); EXPECT_EQ(std::to_string(params_double_9["rii"][3][0]), "0.011022"); EXPECT_EQ(std::to_string(params_double_9["rii"][4][0]), "0.001993"); - // Birnbaum Difference for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["bi"][0][0]), "0.060606"); EXPECT_EQ(std::to_string(params_double_9["bi"][1][0]), "0.011480"); EXPECT_EQ(std::to_string(params_double_9["bi"][2][0]), "0.003494"); From a8204affbb8ee383490636822804287f4d1131bd Mon Sep 17 00:00:00 2001 From: Hugo Nie Date: Mon, 30 Mar 2020 15:09:36 -0400 Subject: [PATCH 05/10] Update FaultTreeUtils.C ref #126 1. make macro #define a method 2. remove all commented sections 3. remove banned keywords --- include/utils/FaultTreeUtils.h | 93 +++++++++++++++++++++++++++++----- 1 file changed, 79 insertions(+), 14 deletions(-) diff --git a/include/utils/FaultTreeUtils.h b/include/utils/FaultTreeUtils.h index 302b11edf4..57df9b3a81 100644 --- a/include/utils/FaultTreeUtils.h +++ b/include/utils/FaultTreeUtils.h @@ -18,22 +18,35 @@ class Quantification; class FaultTree; class CException; - +/** + * String trim method which removes all leading and lagging whitespace in a string + */ std::string trim(const std::string & str); - +/** + * Converts ASCII string to upper case + */ std::string str2Upper(const std::string & str_in, bool trim_input = false); +/** + * Returns interpolated value at x from parallel arrays ( xData, yData ) + * Assumes that xData has at least two elements, is sorted and is strictly monotonic increasing + * boolean argument extrapolate determines behaviour beyond ends of array (if needed) + */ double interpolate(vector_double data, double x, bool extrapolate); - +/** + * Phi(-∞, x) aka N(x) + */ double normalCDF(double x); double normalCDF(double x, double mu, double sigma); double Clip(double a, double min, double max); std::vector genQuantificationRVec(double dataPoint, int n, std::vector rv); -} +} // namespace FTAUtils + +/********************* Fault Tree Definition *********************/ class FTAUtils::FaultTree { @@ -57,63 +70,115 @@ class FTAUtils::FaultTree } }; - + /** + * Constructor for fault tree class + * Default value of root is the first node in file + */ FaultTree(std::string file_name, std::string root = ""); FaultTree(set_string & sets_link, std::map _node_base); - + /** + * Destructor + */ ~FaultTree(); std::string getRoot(); - + /** + * Function returns cut sets at the given point + * NOTE + * If MOCUS ran before this function call, min cut sets will be returned + */ set_string getCutSets(); - + /** + * Builds m-ary fault tree + */ std::map buildTree(Parser parser); + /** + * Computes minimum cut sets based on MOCUS Algorithm + */ set_string computeMinimumCutSets(); std::vector event(std::map); private: + // Hash map for operators std::map _opDict = {{"AND", AND}, {"OR", OR}}; + // Inverse mapping for printing purpose only std::map<_operator_t, std::string> _opDictInv = {{AND, "AND"}, {OR, "OR"}}; std::map _node_d_b; std::string _root; + // Cut sets container for in-place computations set_string _sets; + /** + * Translates string to opeartor + */ _operator_t str2Operator(std::string op); void rmSets(); _node * getNode(std::string name); - + /** + * Recursive call function to flood fill sets by expanding them based on + * operation ALGO + * 1. Iterate through entire list and match for own node's name + * 2. Replace self with children based on operation + * (i) . Replace self with children in same row if AND + * (ii). Replace self with child one per row if OR + * 3. Recurse on each non leaf child + * + * NOTE + * 1. Updates "sets" variable + * 2. Uses std::set 2d array, hence absorption and idempotence properties + * are implicit + */ void cutSetsExpand(_node * node); void removeSubsets(); }; +/************************** Parser Definition **************************/ class FTAUtils::Parser { public: + // Supported formats for parsing enum parseFormatT { FORMAT_CSV, FORMAT_UNDEF }; - + /** + * Constructor for parser class + */ Parser(std::string fileName, parseFormatT format); + /** + * Destructor for parser class + */ ~Parser(); - + /** + * Yields all records, populates the standard structure and returns + * This function acts as an abstract layer to hide different formats + * that might be supported in future + * + * Returns: Array of strings + */ vector_string yieldLines(); - + /** + * Yields a single record, populates the standard structure and returns + * This function acts as an abstract layer to hide different formats + * that might be supported in future + * + * Returns: Array of strings + */ std::vector yieldLine(); private: - + // Handle to file std::ifstream * fileP; }; @@ -124,4 +189,4 @@ class FTAUtils::CException CException(std::string s) : msg(s) {} }; -#endif +#endif // _FAULT_TREE_H From 0c8e507420b1ae95ea469bb76de341b1aad7a13f Mon Sep 17 00:00:00 2001 From: Hugo Nie Date: Mon, 30 Mar 2020 17:15:20 -0400 Subject: [PATCH 06/10] update quanticationUtils.h Ref #126 1. make macro #define a method 2. remove all commented sections 3. remove banned keywords --- include/utils/QuantificationUtils.h | 89 +++++++++++++++++++++++------ 1 file changed, 70 insertions(+), 19 deletions(-) diff --git a/include/utils/QuantificationUtils.h b/include/utils/QuantificationUtils.h index 8f26c6fa5e..ad72c5a7bd 100644 --- a/include/utils/QuantificationUtils.h +++ b/include/utils/QuantificationUtils.h @@ -5,15 +5,16 @@ #include "FaultTreeUtils.h" #include "MastodonUtils.h" +/********************* Quantification Definition *********************/ class FTAUtils::Quantification { public: enum _dist_t { - PE, - EXP, - GAMMA, + PE, // = 0 + EXP, // = 1 + GAMMA, // = 2 WEIBULL, EXTREME, NORM, @@ -26,8 +27,8 @@ class FTAUtils::Quantification enum _analysis_t { - FRAGILITY, - RISK + FRAGILITY, // = 0 + RISK // = 1 }; Quantification(std::map & params_double, @@ -63,6 +64,7 @@ class FTAUtils::Quantification std::vector vec; vec.clear(); + // Remove all INF from vector for (int index = 0; index < vec_in.size(); index++) { if (!std::isinf(vec_in[index])) @@ -73,6 +75,7 @@ class FTAUtils::Quantification int size = vec.size(); + // ASSERT if (!(size > 0)) { fprintf(stderr, @@ -86,6 +89,8 @@ class FTAUtils::Quantification _pe = vec[0]; + // TODO: Remove 0th element from further calculations + // vec.erase( vec.begin() ); sort(vec.begin(), vec.end()); _mean = accumulate(vec.begin(), vec.end(), 0.0) / size; @@ -114,7 +119,8 @@ class FTAUtils::Quantification } }; - + // TODO: 3 of the following are not supported as they need 1 arg rather than 2 + // for their distribution. Remove these 3 if not needed else add support std::map _str2dist = { {"GAMMA", GAMMA}, {"WEIBULL", WEIBULL}, @@ -129,7 +135,11 @@ class FTAUtils::Quantification vector_double _cut_set_prob; std::map> _b_nodes; - + /** + * Generates probability vector for a specified distribution + * Nomenclature of a and b changes with distribution. + * eg., a => mean, b => std for NORM + */ std::vector getProbVector(_dist_t dist, double a, double b, @@ -139,7 +149,9 @@ class FTAUtils::Quantification _analysis_t analysis, bool uc); - + /** + * Parses and floods probabilties of basic elements + */ vector_string beProb(FTAUtils::Parser parser, int n_sample, int seed, @@ -147,7 +159,10 @@ class FTAUtils::Quantification std::vector intmes, bool uncert); - + /** + * Computes probability for a given cut set on a per set basis based on + * pre-flooded basic elem probs + */ vector_double computeCutSetProb(std::set> cut_sets, int n, bool bypass = false, @@ -161,40 +176,76 @@ class FTAUtils::Quantification std::string bypass_key = "", double bypass_value = 0); - + /** + * Computes accumulated probability for the entire cut set + * 3 Digests are computed: + * 0. Min Max + * 1. Upper Bound + * 2. Top Rare + * + * NOTE + * 1. Its better to compute these 2 together as they have common loops + * which can save time + * 2. Sets the vector cutSetProb + */ std::vector * cutSetProbWDigest(std::set> cut_sets, int n, bool only_min_max = false); - + /** + * Wrapper function for gate probability arith + * AND => a * b + * OR => (1 - a) * (1 - b) + */ double getGateProb(double a, double b, bool is_and); - + /** + * Computes cut-set Fussel-Vesely Important measures + */ vector_double minCutIM(std::vector upperBound); - + /** + * For each basic element: + * 1. Calculate #occurrence in minimal cut set (beCount) + * 2. Store row index of all occurrence (beIndex) + */ std::map beIM(std::set> cut_sets, int n, std::vector upper_bound, std::vector & count_v); - + /** + * Translator function for string -> double + */ vector_double linesToDouble(vector_string lines); - + /** + * Function for getting BIN mean values of intensity measure + */ std::vector getBinMeans(double im_lower, double im_upper, int n_bins); - + /** + * Function for interpolating the hazard curve based on range of intensity + * measures + */ std::vector hazInterp(vector_double hazard, std::vector im_bins); - + /** + * Function for top event fragility + */ std::vector fragility(std::set> cut_sets, int n, std::vector im_bins, double & mu, double & sigma); - + /** + * Function for calculating risk by convoluting hazard and fragility + * Assign risk to basic event probabilites + * + * WARNING + * This consumes _b_nodes and then overwrites it + */ void computeRisk(int n, std::vector hazard); }; -#endif +#endif // _QUANTIFICATION_H From 8a6080dae6883b11063d82de48c4f76f693bebdf Mon Sep 17 00:00:00 2001 From: Hugo Nie Date: Mon, 30 Mar 2020 17:44:12 -0400 Subject: [PATCH 07/10] update faultTreeUtils.C Ref #126 --- src/utils/FaultTreeUtils.C | 60 ++++++++++++++++++++++++++++++++++---- 1 file changed, 54 insertions(+), 6 deletions(-) diff --git a/src/utils/FaultTreeUtils.C b/src/utils/FaultTreeUtils.C index 34a8276814..2f79347af4 100644 --- a/src/utils/FaultTreeUtils.C +++ b/src/utils/FaultTreeUtils.C @@ -8,8 +8,10 @@ FTAUtils::FaultTree::FaultTree(std::string file_name, std::string root) FTAUtils::Parser parser = FTAUtils::Parser(file_name, FTAUtils::Parser::FORMAT_CSV); buildTree(parser); + // Override root node if not default if (root != "") { + // ASSERT if (!(getNode(root))) { fprintf(stderr, @@ -34,16 +36,21 @@ FTAUtils::FaultTree::FaultTree(set_string & sets_link, std::map root_set, sub_sets; - root_set.insert("TOP"); + root_set.insert("TOP"); // "TOP", size: 1 _sets.insert(root_set); + // Start with root node to expand on heirarchy cutSetsExpand(getNode("TOP")); - + // Check if first row is empty + // NOTE: Since its std::set, only first row could + // be empty if (_sets.begin()->size() == 0) { _sets.erase(_sets.begin()); @@ -64,17 +71,22 @@ FTAUtils::FaultTree::buildTree(FTAUtils::Parser parser) { line = parser.yieldLine(); + // Stop if no new line to process if (line.size() == 0) break; + // Stash name, operator _node * node = new _node(line[0], str2Operator(line[1])); + // Add children for (int i = 2; i < line.size(); i++) node->_child.push_back(line[i]); + // Stash the first entry as ROOT of tree if (_node_d_b.size() == 0) _root = line[0]; + // Add the newly created node to node lookup hashmap _node_d_b[line[0]] = node; } @@ -85,6 +97,7 @@ FTAUtils::FaultTree::_operator_t FTAUtils::FaultTree::str2Operator(std::string op) { std::string op_s = FTAUtils::str2Upper(op, true); + // ASSERT if (_opDict.count(op_s) == 0) { fprintf(stderr, @@ -101,16 +114,21 @@ FTAUtils::FaultTree::str2Operator(std::string op) set_string FTAUtils::FaultTree::computeMinimumCutSets() { + // Clearing up sets vector for safety rmSets(); + // Add root to set std::set root_set, sub_sets; root_set.insert(getRoot()); _sets.insert(root_set); + // Start with root node to expand on heirarchy cutSetsExpand(getNode(getRoot())); - + // Check if first row is empty + // NOTE: Since its std::set, only first row could + // be empty if (_sets.begin()->size() == 0) { _sets.erase(_sets.begin()); @@ -136,8 +154,10 @@ FTAUtils::FaultTree::removeSubsets() includes(row2_it->begin(), row2_it->end(), row1_it->begin(), row1_it->end()); bool row2_is_subset = includes(row1_it->begin(), row1_it->end(), row2_it->begin(), row2_it->end()); + // Remove row1 if row2 is a subset of row1 by marking it if (row2_is_subset) rm_its.insert(index1); + // Remove row2 if row1 is a subset of row1 by marking it else if (row1_is_subset) rm_its.insert(index2); index2++; @@ -146,6 +166,7 @@ FTAUtils::FaultTree::removeSubsets() } uint64_t index = 0; + // Remove rows marked to be removed for (std::set::iterator it = rm_its.begin(); it != rm_its.end(); ++it) { _sets.erase(next(_sets.begin(), *it - index)); @@ -156,6 +177,7 @@ FTAUtils::FaultTree::removeSubsets() void FTAUtils::FaultTree::cutSetsExpand(_node * node) { + // ASSERT if (!node) { fprintf(stderr, @@ -166,6 +188,7 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) abort(); } + // ASSERT if (node->_child.size() == 0) { fprintf(stderr, @@ -177,21 +200,28 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) abort(); } - + // New rows which have to be appended at end. Adding them might + // result in data hazards set_string mk_rows; + // 1. Iterate through entire list and match for own node's name for (set_string::iterator row_it = _sets.begin(); row_it != _sets.end(); ++row_it) { std::set row(*row_it); + // Search for the element in row std::set::iterator it = find(row.begin(), row.end(), node->_name); if (it != row.end()) { + // Erase self to add children row.erase(it); - + // 2. Replace self with children based on operation + // (i) . Replace self with children in same row if AND + // (ii). Replace self with child one per row if OR switch (node->_op) { case AND: + // Append all children in same row for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) { row.insert(node->_child[c_id]); @@ -199,6 +229,7 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) break; case OR: + // Replicate line and append one child per row for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) { std::set row_set(row); @@ -206,11 +237,13 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) mk_rows.insert(row_set); } + // Clear this row to be removed at end (postprocessing) row.clear(); break; default: { + // ASSERT fprintf(stderr, "[ASSERT] In File: %s, Line: %d => " "Unknown Operator found!.\n", @@ -222,17 +255,22 @@ FTAUtils::FaultTree::cutSetsExpand(_node * node) } } - + // The following 2 lines might result in an empty row to be pushed. + // If we gate insertion with count, we might have a data hazard due + // to wrong increment of iterators. + // A quick fix to this is to check for first row and delete it at end _sets.erase(row_it); _sets.insert(row); } + // Add newly created rows for (set_string::iterator it = mk_rows.begin(); it != mk_rows.end(); ++it) { std::set row(*it); _sets.insert(*it); } + // 3. Recurse on each non leaf child for (uint64_t c_id = 0; c_id < node->_child.size(); c_id++) { _node * child = getNode(node->_child[c_id]); @@ -270,6 +308,7 @@ FTAUtils::FaultTree::getCutSets() FTAUtils::Parser::Parser(std::string fileName, FTAUtils::Parser::parseFormatT format) { + // Assertion on supported parsing formats if (format != FORMAT_CSV) { fprintf(stderr, @@ -298,6 +337,7 @@ FTAUtils::Parser::yieldLines() { line = FTAUtils::Parser::yieldLine(); + // Stop if no new line to process if (line.size() == 0) break; lines.push_back(line); @@ -308,6 +348,7 @@ FTAUtils::Parser::yieldLines() std::vector FTAUtils::Parser::yieldLine() { + // ASSERT if (!(fileP->is_open())) { fprintf(stderr, @@ -321,6 +362,7 @@ FTAUtils::Parser::yieldLine() std::string buffer; std::vector line; + // Get a line that has something except \n if (getline(*fileP, buffer)) { std::string token; @@ -364,8 +406,10 @@ FTAUtils::interpolate(vector_double data, double x, bool extrapolate) { int size = data.size(); + // find left end of interval for interpolation int i = 0; + // special case: beyond right end if (x >= data[size - 2][0]) { i = size - 2; @@ -376,8 +420,10 @@ FTAUtils::interpolate(vector_double data, double x, bool extrapolate) i++; } + // points on either side (unless beyond ends) double xL = data[i][0], yL = data[i][1], xR = data[i + 1][0], yR = data[i + 1][1]; + // if beyond ends of array and not extrapolating if (!extrapolate) { if (x < xL) @@ -386,8 +432,10 @@ FTAUtils::interpolate(vector_double data, double x, bool extrapolate) yL = yR; } + // gradient double dydx = (yR - yL) / (xR - xL); + // linear interpolation return yL + dydx * (x - xL); } From 7adf6186ec88b5fad5da66903223511265dec1a5 Mon Sep 17 00:00:00 2001 From: Hugo Nie Date: Mon, 30 Mar 2020 17:48:19 -0400 Subject: [PATCH 08/10] update QuantificationUtils.C --- src/utils/QuantificationUtils.C | 52 +++++++++++++++++++++++++++++---- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/src/utils/QuantificationUtils.C b/src/utils/QuantificationUtils.C index 4c578b0cb8..e740a5d4fe 100644 --- a/src/utils/QuantificationUtils.C +++ b/src/utils/QuantificationUtils.C @@ -17,6 +17,7 @@ FTAUtils::Quantification::Quantification(std::map & int n_sample, int seed) { + //---------------- SAVE PARAMETERS ------------------------------------------ std::vector params_double1; params_double1.push_back(im_lower); params_double1.push_back(im_upper); @@ -30,44 +31,55 @@ FTAUtils::Quantification::Quantification(std::map & params_bool.insert(std::pair("uncertainty", uncertainty)); params_analysis_t.insert(std::pair("analysis", analysis)); - + //---------------- COMPUTE -------------------------------------------------- + // Construct the fault tree with MOCUS algo for minimal cut sets FTAUtils::FaultTree ft = FTAUtils::FaultTree(events_file, root); Parser parser_events = Parser(events_file, Parser::FORMAT_CSV); vector_string lines_events = parser_events.yieldLines(); set_string cut_sets = ft.getCutSets(); params_string.insert(std::pair("events_files", lines_events)); + // Read basic events file FTAUtils::Parser parser_event_prob = FTAUtils::Parser(events_prob_file, FTAUtils::Parser::FORMAT_CSV); vector_string lines_events_prob; if (analysis == FRAGILITY) { + // Read and interpolate hazard curve FTAUtils::Parser parser_hazard = FTAUtils::Parser(hazard_file, FTAUtils::Parser::FORMAT_CSV); vector_string lines_hazard = parser_hazard.yieldLines(); vector_double hazard = FTAUtils::Quantification::linesToDouble(lines_hazard); params_double.insert(pair_string_double("hazard", hazard)); + // List of Intensity Measure bin values std::vector im_bins = getBinMeans(im_lower, im_upper, n_bins); + // List of Intensity Measure bin values std::vector hazard_freq = hazInterp(hazard, im_bins); + // Dictionary of basic events for fragility input lines_events_prob = beProb(parser_event_prob, n_sample, seed, analysis, im_bins, uncertainty); + // Top event fragility (lognormal parameters) double mu, sigma; fragility(cut_sets, n_bins, im_bins, mu, sigma); + // Dictionary of basic events risk (convoluting fragility and hazard) computeRisk(n_bins, hazard_freq); } else { + // Dictionary of basic events risk (risk inputs) lines_events_prob = beProb(parser_event_prob, n_sample, seed, analysis, std::vector(), uncertainty); } params_string.insert(std::pair("basic_events", lines_events_prob)); - + // Calculate minimal cut sets probability + // Digest cut set probability to unified probability + // Adding 1 to accomodate mean as 0th element std::vector * digest = cutSetProbWDigest(cut_sets, n_sample + 1); vector_double mc_im = minCutIM(digest[1]); Stats s0(digest[0]); @@ -76,6 +88,7 @@ FTAUtils::Quantification::Quantification(std::map & std::vector count_v; std::map be_stats = beIM(cut_sets, n_sample + 1, digest[1], count_v); + // return the FTA top event risk std::vector fta_0; vector_double fta_1; fta_0.push_back(s0._pe); @@ -84,10 +97,12 @@ FTAUtils::Quantification::Quantification(std::map & fta_1.push_back(fta_0); params_double.insert(std::make_pair("fta", fta_1)); + // Return the Risk Reduction Ratio for B1, B2, B3, B4, B5 params_double.insert(pair_string_double("fv", be_stats["fv"])); params_double.insert(pair_string_double("rrr", be_stats["rrr"])); params_double.insert(pair_string_double("rir", be_stats["rir"])); + // Return the Risk Reduction Difference for B1, B2, B3, B4, B5 params_double.insert(pair_string_double("rri", be_stats["rri"])); params_double.insert(pair_string_double("rii", be_stats["rii"])); params_double.insert(pair_string_double("bi", be_stats["bi"])); @@ -97,8 +112,10 @@ std::vector FTAUtils::Quantification::hazInterp(vector_double hazard, std::vector im_bins) { std::vector data; + // Convert all hazard values to log10 for (int index = 0; index < hazard.size(); index++) { + // ASSERT if (hazard[index].size() != 2) { fprintf(stderr, @@ -113,6 +130,7 @@ FTAUtils::Quantification::hazInterp(vector_double hazard, std::vector im hazard[index][1] = log10(hazard[index][1]); } + // Interpolate for (int index = 0; index < im_bins.size(); index++) { data.push_back(pow(10, interpolate(hazard, log10(im_bins[index]), true))); @@ -160,6 +178,7 @@ FTAUtils::genQuantificationRVec(double dataPoint, int n, std::vector rv) { for (int index = 0; index < n; index++) { + // double dataPoint = double(gen); rv.push_back(FTAUtils::Clip(dataPoint, 0.0, 1.0)); } return rv; @@ -181,6 +200,7 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, if (analysis == FRAGILITY) { + // ASSERT if (dist != LNORM) { fprintf(stderr, @@ -214,6 +234,7 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, break; default: { + // ASSERT fprintf(stderr, "[ASSERT] In File: %s, Line: %d => " "Un-supported dist found.\n", @@ -244,6 +265,7 @@ FTAUtils::Quantification::getProbVector(_dist_t dist, break; default: { + // ASSERT fprintf(stderr, "[ASSERT] In File: %s, Line: %d => " "Un-supported dist found.\n", @@ -271,11 +293,13 @@ FTAUtils::Quantification::beProb(FTAUtils::Parser parser, { line = parser.yieldLine(); + // Stop if no new line to process if (line.size() != 0) lines.push_back(line); else break; + // Stash name, probability vector double b = line.size() > 3 ? stod(line[3]) : 0; _b_nodes[line[0]] = getProbVector( _str2dist[line[1]], stod(line[2]), b, n_sample, seed, intmes, analysis, uncert); @@ -289,14 +313,17 @@ FTAUtils::Quantification::cutSetProbWDigest(set_string cut_sets, int n, bool onl { std::vector * digest = new std::vector[3]; + // Digest 2: Generate power set to generate all possible combinations vector_double ps_p; int pow_set_size = pow(2, cut_sets.size()); + // Start from 1 as NULL set is not needed for (int index = 1; index < pow_set_size; index++) { std::set power_set_row; for (int j = 0; j < cut_sets.size(); j++) { - + // Check if jth bit in the index is set + // If set then form row of power set if (index & (1 << j)) { std::set j_elem = *next(cut_sets.begin(), j); @@ -306,12 +333,15 @@ FTAUtils::Quantification::cutSetProbWDigest(set_string cut_sets, int n, bool onl } } } + // For elements taken n at a time, if n is odd, positive else negative bool is_positive = (__builtin_popcount(index) % 2); std::vector cut_prob = cutSetRowProb(power_set_row, n, is_positive); ps_p.push_back(cut_prob); } - + // Digest 0, 1, 2 + // Avoiding an extra function call to interm event like in the python + // counterpart for performance reasons (loop merging) _cut_set_prob = computeCutSetProb(cut_sets, n); for (int index = 0; index < n; index++) { @@ -342,7 +372,9 @@ vector_double FTAUtils::Quantification::computeCutSetProb( set_string cut_sets, int n, bool bypass, std::string bypass_key, double bypass_value) { - + // For each row in cut set, compute: + // 1. For each sample in the generated vector, compute AND gate + // probability analysis for each basic element vector_double quant; for (set_string::iterator row = cut_sets.begin(); row != cut_sets.end(); ++row) { @@ -366,6 +398,7 @@ FTAUtils::Quantification::cutSetRowProb(std::set row, double value = 1; for (std::set::iterator col = row.begin(); col != row.end(); ++col) { + // ASSERT if (_b_nodes.find(*col) == _b_nodes.end()) { fprintf(stderr, @@ -427,9 +460,14 @@ std::vector FTAUtils::Quantification::fragility( set_string cut_sets, int n, std::vector im_bins, double & mu, double & sigma) { + // 1. Calculate TOP event fragility using min-max approach (exact) {digest[0]} std::vector * digest = cutSetProbWDigest(cut_sets, n, true); + // 2. lognormal parameters of TOP Event fragility + // solveLnParams(im_bins, digest[0], mu, sigma); + // Use the maximizelikehood function instead of solver function + // Use SGD for optimization std::vector loc_space = {0.01, 2}; std::vector sca_space = {0.01, 1}; std::vector max_values = MastodonUtils::maximizeLogLikelihood( @@ -437,6 +475,7 @@ FTAUtils::Quantification::fragility( mu = max_values[0]; sigma = max_values[1]; + // TODO: 3. TOP event Fragility plot return digest[0]; } @@ -451,6 +490,7 @@ FTAUtils::Quantification::beIM(set_string cut_sets, bn_it != _b_nodes.end(); ++bn_it) { + // Generate available vector to save computation on per index loop std::vector available; int count = 0; for (set_string::iterator cs_it = cut_sets.begin(); cs_it != cut_sets.end(); ++cs_it) @@ -467,6 +507,7 @@ FTAUtils::Quantification::beIM(set_string cut_sets, vector_double mc_p1 = computeCutSetProb(cut_sets, n, true, bn_it->first, 1); vector_double mc_p0 = computeCutSetProb(cut_sets, n, true, bn_it->first, 0); + // OR it with loop merging for (int index = 0; index < n; index++) { double f0_val = 1; @@ -490,6 +531,7 @@ FTAUtils::Quantification::beIM(set_string cut_sets, rii.push_back(f1_val - upper_bound[index]); bi.push_back(f1_val - f0_val); } + // Calculate stats for current row and stash Stats fv_stats(fv); Stats rrr_stats(rrr); Stats rri_stats(rri); From 128b3d0d8d5f8651175be30c2bad259c2f6c1aff Mon Sep 17 00:00:00 2001 From: Hugo Nie Date: Mon, 30 Mar 2020 17:50:27 -0400 Subject: [PATCH 09/10] update testFaultTreeUtils.C --- unit/src/TestFaultTreeUtils.C | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/unit/src/TestFaultTreeUtils.C b/unit/src/TestFaultTreeUtils.C index 9db7f291e5..ccfc21d31a 100644 --- a/unit/src/TestFaultTreeUtils.C +++ b/unit/src/TestFaultTreeUtils.C @@ -1,58 +1,71 @@ +// MOOSE includes #include "gtest/gtest.h" #include "MooseUtils.h" #include "Conversion.h" +// Boost distribution includes #include "BoostDistribution.h" +// MASTODON includes #include "MastodonUtils.h" +// QUANTIFICATIONUTILS includes #include "FaultTreeUtils.h" typedef FTAUtils::FaultTree::_node ft_node_test; +// Helper to make sure the test tree is correct. void assertTree(std::map tree_node) { + // TOP EXPECT_EQ(tree_node["TOP"]->_name, "TOP"); EXPECT_EQ(tree_node["TOP"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["TOP"]->_child.at(0), "GATE1"); EXPECT_EQ(tree_node["TOP"]->_child.at(1), "GATE2"); + // GATE1 EXPECT_EQ(tree_node["GATE1"]->_name, "GATE1"); EXPECT_EQ(tree_node["GATE1"]->_op, FTAUtils::FaultTree::OR); EXPECT_EQ(tree_node["GATE1"]->_child.at(0), "FT-N/m-1"); EXPECT_EQ(tree_node["GATE1"]->_child.at(1), "FT-N/m-2"); EXPECT_EQ(tree_node["GATE1"]->_child.at(2), "FT-N/m-3"); + // GATE2 EXPECT_EQ(tree_node["GATE2"]->_name, "GATE2"); EXPECT_EQ(tree_node["GATE2"]->_op, FTAUtils::FaultTree::OR); EXPECT_EQ(tree_node["GATE2"]->_child.at(0), "B1"); EXPECT_EQ(tree_node["GATE2"]->_child.at(1), "B3"); EXPECT_EQ(tree_node["GATE2"]->_child.at(2), "B4"); + // FT-N/m-1 EXPECT_EQ(tree_node["FT-N/m-1"]->_name, "FT-N/m-1"); EXPECT_EQ(tree_node["FT-N/m-1"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(0), "GATE3"); EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(1), "B3"); EXPECT_EQ(tree_node["FT-N/m-1"]->_child.at(2), "B5"); + // FT-N/m-2 EXPECT_EQ(tree_node["FT-N/m-2"]->_name, "FT-N/m-2"); EXPECT_EQ(tree_node["FT-N/m-2"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["FT-N/m-2"]->_child.at(0), "GATE3"); EXPECT_EQ(tree_node["FT-N/m-2"]->_child.at(1), "B1"); + // FT-N/m-3 EXPECT_EQ(tree_node["FT-N/m-3"]->_name, "FT-N/m-3"); EXPECT_EQ(tree_node["FT-N/m-3"]->_op, FTAUtils::FaultTree::AND); EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(0), "B3"); EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(1), "B5"); EXPECT_EQ(tree_node["FT-N/m-3"]->_child.at(2), "B1"); + // GATE3 EXPECT_EQ(tree_node["GATE3"]->_name, "GATE3"); EXPECT_EQ(tree_node["GATE3"]->_op, FTAUtils::FaultTree::OR); EXPECT_EQ(tree_node["GATE3"]->_child.at(0), "B2"); EXPECT_EQ(tree_node["GATE3"]->_child.at(1), "B4"); } +// Helper to make sure the mocus is correct. void assertMocus(set_string tree_mocus) { @@ -79,6 +92,7 @@ assertMocus(set_string tree_mocus) } } +// Create a test event list. std::map eventList() { @@ -104,19 +118,26 @@ eventList() for (int i = 0; i < line.size(); i++) { ft_node_test * node = new ft_node_test(line[i][0], op[i]); + // Add children for (int j = 2; j < line[i].size(); j++) node->_child.push_back(line[i][j]); + // Add the newly created node to node lookup hashmap _node_d_b[line[i][0]] = node; } return _node_d_b; } +// For input std::vector file_lists; TEST(FTAUtils, FaultTree) { + // ==================== TestCase for FaultTree Object =================== + // ---------------- Test error message for events input ---------------- + + // ###### File Inputs ###### try { @@ -126,10 +147,14 @@ TEST(FTAUtils, FaultTree) } catch (FTAUtils::CException e) { - + /* + * ------- IO Error ------- + * does not exist + */ EXPECT_EQ(e.msg, "Unable to open file."); } + // ---------------- Test creating tree from list ---------------- std::map tree_node_from_list = eventList(); set_string tree_mocus_from_list; @@ -137,6 +162,7 @@ TEST(FTAUtils, FaultTree) assertTree(tree_node_from_list); assertMocus(tree_mocus_from_list); + // ---------------- Test creating tree from file ---------------- file_lists = {"logic2.txt", ""}; From a94027a0e5c139c0db72f82b9bfb49ebe4ca3b6e Mon Sep 17 00:00:00 2001 From: Hugo Nie Date: Mon, 30 Mar 2020 17:51:47 -0400 Subject: [PATCH 10/10] update QuantificationUtils.C Ref #126 --- unit/src/TestQuantificationUtils.C | 80 ++++++++++++++++++++++++++++-- 1 file changed, 75 insertions(+), 5 deletions(-) diff --git a/unit/src/TestQuantificationUtils.C b/unit/src/TestQuantificationUtils.C index 9920d089f1..2295e91c84 100644 --- a/unit/src/TestQuantificationUtils.C +++ b/unit/src/TestQuantificationUtils.C @@ -1,13 +1,18 @@ +// MOOSE includes #include "gtest/gtest.h" #include "MooseUtils.h" #include "Conversion.h" +// Boost distribution includes #include "BoostDistribution.h" +// MASTODON includes #include "MastodonUtils.h" +// QUANTIFICATIONUTILS includes #include "QuantificationUtils.h" +// For input std::vector file_lists_q; std::vector IM; int nbins; @@ -15,6 +20,7 @@ bool uncertainty; int nsamp; int seed; +// For output std::map>> params_double_1, params_double_2, params_double_3, params_double_4, params_double_5, params_double_6, params_double_7, params_double_8, params_double_9; @@ -33,9 +39,18 @@ std::map params_analysis_t_1 std::map params_int_1, params_int_2, params_int_3, params_int_4, params_int_5, params_int_6, params_int_7, params_int_8, params_int_9; +// TestCase for Quantification Object TEST(FTAUtils, Quantification) { + // =========================================================================== + // =============================== Test Inputs =============================== + // =========================================================================== + + // +++++++++++++++++++++++ Input values +++++++++++++++++++++++ + + // ############### File Inputs ############### + try { file_lists_q = { @@ -58,7 +73,10 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - + /* + * ------- IO Error ------- + * does not exist + */ EXPECT_EQ(e.msg, "Unable to open file."); } @@ -83,10 +101,14 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - + /* + * ------- Type Error ------- + * must be a filename or a list + */ EXPECT_EQ(e.msg, "Unable to open file."); } + // ############## Inputs for Fragility ############## file_lists_q = {"logic1.txt", "logic1_bas_events_LNORM.txt", "hazard.txt"}; IM = {0.1, 4}; @@ -104,6 +126,7 @@ TEST(FTAUtils, Quantification) IM[1], nbins); + // >>>>>>>> logic Value Check <<<<<<<< std::vector> event_files_matrix3{{"TE", "OR", "IE3", "IE4"}, {"IE4", "OR", "C4"}, @@ -112,6 +135,7 @@ TEST(FTAUtils, Quantification) {"IE1", "OR", "C1"}}; EXPECT_EQ(params_string_3["events_files"], event_files_matrix3); + // >>>>>>>> Basic Events Value Check <<<<<<<< std::vector> basic_events_matrix3{{"C1", "LNORM", "1.88", "0.5"}, {"C2", "LNORM", "3.78", "0.79"}, @@ -120,9 +144,11 @@ TEST(FTAUtils, Quantification) std::vector> basic_events_3 = params_string_3["basic_events"]; EXPECT_EQ(basic_events_3, basic_events_matrix3); + // >>>>>>>> antype Value Check <<<<<<<< EXPECT_EQ(params_analysis_t_3["analysis"], FTAUtils::Quantification::FRAGILITY); + // >>>>>>>> hazard Value Check <<<<<<<< std::vector> matrix_hazard{{0.0608, 0.01}, {0.2124, 0.001}, @@ -133,13 +159,16 @@ TEST(FTAUtils, Quantification) std::vector> hazard = params_double_3["hazard"]; EXPECT_EQ(hazard, matrix_hazard); + // >>>>>>>> imrange Value Check <<<<<<<< EXPECT_EQ(params_double_3["IM"][0][0], 0.1); EXPECT_EQ(params_double_3["IM"][0][1], 4); + // >>>>>>>> nbins Value Check <<<<<<<< EXPECT_EQ(params_int_3["n_bins"], 15); + // >>>>>>>> nbins Value Check <<<<<<<< try { @@ -165,10 +194,14 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - + /* + * ------- Value Error ------- + * The supplied value of nbins must be a +ve integer + */ EXPECT_EQ(e.msg, "ValueError"); } + // ############## Inputs for Risk analysis (not fragility) ############## file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt"}; @@ -181,6 +214,7 @@ TEST(FTAUtils, Quantification) file_lists_q[1], FTAUtils::Quantification::RISK); + // >>>>>>>> logic Value Check <<<<<<<< std::vector> event_files_matrix5{ {"TOP", "AND", "GATE1", "GATE2"}, @@ -193,6 +227,7 @@ TEST(FTAUtils, Quantification) EXPECT_EQ(params_string_5["events_files"], event_files_matrix5); + // >>>>>>>> Basic Events Value Check <<<<<<<< std::vector> basic_events_matrix5{{"B1", "PE", "0.01"}, {"B2", "PE", "0.02"}, @@ -204,18 +239,23 @@ TEST(FTAUtils, Quantification) EXPECT_EQ(basic_events_5, basic_events_matrix5); + // >>>>>>>> antype Value Check <<<<<<<< EXPECT_NE(params_analysis_t_5["analysis"], FTAUtils::Quantification::FRAGILITY); + // >>>>>>>> uncertainty Value Check <<<<<<<< EXPECT_EQ(params_bool_5["uncertainty"], false); + // >>>>>>>> nsamp Value Check <<<<<<<< EXPECT_EQ(params_int_5["nsamp"], 1); + // >>>>>>>> seed Value Check <<<<<<<< EXPECT_EQ(params_int_5["seed"], 0); + // ###### Testing for input errors making sure parameters are input correctly ###### file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; @@ -241,16 +281,21 @@ TEST(FTAUtils, Quantification) nsamp, seed); + // >>>>>>>> uncertainty Value Check <<<<<<<< EXPECT_EQ(params_bool_6["uncertainty"], true); + // >>>>>>>> nsamp Value Check <<<<<<<< EXPECT_EQ(params_int_6["nsamp"], 1000); + // >>>>>>>> seed Value Check <<<<<<<< EXPECT_EQ(params_int_6["seed"], 436546754); + // ############## Testing for input type errors and value errors ############## + // >>>>>>>> nsamp Value Check <<<<<<<< try { @@ -284,10 +329,14 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - + /* + * ------- Value Error ------- + * The supplied value of nsamp must be a +ve integer. + */ EXPECT_EQ(e.msg, "ValueError"); } + // >>>>>>>> seed Value Check <<<<<<<< try { @@ -320,11 +369,20 @@ TEST(FTAUtils, Quantification) } catch (FTAUtils::CException e) { - + /* + * ------- Value Error ------- + * The supplied value of seed must be a +ve integer. + */ EXPECT_EQ(e.msg, "ValueError"); } + // =========================================================================== + // ============================== Test TOP Risk ============================== + // =========================================================================== + + // ############ Function for asserting FTA top event risk ############ + // ++++++++++ Check FTA top event risk ++++++++++ file_lists_q = {"logic2.txt", "logic2_bas_events_PE.txt", "hazard.txt"}; IM = {0.1, 4}; @@ -342,45 +400,57 @@ TEST(FTAUtils, Quantification) IM[1], nbins); + // >>>>>>>> Risk Value Check <<<<<<<< + // min-max EXPECT_EQ(std::to_string(params_double_9["fta"][0][0]), "0.000694"); + // upper bound EXPECT_EQ(std::to_string(params_double_9["fta"][0][1]), "0.000705"); + // rare event EXPECT_EQ(std::to_string(params_double_9["fta"][0][2]), "0.000705"); + // >>>>>>>> IMratio Value Check <<<<<<<< + // Fussell-Vesely for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["fv"][0][0]), "0.872395"); EXPECT_EQ(std::to_string(params_double_9["fv"][1][0]), "0.326300"); EXPECT_EQ(std::to_string(params_double_9["fv"][2][0]), "0.148963"); EXPECT_EQ(std::to_string(params_double_9["fv"][3][0]), "0.652584"); EXPECT_EQ(std::to_string(params_double_9["fv"][4][0]), "0.148963"); + // Risk Reduction Ratio for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rrr"][0][0]), "7.831866"); EXPECT_EQ(std::to_string(params_double_9["rrr"][1][0]), "1.483999"); EXPECT_EQ(std::to_string(params_double_9["rrr"][2][0]), "1.174913"); EXPECT_EQ(std::to_string(params_double_9["rrr"][3][0]), "2.877066"); EXPECT_EQ(std::to_string(params_double_9["rrr"][4][0]), "1.174913"); + // Risk Increase Ratio for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rir"][0][0]), "86.111103"); EXPECT_EQ(std::to_string(params_double_9["rir"][1][0]), "16.960273"); EXPECT_EQ(std::to_string(params_double_9["rir"][2][0]), "5.808755"); EXPECT_EQ(std::to_string(params_double_9["rir"][3][0]), "16.637742"); EXPECT_EQ(std::to_string(params_double_9["rir"][4][0]), "3.826894"); + // >>>>>>>> IMdiff Value Check <<<<<<<< + // Risk Reduction Difference for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rri"][0][0]), "0.000615"); EXPECT_EQ(std::to_string(params_double_9["rri"][1][0]), "0.000230"); EXPECT_EQ(std::to_string(params_double_9["rri"][2][0]), "0.000105"); EXPECT_EQ(std::to_string(params_double_9["rri"][3][0]), "0.000460"); EXPECT_EQ(std::to_string(params_double_9["rri"][4][0]), "0.000105"); + // Risk Increase Difference for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["rii"][0][0]), "0.059991"); EXPECT_EQ(std::to_string(params_double_9["rii"][1][0]), "0.011250"); EXPECT_EQ(std::to_string(params_double_9["rii"][2][0]), "0.003389"); EXPECT_EQ(std::to_string(params_double_9["rii"][3][0]), "0.011022"); EXPECT_EQ(std::to_string(params_double_9["rii"][4][0]), "0.001993"); + // Birnbaum Difference for B1, B2, B3, B4, B5 EXPECT_EQ(std::to_string(params_double_9["bi"][0][0]), "0.060606"); EXPECT_EQ(std::to_string(params_double_9["bi"][1][0]), "0.011480"); EXPECT_EQ(std::to_string(params_double_9["bi"][2][0]), "0.003494");