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,228 @@ 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
     {
+      /// 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<ErrorAggregate> vars{};
+      /// Aggregate of the combined error value at each time step
+      TemporalNormAggregate              total_error{"Total"};
+      /// Aggregate error for each variable
+      std::vector<TemporalNormAggregate> var_errors{};
 
+      /**
+       * @brief Construct with variable labels
+       */
       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()
+      /**
+       * @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 <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