Use relative errors by default when comparing output files

GridKit/Testing/AggregateErrors.hpp

@@ -2,6 +2,8 @@
 
 #include <cmath>
 #include <format>
+#include <limits>
+#include <memory>
 #include <vector>
 
 #include <GridKit/Testing/OutputAtTime.hpp>
@@ -11,96 +13,196 @@ namespace GridKit
   namespace Testing
   {
 
+    inline constexpr double DEFAULT_ABS_ERROR_THRESHOLD =
+        std::numeric_limits<double>::epsilon();
+
     /**
-     * @brief Aggregate error data for a single variable
+     * @brief Aggregate norm data for a single variable over time
      */
-    struct ErrorAggregate
+    struct TemporalNormAggregate
     {
       std::string label;
-      double      max_error{0.0};
-      double      max_error_time{0.0};
-      double      error_norm_L2{0.0};
+      double      max{0.0};
+      double      max_time{0.0};
+      double      L2{0.0};
 
       void push(double err, double time)
       {
-        if (err > max_error)
+        if (err > max)
         {
-          max_error      = err;
-          max_error_time = time;
+          max      = err;
+          max_time = time;
         }
-        error_norm_L2 += err * err;
+        L2 += err * err;
       }
 
       void wrap()
       {
-        error_norm_L2 = std::sqrt(error_norm_L2);
+        L2 = std::sqrt(L2);
+      }
+
+      void scale(const TemporalNormAggregate& ref, double threshold)
+      {
+        if (ref.max > threshold)
+        {
+          max /= ref.max;
+        }
+        if (ref.L2 > threshold)
+        {
+          L2 /= ref.L2;
+        }
       }
 
       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, max_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
-     *
-     * @note The "total" 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<ErrorAggregate> vars{};
+      TemporalNormAggregate              total_error{"Total"};
+      std::vector<TemporalNormAggregate> var_errors{};
 
       template <typename C>
       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()
+      virtual void wrap()
       {
-        for (auto& agg : vars)
+        for (auto& agg : var_errors)
         {
           agg.wrap();
         }
-        total.wrap();
+        total_error.wrap();
       }
 
+      virtual void push(const OutputAtTime&, const OutputAtTime&) = 0;
+
       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;
+
+    /**
+     * @brief A set of TemporalNormAggregate for the error in each variable plus
+     * a TemporalNormAggregate representing total (combined) error
+     *
+     * @note The "total_error" aggregate is based on the L-infinity norm of the local
+     * error of variables at a given time step.
+     */
+    template <typename error_type = RelativeError>
+    struct ErrorSetImpl;
+
+    template <>
+    struct ErrorSetImpl<RelativeError> : ErrorSet
+    {
+      template <typename C>
+      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<TemporalNormAggregate> ref_norms_{};
+    };
+
+    template <>
+    struct ErrorSetImpl<AbsoluteError> : 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 <typename C>
+    std::unique_ptr<ErrorSet> makeErrorSet(
+        ErrorType type,
+        const C&  labels,
+        double    abs_threshold = DEFAULT_ABS_ERROR_THRESHOLD)
+    {
+      switch (type)
+      {
+      case ErrorType::RELATIVE:
+        return std::make_unique<ErrorSetImpl<RelativeError>>(labels, abs_threshold);
+      case ErrorType::ABSOLUTE:
+        return std::make_unique<ErrorSetImpl<AbsoluteError>>(labels);
+      default:
+        throw std::runtime_error("Invalid error type");
+      }
+    }
+
   } // namespace Testing
 } // namespace GridKit

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<ErrorSet> 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<double>(line_b, ',')());
-        OutputAtTime d_b(Tokenizer<double>(line_a, ',')());
+        OutputAtTime d_test(Tokenizer<double>(line_test, ',')());
+        OutputAtTime d_ref(Tokenizer<double>(line_ref, ',')());
 
-        err.push(d_a - d_b);
+        err->push(d_test, d_ref);
       }
 
-      err.wrap();
+      err->wrap();
 
       return err;
     }

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<ErrorSet> 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

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<std::string>();
+        auto type_wrap = enum_cast<ErrorType>(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);
     }
 
     /**
@@ -220,16 +243,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 < study_data.error_tol;
 
         if (print_results)
         {