diff --git a/GridKit/Testing/AggregateErrors.hpp b/GridKit/Testing/AggregateErrors.hpp index 2c03722eb..3a112fa96 100644 --- a/GridKit/Testing/AggregateErrors.hpp +++ b/GridKit/Testing/AggregateErrors.hpp @@ -2,6 +2,8 @@ #include #include +#include +#include #include #include @@ -11,96 +13,228 @@ namespace GridKit namespace Testing { + inline constexpr double DEFAULT_ABS_ERROR_THRESHOLD = + std::numeric_limits::epsilon(); + /** - * @brief Aggregate error data for a single variable + * @brief Aggregate norm data for a single variable over time */ - struct ErrorAggregate + struct TemporalNormAggregate { + /// Name of variable std::string label; - double max_error{0.0}; - double max_error_time{0.0}; - double error_norm_L2{0.0}; + /// L-inf norm + double max_value{0.0}; + /// Time at which max value occurred + double max_value_time{0.0}; + /// L2 norm + double L2{0.0}; - void push(double err, double time) + /** + * @brief Add a new value for the variable to be aggregated + */ + void push(double val, double time) { - if (err > max_error) + if (val > max_value) { - max_error = err; - max_error_time = time; + max_value = val; + max_value_time = time; } - error_norm_L2 += err * err; + L2 += val * val; } + /** + * @brief Finalize the calculation(s) + */ void wrap() { - error_norm_L2 = std::sqrt(error_norm_L2); + L2 = std::sqrt(L2); + } + + /** + * @brief Scale by a reference value (only if the current value is above + * the given threshold (used internally for relative errors) + */ + void scale(const TemporalNormAggregate& ref, double threshold) + { + if (ref.max_value > threshold) + { + max_value /= ref.max_value; + } + if (ref.L2 > threshold) + { + L2 /= ref.L2; + } } + /** + * @brief Pretty-print label and values to output stream + */ std::ostream& display( std::ostream& os = std::cout, const std::string& indent = " ") const { os << indent << label << ":\n" << indent << indent << "max : " - << std::format("{:.6e} (at time {:.3e})", max_error, max_error_time) + << std::format("{:.6e} (at time {:.3e})", max_value, max_value_time) << '\n' << indent << indent << "L2-norm : " - << std::format("{:.6e}", error_norm_L2) << '\n'; + << std::format("{:.6e}", L2) << '\n'; return os; } }; /** - * @brief A set of ErrorAggregate for each variable plus a total (combined) - * ErrorAggregate + * @brief A set of aggregate norms (represented with TemporalNormAggregate) + * for the error in each variable plus one for the total (combined) error * - * @note The "total" aggregate is based on the L-infinity norm of the local - * error of variables at a given time step. + * @note The "total_error" aggregate is based on the L-infinity norm of the + * local error of variables at a given time step. */ struct ErrorSet { - ErrorAggregate total{"Total"}; - std::vector vars{}; + /// Aggregate of the combined error value at each time step + TemporalNormAggregate total_error{"Total"}; + /// Aggregate error for each variable + std::vector var_errors{}; + /** + * @brief Construct with variable labels + */ template explicit ErrorSet(const C& labels) - : vars(std::size(labels)) + : var_errors(std::size(labels)) { for (std::size_t i = 0; i < std::size(labels); ++i) { - vars[i].label = labels[i]; + var_errors[i].label = labels[i]; } } - void push(const OutputAtTime& err) + virtual ~ErrorSet() { - for (std::size_t i = 0; i < vars.size(); ++i) - { - vars[i].push(std::abs(err.data[i]), err.t); - } - total.push(lInfNorm(err), err.t); } - void wrap() + /** + * @brief Finalize the calculations + */ + virtual void wrap() { - for (auto& agg : vars) + for (auto& agg : var_errors) { agg.wrap(); } - total.wrap(); + total_error.wrap(); } + /** + * @brief Take the error between the two output parameters for each + * variable and add to the aggregate + */ + virtual void push(const OutputAtTime&, const OutputAtTime&) = 0; + + /** + * @brief Pretty-print the set of errors for each variable and total + */ std::ostream& display(std::ostream& os = std::cout) const { std::string indent{" "}; os << "Error Set:\n"; - for (const auto& var : vars) + for (const auto& var : var_errors) { var.display(os, indent); } - total.display(os, indent); + total_error.display(os, indent); return os; } }; + enum class ErrorType + { + RELATIVE, + ABSOLUTE + }; + + struct RelativeError; + struct AbsoluteError; + + template + struct ErrorSetImpl; + + template <> + struct ErrorSetImpl : ErrorSet + { + template + explicit ErrorSetImpl( + const C& labels, + double abs_threshold = DEFAULT_ABS_ERROR_THRESHOLD) + : ErrorSet(labels), + abs_threshold_(abs_threshold), + ref_norms_(std::size(labels)) + { + } + + void push(const OutputAtTime& test, const OutputAtTime& ref) override + { + auto err = test - ref; + for (std::size_t i = 0; i < var_errors.size(); ++i) + { + var_errors[i].push(std::abs(err.data[i]), err.t); + ref_norms_[i].push(std::abs(ref.data[i]), ref.t); + } + total_error.push(lInfNorm(err), err.t); + ref_total_norm_.push(lInfNorm(ref), ref.t); + } + + void wrap() override + { + ErrorSet::wrap(); + for (std::size_t i = 0; i < var_errors.size(); ++i) + { + ref_norms_[i].wrap(); + var_errors[i].scale(ref_norms_[i], abs_threshold_); + } + ref_total_norm_.wrap(); + total_error.scale(ref_total_norm_, abs_threshold_); + } + + private: + double abs_threshold_{}; + TemporalNormAggregate ref_total_norm_{}; + std::vector ref_norms_{}; + }; + + template <> + struct ErrorSetImpl : ErrorSet + { + using ErrorSet::ErrorSet; + + void push(const OutputAtTime& test, const OutputAtTime& ref) override + { + auto err = test - ref; + for (std::size_t i = 0; i < var_errors.size(); ++i) + { + var_errors[i].push(std::abs(err.data[i]), err.t); + } + total_error.push(lInfNorm(err), err.t); + } + }; + + template + std::unique_ptr makeErrorSet( + ErrorType type, + const C& labels, + double abs_threshold = DEFAULT_ABS_ERROR_THRESHOLD) + { + switch (type) + { + case ErrorType::RELATIVE: + return std::make_unique>(labels, abs_threshold); + case ErrorType::ABSOLUTE: + return std::make_unique>(labels); + default: + throw std::runtime_error("Invalid error type"); + } + } + } // namespace Testing } // namespace GridKit diff --git a/GridKit/Testing/CSV.cpp b/GridKit/Testing/CSV.cpp index 22baefddb..aa161216f 100644 --- a/GridKit/Testing/CSV.cpp +++ b/GridKit/Testing/CSV.cpp @@ -35,50 +35,54 @@ namespace GridKit return ifs; } - ErrorSet compareCSV(const std::string& f_a, const std::string& f_b) + std::unique_ptr compareCSV(const std::string& f_test, + const std::string& f_ref, + ErrorType error_type, + double abs_threshold) { // Open files and read labels - auto ifs_a = checkOpenFile(f_a); - auto ifs_b = checkOpenFile(f_b); + auto ifs_test = checkOpenFile(f_test); + auto ifs_ref = checkOpenFile(f_ref); - std::string line_a; - std::getline(ifs_a, line_a); - auto labels_a = Tokenizer<>(line_a, ',')(); + std::string line_test; + std::getline(ifs_test, line_test); + auto labels_a = Tokenizer<>(line_test, ',')(); - std::string line_b; - std::getline(ifs_b, line_b); - auto labels_b = Tokenizer<>(line_b, ',')(); + std::string line_ref; + std::getline(ifs_ref, line_ref); + auto labels_b = Tokenizer<>(line_ref, ',')(); if (labels_a.size() != labels_b.size()) { - Log::error() << "Files \"" << f_a << "\" and \"" << f_b + Log::error() << "Files \"" << f_test << "\" and \"" << f_ref << "\" have different number of variables." << std::endl; } // Create error set - auto err = ErrorSet(std::span{next(begin(labels_a)), end(labels_a)}); + auto var_labels = std::span{next(begin(labels_a)), end(labels_a)}; + auto err = makeErrorSet(error_type, var_labels, abs_threshold); - while (ifs_a && ifs_b) + while (ifs_test && ifs_ref) { - std::getline(ifs_a, line_a); - std::getline(ifs_b, line_b); - if (!(ifs_a && ifs_b)) + std::getline(ifs_test, line_test); + std::getline(ifs_ref, line_ref); + if (!(ifs_test && ifs_ref)) { - if (ifs_a || ifs_b) + if (ifs_test || ifs_ref) { - Log::error() << "Files \"" << f_a << "\" and \"" << f_b + Log::error() << "Files \"" << f_test << "\" and \"" << f_ref << "\" have different lengths." << std::endl; } break; } - OutputAtTime d_a(Tokenizer(line_b, ',')()); - OutputAtTime d_b(Tokenizer(line_a, ',')()); + OutputAtTime d_test(Tokenizer(line_test, ',')()); + OutputAtTime d_ref(Tokenizer(line_ref, ',')()); - err.push(d_a - d_b); + err->push(d_test, d_ref); } - err.wrap(); + err->wrap(); return err; } diff --git a/GridKit/Testing/CSV.hpp b/GridKit/Testing/CSV.hpp index 462f31049..dadf5c99a 100644 --- a/GridKit/Testing/CSV.hpp +++ b/GridKit/Testing/CSV.hpp @@ -20,7 +20,11 @@ namespace GridKit * * @note This assumes a "time" variable is always in the first column. */ - ErrorSet compareCSV(const std::string& f_a, const std::string& f_b); + std::unique_ptr compareCSV( + const std::string& test_file, + const std::string& reference_file, + ErrorType error_type = ErrorType::RELATIVE, + double abs_threshold = DEFAULT_ABS_ERROR_THRESHOLD); } // namespace Testing } // namespace GridKit diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index 2930013e1..cb60b93c6 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -62,6 +62,10 @@ namespace GridKit fs::path reference_file; /// Error tolerance (between output file and reference file) double error_tol; + /// Type of total error (relative or absolute) + Testing::ErrorType error_type; + /// Smallest value at which to scale for relative error + double abs_err_threshold; /// Instance of model data SystemModelData<> model_data; }; @@ -107,6 +111,25 @@ namespace GridKit } c.error_tol = j.value("error_tolerance", 1.0e-4); + + using ErrorType = Testing::ErrorType; + if (j.contains("error_type")) + { + auto type_str = j.at("error_type").get(); + auto type_wrap = enum_cast(type_str, case_insensitive); + if (!type_wrap.has_value()) + { + Log::error() << "Invalid error type \"" << type_str << "\"; " + << "must be either \"relative\" or \"absolute\""; + } + c.error_type = type_wrap.value(); + } + else + { + c.error_type = ErrorType::RELATIVE; + } + + c.abs_err_threshold = j.value("abs_err_threshold", Testing::DEFAULT_ABS_ERROR_THRESHOLD); } /** @@ -223,16 +246,19 @@ namespace GridKit const auto& ref_file = study_data.reference_file; if (!out_file.empty() && !ref_file.empty()) { - auto errorSet = Testing::compareCSV(out_file, ref_file); + auto errorSet = Testing::compareCSV(out_file, + ref_file, + study_data.error_type, + study_data.abs_err_threshold); // Print the errors if (print_results) { - errorSet.display(); + errorSet->display(); } // Check against specified tolerance - status *= errorSet.total.max_error < study_data.error_tol; + status *= errorSet->total_error.max_value < study_data.error_tol; if (print_results) { diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 7f7d0f902..d6886267c 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -5,12 +5,14 @@ Name | Value ---------------------|------------------------------------------------------- `system_model_file` | Path to the system model file[^1] - `dt` | A floating point value for time step size - `tmax` | A floating point value for max time + `dt` | A floating-point value for time step size + `tmax` | A floating-point value for max time `events` | An array of event groups (see [Events](#events) below) - `output_file` | Path to output (CSV) file - `reference_file` | A string containing the name of the case - `error_tolerance` | A string containing the name of the case + `output_file` | Path to output (CSV) file (optional) + `reference_file` | A string containing the name of the case (optional) + `error_type` | One of { "relative" (default), "absolute" } + `error_tolerance` | A floating-point value for highest allowable total error (default: 1.0e-4) + `abs_err_threshold` | A floating-point value for the smallest value at which to scale relative error (default: machine epsilon for double-precision) [^1]: See system model [case format](../../GridKit/Model/PhasorDynamics/INPUT_FORMAT.md) diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.solver.json b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.solver.json index 37d0398ec..fa4d1f10f 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.solver.json +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.solver.json @@ -8,5 +8,6 @@ ], "output_file": "ThreeBus_six_cycle_fault_bus1.csv", "reference_file": "ThreeBusBasic.ref.csv", + "error_type": "absolute", "error_tolerance": 1e-4 } diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp index fe00fe68e..3b0ab5b9f 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp @@ -112,11 +112,11 @@ int main(int argc, const char* argv[]) auto errorSet = GridKit::Testing::compareCSV(out_file, ref_file); // Print the errors - errorSet.display(); + errorSet->display(); auto error_allowed = args["tolerance"].as(); - status *= errorSet.total.max_error < error_allowed; + status *= errorSet->total_error.max_value < error_allowed; status.report(); }