diff --git a/include/utils/FaultTreeUtils.h b/include/utils/FaultTreeUtils.h new file mode 100644 index 0000000000..57df9b3a81 --- /dev/null +++ b/include/utils/FaultTreeUtils.h @@ -0,0 +1,192 @@ +#ifndef _FAULT_TREE_H +#define _FAULT_TREE_H + +#include +#include +#include +#include + +typedef std::vector> vector_double; +typedef std::vector> vector_string; +typedef std::set> set_string; +typedef std::pair pair_string_double; + +namespace FTAUtils +{ +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); + +/** + * 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 +{ +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; + } + }; + + /** + * 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; +}; + +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..ad72c5a7bd --- /dev/null +++ b/include/utils/QuantificationUtils.h @@ -0,0 +1,251 @@ +#ifndef _QUANTIFICATION_H +#define _QUANTIFICATION_H + +#include +#include "FaultTreeUtils.h" +#include "MastodonUtils.h" + +/********************* Quantification Definition *********************/ + +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 + 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]; + + // 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])); + } + }; + + // 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}, + {"EXTREME", EXTREME}, + {"NORM", NORM}, + {"LNORM", LNORM}, + {"PE", PE}, + {"CAUCHY", CAUCHY}, + {"FISHER_F", FISHER_F}, + }; + + 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, + int n, + int seed, + std::vector im, + _analysis_t analysis, + bool uc); + + /** + * 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); + + /** + * 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 diff --git a/src/utils/FaultTreeUtils.C b/src/utils/FaultTreeUtils.C new file mode 100644 index 0000000000..2f79347af4 --- /dev/null +++ b/src/utils/FaultTreeUtils.C @@ -0,0 +1,456 @@ +#include "FaultTreeUtils.h" +#include "MastodonUtils.h" + +typedef FTAUtils::FaultTree::_node ft_node; + +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, + "[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(); +} + +FTAUtils::FaultTree::FaultTree(set_string & sets_link, std::map _node_base) +{ + 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; +} + +FTAUtils::FaultTree::~FaultTree() {} + +std::map +FTAUtils::FaultTree::buildTree(FTAUtils::Parser parser) +{ + 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; +} + +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, + "[ASSERT] In File: %s, Line: %d => " + "Illegal Operator found: %s.\n", + __FILE__, + __LINE__, + op.c_str()); + abort(); + } + return _opDict[op_s]; +} + +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()); + } + + removeSubsets(); + + return _sets; +} + +void +FTAUtils::FaultTree::removeSubsets() +{ + std::set rm_its; + + uint64_t index1 = 0; + for (set_string::iterator row1_it = _sets.begin(); row1_it != _sets.end(); ++row1_it) + { + uint64_t index2 = index1 + 1; + 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()); + 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++; + } +} + +void +FTAUtils::FaultTree::cutSetsExpand(_node * node) +{ + // 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 + 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]); + } + 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 + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Unknown Operator found!.\n", + __FILE__, + __LINE__); + abort(); + } + 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 (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]); + if (child) + cutSetsExpand(child); + } +} + +void +FTAUtils::FaultTree::rmSets() +{ + for (set_string::iterator row_it = _sets.begin(); row_it != _sets.end(); ++row_it) + { + _sets.erase(row_it); + } +} + +ft_node * +FTAUtils::FaultTree::getNode(std::string name) +{ + return (_node_d_b.count(name) != 0) ? _node_d_b[name] : NULL; +} + +std::string +FTAUtils::FaultTree::getRoot() +{ + return _root; +} + +set_string +FTAUtils::FaultTree::getCutSets() +{ + return _sets; +} + +FTAUtils::Parser::Parser(std::string fileName, FTAUtils::Parser::parseFormatT format) +{ + // Assertion on supported parsing formats + 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); + + if (!fileP->is_open()) + throw FTAUtils::CException("Unable to open file."); +} + +FTAUtils::Parser::~Parser() {} + +vector_string +FTAUtils::Parser::yieldLines() +{ + vector_string 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; +} + +std::vector +FTAUtils::Parser::yieldLine() +{ + // 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; + + // 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; +} + +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)); +} + +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; +} + +double +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; + } + 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); +} + +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..e740a5d4fe --- /dev/null +++ b/src/utils/QuantificationUtils.C @@ -0,0 +1,580 @@ +#include "QuantificationUtils.h" + +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); + + vector_double params_double2; + params_double2.push_back(params_double1); + 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)); + 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]); + 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); + + // return the FTA top event risk + std::vector fta_0; + 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)); + + // 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"])); +} + +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, + "[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]); + } + + // Interpolate + for (int index = 0; index < im_bins.size(); index++) + { + data.push_back(pow(10, interpolate(hazard, log10(im_bins[index]), true))); + } + return data; +} + +std::vector +FTAUtils::Quantification::getBinMeans(double im_lower, double im_upper, int n_bins) +{ + 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; +} + +vector_double +FTAUtils::Quantification::linesToDouble(vector_string lines) +{ + vector_double 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; +} + +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, + double b, + int n, + int seed, + std::vector im, + _analysis_t analysis, + bool uc) +{ + std::vector rv; + + std::default_random_engine gen(seed); + + if (analysis == FRAGILITY) + { + // 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)); + } + } + else + { + if (!uc) + { + switch (dist) + { + case PE: + { + rv = FTAUtils::genQuantificationRVec(double(a), 2, rv); + } + break; + case NORM: + { + rv = FTAUtils::genQuantificationRVec(double(a), 2, rv); + } + break; + default: + { + // ASSERT + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Un-supported dist found.\n", + __FILE__, + __LINE__); + abort(); + } + } + } + else + { + switch (dist) + { + case PE: + { + rv = FTAUtils::genQuantificationRVec(double(a), n, rv); + } + break; + case NORM: + { + double rv_norm = FTAUtils::Clip(a, 0, 1); + rv.push_back(rv_norm); + + std::normal_distribution d(a, b); + + rv = FTAUtils::genQuantificationRVec(d(gen), n, rv); + } + break; + default: + { + // ASSERT + fprintf(stderr, + "[ASSERT] In File: %s, Line: %d => " + "Un-supported dist found.\n", + __FILE__, + __LINE__); + abort(); + } + } + } + } + return rv; +} + +vector_string +FTAUtils::Quantification::beProb(FTAUtils::Parser parser, + int n_sample, + int seed, + _analysis_t analysis, + std::vector intmes, + bool uncert) +{ + std::vector line; + vector_string 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; +} + +std::vector * +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 + 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); + 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 (vector_double::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 (vector_double::iterator it = ps_p.begin(); it != ps_p.end(); ++it) + { + min_max += (*it)[index]; + } + digest[0].push_back(min_max); + } + + return digest; +} + +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) + { + quant.push_back(cutSetRowProb(*row, n, true, bypass, bypass_key, bypass_value)); + } + return quant; +} + +std::vector +FTAUtils::Quantification::cutSetRowProb(std::set row, + int n, + bool sign_positive, + bool bypass, + std::string bypass_key, + double bypass_value) +{ + 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 + 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); + } + prob.push_back(value * sign_m); + } + return prob; +} + +double +FTAUtils::Quantification::getGateProb(double a, double b, bool is_and) +{ + return is_and ? a * b : (1.0 - a) * (1.0 - b); +} + +vector_double +FTAUtils::Quantification::minCutIM(std::vector upper_bound) +{ + vector_double 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; +} + +void +FTAUtils::Quantification::computeRisk(int n, std::vector hazard) +{ + 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); + } +} + +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( + 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]; +} + +std::map +FTAUtils::Quantification::beIM(set_string cut_sets, + int n, + std::vector upper_bound, + std::vector & count_v) +{ + 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 (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; + available.push_back(is_a); + } + count_v.push_back(count); + if (count != 0) + { + std::vector fv, rrr, rir, rri, rii, bi; + + 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; + 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..ccfc21d31a --- /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" + +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) +{ + + set_string::iterator it_tree_mocus; + std::set tree_mocus_children, set_mocus, set_gold; + std::set::iterator it_tree_mocus_children; + + set_string matrix_gold = { + {"B2", "B1"}, {"B4", "B1"}, {"B1", "B5", "B3"}, {"B5", "B2", "B3"}, {"B4", "B5", "B3"}}; + + set_string::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++) + { + 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 + { + file_lists = {"not_a_valid_filename.txt", ""}; + + FTAUtils::FaultTree(file_lists[0], file_lists[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(); + 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 = {"logic2.txt", ""}; + + 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); + 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 new file mode 100644 index 0000000000..2295e91c84 --- /dev/null +++ b/unit/src/TestQuantificationUtils.C @@ -0,0 +1,459 @@ +// 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; +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_q = { + "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_q[0], + file_lists_q[1], + FTAUtils::Quantification::RISK, + file_lists_q[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_q = {"", "", ""}; + 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_q[0], + file_lists_q[1], + FTAUtils::Quantification::RISK, + file_lists_q[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_q = {"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_q[0], + file_lists_q[1], + FTAUtils::Quantification::FRAGILITY, + file_lists_q[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_q = {"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_q[0], + file_lists_q[1], + FTAUtils::Quantification::FRAGILITY, + file_lists_q[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_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_q[0], + file_lists_q[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_q = {"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_q[0], + file_lists_q[1], + FTAUtils::Quantification::RISK, + file_lists_q[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_q = {"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_q[0], + file_lists_q[1], + FTAUtils::Quantification::RISK, + file_lists_q[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_q = {"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_q[0], + file_lists_q[1], + FTAUtils::Quantification::RISK, + file_lists_q[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_q = {"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_q[0], + file_lists_q[1], + FTAUtils::Quantification::RISK, + file_lists_q[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"); +}