FESADev/include/fesa/Results/ReferenceComparison.hpp

#pragma once

#include "fesa/Results/ResultModel.hpp"
#include "fesa/Util/Util.hpp"

#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <istream>
#include <limits>
#include <map>
#include <sstream>
#include <string>
#include <vector>

namespace fesa {

struct CsvDisplacementRow {
  GlobalId node_id = 0;
  std::array<Real, 6> values{};
};

struct CsvDisplacementTable {
  std::map<GlobalId, CsvDisplacementRow> rows;
  std::vector<Diagnostic> diagnostics;
};

inline std::vector<std::string> displacementCsvRequiredColumns() {
  return {"Node Label", "U-U1", "U-U2", "U-U3", "UR-UR1", "UR-UR2", "UR-UR3"};
}

inline CsvDisplacementTable loadDisplacementCsvFromStream(std::istream& input, const std::string& source_name) {
  CsvDisplacementTable table;
  std::string line;
  if (!std::getline(input, line)) {
    table.diagnostics.push_back({Severity::Error, "FESA-CSV-EMPTY", "Displacement CSV is empty", {source_name, 1, ""}});
    return table;
  }
  const std::vector<std::string> required = displacementCsvRequiredColumns();
  std::vector<std::string> headers = splitCsv(line);
  std::map<std::string, std::size_t> column;
  for (std::size_t i = 0; i < headers.size(); ++i) {
    column[trim(headers[i])] = i;
  }
  for (const std::string& name : required) {
    if (column.count(name) == 0) {
      table.diagnostics.push_back({Severity::Error, "FESA-CSV-MISSING-COLUMN", "Missing CSV column: " + name, {source_name, 1, ""}});
    }
  }
  if (hasError(table.diagnostics)) {
    return table;
  }
  LocalIndex line_number = 1;
  while (std::getline(input, line)) {
    ++line_number;
    if (trim(line).empty()) {
      continue;
    }
    std::vector<std::string> fields = splitCsv(line);
    auto get = [&](const std::string& name) -> std::string {
      const std::size_t index = column[name];
      return index < fields.size() ? fields[index] : "";
    };
    auto node_id = parseInt64(get("Node Label"));
    if (!node_id) {
      table.diagnostics.push_back({Severity::Error, "FESA-CSV-NODE", "Invalid node label", {source_name, line_number, ""}});
      continue;
    }
    if (table.rows.count(*node_id) != 0) {
      table.diagnostics.push_back({Severity::Error, "FESA-CSV-DUPLICATE-NODE", "Duplicate node label", {source_name, line_number, ""}});
      continue;
    }
    CsvDisplacementRow row;
    row.node_id = *node_id;
    for (std::size_t i = 0; i < 6; ++i) {
      auto value = parseReal(get(required[i + 1]));
      if (!value) {
        table.diagnostics.push_back({Severity::Error, "FESA-CSV-NUMERIC", "Invalid displacement value", {source_name, line_number, ""}});
        value = 0.0;
      }
      row.values[i] = *value;
    }
    table.rows[*node_id] = row;
  }
  return table;
}

inline CsvDisplacementTable loadDisplacementCsvFromString(const std::string& text, const std::string& source_name = "<memory>") {
  std::istringstream input(text);
  return loadDisplacementCsvFromStream(input, source_name);
}

inline CsvDisplacementTable loadDisplacementCsv(const std::string& path) {
  std::ifstream input(path);
  if (!input.good()) {
    CsvDisplacementTable table;
    table.diagnostics.push_back({Severity::Error, "FESA-CSV-READ", "Could not read displacement CSV", {path, 0, ""}});
    return table;
  }
  return loadDisplacementCsvFromStream(input, path);
}

struct ComparisonOptions {
  Real abs_tol = 1.0e-12;
  Real rel_tol = 1.0e-5;
  Real reference_scale = 1.0;
};

struct ComparisonResult {
  bool pass = false;
  Real max_abs_error = 0.0;
  Real max_rel_error = 0.0;
  std::vector<Diagnostic> diagnostics;
};

inline ComparisonResult compareDisplacements(const FieldOutput& actual,
                                             const CsvDisplacementTable& expected,
                                             ComparisonOptions options = {}) {
  ComparisonResult result;
  result.diagnostics = expected.diagnostics;
  if (hasError(result.diagnostics)) {
    return result;
  }
  if (actual.name != "U") {
    result.diagnostics.push_back({Severity::Error, "FESA-COMPARE-FIELD-NAME", "Expected FESA displacement field named U", {}});
  }
  if (actual.component_labels != displacementComponentLabels()) {
    result.diagnostics.push_back({Severity::Error,
                                  "FESA-COMPARE-COMPONENT-LABELS",
                                  "FESA U field component labels must be UX,UY,UZ,RX,RY,RZ",
                                  {}});
  }
  if (actual.position != "NODAL" || actual.entity_type != "node" || actual.basis != "GLOBAL") {
    result.diagnostics.push_back({Severity::Error, "FESA-COMPARE-FIELD-METADATA", "FESA U field must be nodal values in the global basis", {}});
  }
  if (actual.entity_ids.size() != actual.values.size()) {
    result.diagnostics.push_back({Severity::Error, "FESA-COMPARE-FIELD-SIZE", "FESA U field entity/value counts differ", {}});
  }
  std::map<GlobalId, std::array<Real, 6>> actual_by_node;
  const std::size_t actual_count = std::min(actual.entity_ids.size(), actual.values.size());
  for (std::size_t i = 0; i < actual_count; ++i) {
    if (actual_by_node.count(actual.entity_ids[i]) != 0) {
      result.diagnostics.push_back(
          {Severity::Error, "FESA-COMPARE-DUPLICATE-ACTUAL", "FESA U field contains duplicate node " + std::to_string(actual.entity_ids[i]), {}});
      continue;
    }
    actual_by_node[actual.entity_ids[i]] = actual.values[i];
  }
  for (const auto& [node_id, row] : expected.rows) {
    auto actual_it = actual_by_node.find(node_id);
    if (actual_it == actual_by_node.end()) {
      result.diagnostics.push_back({Severity::Error, "FESA-COMPARE-MISSING-ACTUAL", "FESA U field is missing node " + std::to_string(node_id), {}});
      continue;
    }
    for (std::size_t component = 0; component < 6; ++component) {
      const Real expected_value = row.values[component];
      const Real actual_value = actual_it->second[component];
      const Real abs_error = std::fabs(actual_value - expected_value);
      const Real scale = std::max(std::fabs(expected_value), std::fabs(options.reference_scale));
      const Real rel_error = scale > 0.0 ? abs_error / scale : (abs_error == 0.0 ? 0.0 : std::numeric_limits<Real>::infinity());
      result.max_abs_error = std::max(result.max_abs_error, abs_error);
      result.max_rel_error = std::max(result.max_rel_error, rel_error);
      if (!(abs_error <= options.abs_tol || rel_error <= options.rel_tol)) {
        const std::string component_label = displacementComponentLabels()[component];
        result.diagnostics.push_back({Severity::Error,
                                      "FESA-COMPARE-TOLERANCE",
                                      "Displacement comparison failed at node " + std::to_string(node_id) + " component " + component_label,
                                      {}});
      }
    }
  }
  result.pass = !hasError(result.diagnostics);
  return result;
}

}  // namespace fesa