#include "mitc4_reference_comparison.hpp"

#include <hdf5.h>

#include <algorithm>
#include <array>
#include <cerrno>
#include <cctype>
#include <charconv>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <locale>
#include <map>
#include <sstream>
#include <stdexcept>
#include <string>
#include <system_error>
#include <utility>
#include <vector>

namespace fesa::test {
namespace {

constexpr const char* kDisplacementPath =
    "/steps/Step-1/frames/0/nodal/displacement";
constexpr const char* kInternalFormulation = "FESA-MITC4";
constexpr const char* kIntegrationRule =
    "2x2x2-gauss; mitc4-edge-midpoint-shear";
constexpr double kFixedAbsoluteTolerance = 1.0e-5;
constexpr std::array<const char*, 6> kComponents{
    "U1", "U2", "U3", "UR1", "UR2", "UR3"};
const std::vector<std::string> kExpectedHeader{
    "Part Instance Name", "Node Label", "U-U1", "U-U2", "U-U3",
    "UR-UR1", "UR-UR2", "UR-UR3"};

class ComparisonFailure final : public std::runtime_error {
public:
    ComparisonFailure(std::string code, std::string message)
        : std::runtime_error{std::move(message)}, code_{std::move(code)} {}

    const std::string& code() const noexcept { return code_; }

private:
    std::string code_;
};

[[noreturn]] void fail(const std::string& code, const std::string& message) {
    throw ComparisonFailure{code, message};
}

Status failureStatus(
    const std::string& caseId,
    const std::string& code,
    const std::string& message) {
    return Status::failure(
        FailureCategory::model,
        {{Severity::error, code, {}, "", caseId, message}});
}

std::string trim(const std::string& value) {
    const auto isSpace = [](const unsigned char character) {
        return std::isspace(character) != 0;
    };
    const auto begin = std::find_if_not(
        value.begin(), value.end(), [&](const char character) {
            return isSpace(static_cast<unsigned char>(character));
        });
    const auto end = std::find_if_not(
        value.rbegin(), value.rend(), [&](const char character) {
            return isSpace(static_cast<unsigned char>(character));
        }).base();
    return begin < end ? std::string{begin, end} : std::string{};
}

std::vector<std::string> splitCsvLine(const std::string& line) {
    std::vector<std::string> fields;
    std::size_t start = 0U;
    while (true) {
        const std::size_t comma = line.find(',', start);
        fields.push_back(trim(line.substr(start, comma - start)));
        if (comma == std::string::npos) {
            break;
        }
        start = comma + 1U;
    }
    return fields;
}

std::int64_t parsePositiveLabel(const std::string& field) {
    std::int64_t value = 0;
    const char* const begin = field.data();
    const char* const end = begin + field.size();
    const auto parsed = std::from_chars(begin, end, value);
    if (parsed.ec != std::errc{} || parsed.ptr != end || value <= 0) {
        fail("schema-mismatch", "A source-node label is invalid.");
    }
    return value;
}

double parseFiniteDouble(const std::string& field) {
    if (field.empty()) {
        fail("schema-mismatch", "A displacement CSV numeric field is empty.");
    }
    errno = 0;
    char* end = nullptr;
    const double value = std::strtod(field.c_str(), &end);
    if (errno == ERANGE || end == field.c_str() || end == nullptr ||
        *end != '\0' || !std::isfinite(value)) {
        fail(
            "schema-mismatch",
            "A displacement CSV numeric field is invalid or nonfinite.");
    }
    return value;
}

std::string uppercaseAscii(std::string value) {
    std::transform(
        value.begin(), value.end(), value.begin(), [](const char character) {
            return character >= 'a' && character <= 'z'
                       ? static_cast<char>(character - 'a' + 'A')
                       : character;
        });
    return value;
}

struct IdentityKey {
    std::string instanceName;
    std::int64_t sourceNodeLabel;

    bool operator<(const IdentityKey& other) const {
        const std::string normalizedInstance = uppercaseAscii(instanceName);
        const std::string normalizedOther = uppercaseAscii(other.instanceName);
        if (normalizedInstance != normalizedOther) {
            return normalizedInstance < normalizedOther;
        }
        return sourceNodeLabel < other.sourceNodeLabel;
    }
};

struct WideRow {
    IdentityKey identity;
    std::array<double, 6> values;
};

std::vector<WideRow> readReferenceCsv(const std::filesystem::path& path) {
    std::ifstream stream{path};
    if (!stream) {
        fail(
            "needs-reference-artifacts",
            "The declared displacement CSV is missing or unreadable.");
    }
    std::string line;
    if (!std::getline(stream, line)) {
        fail("schema-mismatch", "The declared displacement CSV is empty.");
    }
    if (!line.empty() && line.back() == '\r') {
        line.pop_back();
    }
    if (splitCsvLine(line) != kExpectedHeader) {
        fail(
            "schema-mismatch",
            "The displacement CSV header does not match the six-component contract.");
    }

    std::vector<WideRow> rows;
    std::map<IdentityKey, std::size_t> identities;
    while (std::getline(stream, line)) {
        if (!line.empty() && line.back() == '\r') {
            line.pop_back();
        }
        if (line.empty()) {
            fail("schema-mismatch", "Blank displacement CSV rows are not allowed.");
        }
        const auto fields = splitCsvLine(line);
        if (fields.size() != kExpectedHeader.size() || fields[0U].empty()) {
            fail("schema-mismatch", "A displacement CSV row has invalid schema.");
        }
        WideRow row{};
        row.identity.instanceName = fields[0U];
        row.identity.sourceNodeLabel = parsePositiveLabel(fields[1U]);
        for (std::size_t component = 0U;
             component < row.values.size();
             ++component) {
            row.values[component] = parseFiniteDouble(fields[component + 2U]);
        }
        if (!identities.emplace(row.identity, rows.size()).second) {
            fail("schema-mismatch", "A displacement CSV row identity is duplicated.");
        }
        rows.push_back(std::move(row));
    }
    if (rows.empty()) {
        fail("schema-mismatch", "The displacement CSV contains no data rows.");
    }
    return rows;
}

class Hdf5Handle {
public:
    using Closer = herr_t (*)(hid_t);

    Hdf5Handle() = default;
    Hdf5Handle(const hid_t value, Closer closer)
        : value_{value}, closer_{closer} {}
    Hdf5Handle(const Hdf5Handle&) = delete;
    Hdf5Handle& operator=(const Hdf5Handle&) = delete;
    Hdf5Handle(Hdf5Handle&& other) noexcept
        : value_{other.value_}, closer_{other.closer_} {
        other.value_ = -1;
        other.closer_ = nullptr;
    }
    ~Hdf5Handle() {
        if (value_ >= 0 && closer_ != nullptr) {
            (void)closer_(value_);
        }
    }

    hid_t get() const noexcept { return value_; }

private:
    hid_t value_{-1};
    Closer closer_{nullptr};
};

class Hdf5ErrorSilencer {
public:
    Hdf5ErrorSilencer() {
        if (H5Eget_auto2(H5E_DEFAULT, &callback_, &clientData_) >= 0 &&
            H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr) >= 0) {
            active_ = true;
        }
    }
    Hdf5ErrorSilencer(const Hdf5ErrorSilencer&) = delete;
    Hdf5ErrorSilencer& operator=(const Hdf5ErrorSilencer&) = delete;
    ~Hdf5ErrorSilencer() {
        if (active_) {
            (void)H5Eset_auto2(H5E_DEFAULT, callback_, clientData_);
        }
    }

private:
    H5E_auto2_t callback_{nullptr};
    void* clientData_{nullptr};
    bool active_{false};
};

class Hdf5VlenReclaimer {
public:
    Hdf5VlenReclaimer(
        const hid_t memoryType,
        const hid_t dataSpace,
        void* const data) noexcept
        : memoryType_{memoryType}, dataSpace_{dataSpace}, data_{data} {}
    Hdf5VlenReclaimer(const Hdf5VlenReclaimer&) = delete;
    Hdf5VlenReclaimer& operator=(const Hdf5VlenReclaimer&) = delete;
    ~Hdf5VlenReclaimer() {
        if (active_) {
            (void)H5Dvlen_reclaim(
                memoryType_, dataSpace_, H5P_DEFAULT, data_);
        }
    }

    void reclaim() {
        active_ = false;
        if (H5Dvlen_reclaim(memoryType_, dataSpace_, H5P_DEFAULT, data_) < 0) {
            fail("schema-mismatch", "Unable to reclaim HDF5 variable strings.");
        }
    }

private:
    hid_t memoryType_;
    hid_t dataSpace_;
    void* data_;
    bool active_{true};
};

hid_t requireId(const hid_t value, const char* message) {
    if (value < 0) {
        fail("schema-mismatch", message);
    }
    return value;
}

void requireHdf(const herr_t value, const char* message) {
    if (value < 0) {
        fail("schema-mismatch", message);
    }
}

Hdf5Handle makeUtf8StringType() {
    Hdf5Handle type{
        requireId(H5Tcopy(H5T_C_S1), "Unable to copy an HDF5 string type."),
        H5Tclose};
    requireHdf(
        H5Tset_size(type.get(), H5T_VARIABLE),
        "Unable to define an HDF5 variable string type.");
    requireHdf(
        H5Tset_cset(type.get(), H5T_CSET_UTF8),
        "Unable to define an HDF5 UTF-8 string type.");
    return type;
}

Hdf5Handle openDataset(const hid_t file, const char* path) {
    return {
        requireId(
            H5Dopen2(file, path, H5P_DEFAULT),
            "A required HDF5 dataset is missing."),
        H5Dclose};
}

std::vector<hsize_t> dimensions(const hid_t dataset) {
    Hdf5Handle space{
        requireId(H5Dget_space(dataset), "Unable to inspect HDF5 dimensions."),
        H5Sclose};
    const int rank = H5Sget_simple_extent_ndims(space.get());
    if (rank < 0) {
        fail("schema-mismatch", "Unable to inspect HDF5 rank.");
    }
    std::vector<hsize_t> result(static_cast<std::size_t>(rank));
    if (rank > 0) {
        requireHdf(
            H5Sget_simple_extent_dims(space.get(), result.data(), nullptr),
            "Unable to inspect HDF5 extents.");
    }
    return result;
}

std::string readStringAttribute(const hid_t object, const char* name) {
    Hdf5Handle attribute{
        requireId(
            H5Aopen(object, name, H5P_DEFAULT),
            "A required HDF5 string attribute is missing."),
        H5Aclose};
    Hdf5Handle type{
        requireId(H5Aget_type(attribute.get()), "Unable to inspect an attribute."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_STRING ||
        H5Tis_variable_str(type.get()) <= 0 ||
        H5Tget_cset(type.get()) != H5T_CSET_UTF8) {
        fail("schema-mismatch", "An HDF5 string attribute has the wrong type.");
    }
    char* raw = nullptr;
    requireHdf(
        H5Aread(attribute.get(), type.get(), &raw),
        "Unable to read an HDF5 string attribute.");
    if (raw == nullptr) {
        fail("schema-mismatch", "An HDF5 string attribute is null.");
    }
    const std::string result{raw};
    requireHdf(H5free_memory(raw), "Unable to free HDF5 attribute memory.");
    return result;
}

std::uint64_t readUint64Attribute(const hid_t object, const char* name) {
    Hdf5Handle attribute{
        requireId(
            H5Aopen(object, name, H5P_DEFAULT),
            "A required HDF5 integer attribute is missing."),
        H5Aclose};
    Hdf5Handle type{
        requireId(H5Aget_type(attribute.get()), "Unable to inspect an attribute."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_INTEGER ||
        H5Tget_size(type.get()) != sizeof(std::uint64_t) ||
        H5Tget_sign(type.get()) != H5T_SGN_NONE) {
        fail("schema-mismatch", "An HDF5 integer attribute has the wrong type.");
    }
    std::uint64_t result = 0U;
    requireHdf(
        H5Aread(attribute.get(), H5T_NATIVE_UINT64, &result),
        "Unable to read an HDF5 integer attribute.");
    return result;
}

void requireStringAttribute(
    const hid_t object, const char* name, const std::string& expected) {
    if (readStringAttribute(object, name) != expected) {
        fail("schema-mismatch", "An HDF5 string attribute has the wrong value.");
    }
}

void requireExactCompoundMembers(
    const hid_t dataset, const std::vector<const char*>& expected) {
    Hdf5Handle type{
        requireId(H5Dget_type(dataset), "Unable to inspect an HDF5 compound type."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_COMPOUND ||
        H5Tget_nmembers(type.get()) != static_cast<int>(expected.size())) {
        fail("schema-mismatch", "An HDF5 compound schema has the wrong member count.");
    }
    for (std::size_t index = 0U; index < expected.size(); ++index) {
        char* name = H5Tget_member_name(type.get(), static_cast<unsigned>(index));
        if (name == nullptr) {
            fail("schema-mismatch", "Unable to inspect an HDF5 compound member.");
        }
        const std::string actual{name};
        requireHdf(H5free_memory(name), "Unable to free HDF5 member memory.");
        if (actual != expected[index]) {
            fail("schema-mismatch", "An HDF5 compound member is out of contract order.");
        }
    }
}

struct NodeReadRow {
    std::uint64_t internalNodeId;
    char* instanceName;
    char* sourceLabel;
    double coordinates[3];
};

std::vector<WideRow> readNodes(const hid_t file) {
    auto dataset = openDataset(file, "/model/nodes");
    const auto shape = dimensions(dataset.get());
    if (shape.size() != 1U || shape[0U] == 0U ||
        shape[0U] > static_cast<hsize_t>((std::numeric_limits<std::size_t>::max)())) {
        fail("schema-mismatch", "The HDF5 node dataset has an invalid shape.");
    }
    requireExactCompoundMembers(
        dataset.get(),
        {"internal_node_id", "instance_name", "source_label", "coordinates"});
    requireStringAttribute(dataset.get(), "coordinate_system", "global-cartesian");
    requireStringAttribute(dataset.get(), "units_label", "length");

    auto stringType = makeUtf8StringType();
    const hsize_t coordinateDimensions[] = {3U};
    Hdf5Handle coordinates{
        requireId(
            H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinateDimensions),
            "Unable to create node coordinate memory type."),
        H5Tclose};
    Hdf5Handle memoryType{
        requireId(
            H5Tcreate(H5T_COMPOUND, sizeof(NodeReadRow)),
            "Unable to create node memory type."),
        H5Tclose};
    requireHdf(
        H5Tinsert(
            memoryType.get(), "internal_node_id",
            HOFFSET(NodeReadRow, internalNodeId), H5T_NATIVE_UINT64),
        "Unable to define node ID memory field.");
    requireHdf(
        H5Tinsert(
            memoryType.get(), "instance_name", HOFFSET(NodeReadRow, instanceName),
            stringType.get()),
        "Unable to define node instance memory field.");
    requireHdf(
        H5Tinsert(
            memoryType.get(), "source_label", HOFFSET(NodeReadRow, sourceLabel),
            stringType.get()),
        "Unable to define node label memory field.");
    requireHdf(
        H5Tinsert(
            memoryType.get(), "coordinates", HOFFSET(NodeReadRow, coordinates),
            coordinates.get()),
        "Unable to define node coordinate memory field.");

    const std::size_t count = static_cast<std::size_t>(shape[0U]);
    std::vector<NodeReadRow> raw(count);
    requireHdf(
        H5Dread(
            dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL, H5P_DEFAULT,
            raw.data()),
        "Unable to read HDF5 node identities.");
    Hdf5Handle space{
        requireId(H5Dget_space(dataset.get()), "Unable to inspect node space."),
        H5Sclose};
    Hdf5VlenReclaimer reclaimer{memoryType.get(), space.get(), raw.data()};
    std::vector<WideRow> rows;
    rows.reserve(count);
    std::map<IdentityKey, std::size_t> identities;
    for (std::size_t index = 0U; index < count; ++index) {
        if (raw[index].internalNodeId != index || raw[index].instanceName == nullptr ||
            raw[index].sourceLabel == nullptr || raw[index].instanceName[0] == '\0' ||
            !std::all_of(
                std::begin(raw[index].coordinates),
                std::end(raw[index].coordinates),
                [](const double value) { return std::isfinite(value); })) {
            fail("schema-mismatch", "An HDF5 node row has invalid identity or values.");
        }
        WideRow row{};
        row.identity.instanceName = raw[index].instanceName;
        row.identity.sourceNodeLabel = parsePositiveLabel(raw[index].sourceLabel);
        if (!identities.emplace(row.identity, rows.size()).second) {
            fail("schema-mismatch", "An HDF5 node identity is duplicated.");
        }
        rows.push_back(std::move(row));
    }
    reclaimer.reclaim();
    return rows;
}

struct ElementIdentityReadRow {
    char* sourceElementType;
    char* internalFormulation;
};

void requireElementIdentity(
    const hid_t file, const std::string& expectedSourceType) {
    auto dataset = openDataset(file, "/model/elements");
    const auto shape = dimensions(dataset.get());
    if (shape.size() != 1U || shape[0U] == 0U ||
        shape[0U] > static_cast<hsize_t>((std::numeric_limits<std::size_t>::max)())) {
        fail("schema-mismatch", "The HDF5 element dataset has an invalid shape.");
    }
    requireExactCompoundMembers(
        dataset.get(),
        {"internal_element_id", "instance_name", "source_label",
         "source_element_type", "internal_formulation", "node_internal_ids",
         "shell_section_internal_id", "material_internal_id"});
    requireStringAttribute(dataset.get(), "formulation", kInternalFormulation);

    auto stringType = makeUtf8StringType();
    Hdf5Handle memoryType{
        requireId(
            H5Tcreate(H5T_COMPOUND, sizeof(ElementIdentityReadRow)),
            "Unable to create element identity memory type."),
        H5Tclose};
    requireHdf(
        H5Tinsert(
            memoryType.get(), "source_element_type",
            HOFFSET(ElementIdentityReadRow, sourceElementType), stringType.get()),
        "Unable to define source element type memory field.");
    requireHdf(
        H5Tinsert(
            memoryType.get(), "internal_formulation",
            HOFFSET(ElementIdentityReadRow, internalFormulation), stringType.get()),
        "Unable to define formulation memory field.");
    const std::size_t count = static_cast<std::size_t>(shape[0U]);
    std::vector<ElementIdentityReadRow> rows(count);
    requireHdf(
        H5Dread(
            dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL, H5P_DEFAULT,
            rows.data()),
        "Unable to read HDF5 element identity.");
    Hdf5Handle space{
        requireId(H5Dget_space(dataset.get()), "Unable to inspect element space."),
        H5Sclose};
    Hdf5VlenReclaimer reclaimer{memoryType.get(), space.get(), rows.data()};
    for (const auto& row : rows) {
        if (row.sourceElementType == nullptr || row.internalFormulation == nullptr ||
            row.sourceElementType != expectedSourceType ||
            row.internalFormulation != std::string{kInternalFormulation}) {
            fail(
                "schema-mismatch",
                "The HDF5 source element type or internal formulation is invalid.");
        }
    }
    reclaimer.reclaim();
}

std::vector<double> readDisplacement(
    const hid_t file, const std::size_t nodeCount) {
    auto dataset = openDataset(file, kDisplacementPath);
    if (dimensions(dataset.get()) !=
        std::vector<hsize_t>{static_cast<hsize_t>(nodeCount), 6U}) {
        fail("schema-mismatch", "The HDF5 displacement dataset has the wrong shape.");
    }
    Hdf5Handle type{
        requireId(H5Dget_type(dataset.get()), "Unable to inspect displacement type."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_FLOAT ||
        H5Tget_size(type.get()) != sizeof(double) ||
        H5Tequal(type.get(), H5T_IEEE_F64LE) <= 0) {
        fail("schema-mismatch", "The HDF5 displacement dataset is not float64 LE.");
    }
    requireStringAttribute(dataset.get(), "component_names", "UX,UY,UZ,URX,URY,URZ");
    requireStringAttribute(
        dataset.get(), "component_unit_dimensions",
        "length,length,length,radian,radian,radian");
    requireStringAttribute(dataset.get(), "coordinate_system", "global-cartesian");
    requireStringAttribute(dataset.get(), "location", "nodal");
    requireStringAttribute(dataset.get(), "step_name", "Step-1");
    if (readUint64Attribute(dataset.get(), "frame_index") != 0U) {
        fail("schema-mismatch", "The HDF5 displacement frame identity is invalid.");
    }
    std::vector<double> values(nodeCount * kComponents.size());
    if (!values.empty()) {
        requireHdf(
            H5Dread(
                dataset.get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                H5P_DEFAULT, values.data()),
            "Unable to read HDF5 displacement values.");
    }
    if (!std::all_of(values.begin(), values.end(), [](const double value) {
            return std::isfinite(value);
        })) {
        fail("schema-mismatch", "An HDF5 displacement value is nonfinite.");
    }
    return values;
}

struct Hdf5Projection {
    std::vector<WideRow> rows;
    std::string sourceElementType;
    std::string internalFormulation;
    std::string integrationRule;
};

Hdf5Projection readHdf5(const Mitc4ReferenceCase& referenceCase) {
    Hdf5ErrorSilencer silencer;
    Hdf5Handle file{
        requireId(
            H5Fopen(
                referenceCase.resultsHdf5Path.string().c_str(), H5F_ACC_RDONLY,
                H5P_DEFAULT),
            "The authoritative HDF5 output cannot be opened."),
        H5Fclose};
    Hdf5Handle metadata{
        requireId(
            H5Gopen2(file.get(), "/metadata", H5P_DEFAULT),
            "The HDF5 metadata group is missing."),
        H5Gclose};
    if (readUint64Attribute(metadata.get(), "schema_version") != 0U ||
        readUint64Attribute(metadata.get(), "frame_index") != 0U) {
        fail("schema-mismatch", "The HDF5 schema or frame version is invalid.");
    }
    requireStringAttribute(metadata.get(), "feature_id", "linear-static-mitc4-shell");
    requireStringAttribute(metadata.get(), "step_name", "Step-1");
    requireStringAttribute(
        metadata.get(), "internal_formulation", kInternalFormulation);
    requireStringAttribute(metadata.get(), "integration_rule", kIntegrationRule);
    const std::string normalizedInput =
        std::filesystem::absolute(referenceCase.inputPath)
            .lexically_normal()
            .generic_u8string();
    const std::string sourceIdentity =
        readStringAttribute(metadata.get(), "source_input_identity");
    if (sourceIdentity.find("path=" + normalizedInput + ";content_identity=") != 0U) {
        fail(
            "schema-mismatch",
            "The HDF5 source input identity does not match the declared input.");
    }
    auto rows = readNodes(file.get());
    requireElementIdentity(file.get(), referenceCase.expectedSourceElementType);
    const auto values = readDisplacement(file.get(), rows.size());
    for (std::size_t row = 0U; row < rows.size(); ++row) {
        std::copy_n(
            values.begin() + static_cast<std::ptrdiff_t>(row * kComponents.size()),
            kComponents.size(), rows[row].values.begin());
    }
    return {
        std::move(rows), referenceCase.expectedSourceElementType,
        kInternalFormulation, kIntegrationRule};
}

void requireArtifacts(const Mitc4ReferenceCase& referenceCase) {
    if (referenceCase.caseId.empty() ||
        (referenceCase.expectedSourceElementType != "S4" &&
         referenceCase.expectedSourceElementType != "S4R")) {
        fail("schema-mismatch", "The MITC4 reference case identity is invalid.");
    }
    std::error_code error;
    for (const auto* path : {
             &referenceCase.inputPath,
             &referenceCase.displacementCsvPath,
             &referenceCase.resultsHdf5Path}) {
        if (!std::filesystem::is_regular_file(*path, error) || error) {
            fail(
                "needs-reference-artifacts",
                "A declared MITC4 input, displacement CSV, or HDF5 file is missing.");
        }
    }
}

std::string finiteText(const double value) {
    std::ostringstream stream;
    stream.imbue(std::locale::classic());
    stream << std::setprecision(std::numeric_limits<double>::max_digits10)
           << value;
    return stream.str();
}

std::string jsonEscape(const std::string& value) {
    std::ostringstream stream;
    for (const unsigned char character : value) {
        switch (character) {
        case '"': stream << "\\\""; break;
        case '\\': stream << "\\\\"; break;
        case '\b': stream << "\\b"; break;
        case '\f': stream << "\\f"; break;
        case '\n': stream << "\\n"; break;
        case '\r': stream << "\\r"; break;
        case '\t': stream << "\\t"; break;
        default:
            if (character < 0x20U) {
                stream << "\\u00" << std::hex << std::setw(2)
                       << std::setfill('0') << static_cast<unsigned>(character)
                       << std::dec << std::setfill(' ');
            } else {
                stream << static_cast<char>(character);
            }
            break;
        }
    }
    return stream.str();
}

} // namespace

Result<Mitc4ComparisonReport> Mitc4ReferenceComparison::compare(
    const Mitc4ReferenceCase& referenceCase) {
    try {
        requireArtifacts(referenceCase);
        const auto referenceRows =
            readReferenceCsv(referenceCase.displacementCsvPath);
        auto hdf5 = readHdf5(referenceCase);

        std::map<IdentityKey, const WideRow*> referenceByIdentity;
        for (const auto& row : referenceRows) {
            referenceByIdentity.emplace(row.identity, &row);
        }
        if (referenceRows.size() != hdf5.rows.size()) {
            fail(
                "schema-mismatch",
                "The HDF5 and CSV source-node row counts are not equal.");
        }
        for (const auto& row : hdf5.rows) {
            if (referenceByIdentity.find(row.identity) == referenceByIdentity.end()) {
                fail(
                    "schema-mismatch",
                    "The HDF5 and CSV source-node identities are not equal.");
            }
        }

        std::array<double, 6> scales{};
        for (const auto& row : referenceRows) {
            for (std::size_t component = 0U;
                 component < scales.size();
                 ++component) {
                scales[component] =
                    (std::max)(scales[component], std::abs(row.values[component]));
            }
        }

        Mitc4ComparisonReport report{};
        report.caseId = referenceCase.caseId;
        report.sourceElementType = std::move(hdf5.sourceElementType);
        report.internalFormulation = std::move(hdf5.internalFormulation);
        report.integrationRule = std::move(hdf5.integrationRule);
        report.passed = true;
        report.rows.reserve(hdf5.rows.size() * kComponents.size());
        report.vectorMetrics.reserve(hdf5.rows.size());
        std::array<double, 6> errorNorms{};
        std::array<double, 6> maximumErrors{};
        std::array<double, 6> maximumNormalized{};
        std::array<std::size_t, 6> worstRows{};
        double globalWorstNormalized = -1.0;

        for (const auto& hdf5Row : hdf5.rows) {
            const auto reference = referenceByIdentity.find(hdf5Row.identity);
            if (reference == referenceByIdentity.end()) {
                fail("schema-mismatch", "A projected reference row is missing.");
            }
            std::array<double, 6> errors{};
            for (std::size_t component = 0U;
                 component < kComponents.size();
                 ++component) {
                const double tolerance = kFixedAbsoluteTolerance;
                const double absoluteError = std::abs(
                    hdf5Row.values[component] -
                    reference->second->values[component]);
                const double normalizedError = absoluteError / tolerance;
                if (!std::isfinite(tolerance) || !(tolerance > 0.0) ||
                    !std::isfinite(absoluteError) ||
                    !std::isfinite(normalizedError)) {
                    fail(
                        "schema-mismatch",
                        "A finite row produced a nonfinite comparison metric.");
                }
                const bool blocking = component < 3U;
                const bool withinTolerance = absoluteError <= tolerance;
                const std::size_t rowIndex = report.rows.size();
                report.rows.push_back({
                    referenceCase.caseId,
                    hdf5Row.identity.instanceName,
                    hdf5Row.identity.sourceNodeLabel,
                    kComponents[component],
                    hdf5Row.values[component],
                    reference->second->values[component],
                    absoluteError,
                    tolerance,
                    normalizedError,
                    blocking,
                    withinTolerance});
                errors[component] = absoluteError;
                errorNorms[component] =
                    std::hypot(errorNorms[component], absoluteError);
                if (normalizedError > maximumNormalized[component]) {
                    maximumNormalized[component] = normalizedError;
                    maximumErrors[component] = absoluteError;
                    worstRows[component] = rowIndex;
                }
                if (normalizedError > globalWorstNormalized) {
                    globalWorstNormalized = normalizedError;
                    report.worstRow = rowIndex;
                }
                if (blocking && !withinTolerance) {
                    report.passed = false;
                } else if (!blocking && !withinTolerance) {
                    const std::string message =
                        "case=" + referenceCase.caseId + ";instance=" +
                        hdf5Row.identity.instanceName + ";node=" +
                        std::to_string(hdf5Row.identity.sourceNodeLabel) +
                        ";component=" + kComponents[component] +
                        ";absolute_error=" + finiteText(absoluteError) +
                        ";tolerance=" + finiteText(tolerance);
                    report.warnings.push_back({
                        "rotation-reference-exceedance", rowIndex, message});
                }
            }
            report.vectorMetrics.push_back({
                hdf5Row.identity.instanceName,
                hdf5Row.identity.sourceNodeLabel,
                std::hypot(errors[0U], errors[1U], errors[2U]),
                std::hypot(errors[3U], errors[4U], errors[5U])});
        }

        report.metrics.reserve(kComponents.size());
        const double rowCount = static_cast<double>(hdf5.rows.size());
        for (std::size_t component = 0U;
             component < kComponents.size();
             ++component) {
            report.metrics.push_back({
                kComponents[component],
                scales[component],
                kFixedAbsoluteTolerance,
                maximumErrors[component],
                maximumNormalized[component],
                errorNorms[component] / std::sqrt(rowCount),
                errorNorms[component],
                worstRows[component]});
        }
        return Result<Mitc4ComparisonReport>::success(std::move(report));
    } catch (const ComparisonFailure& exception) {
        return Result<Mitc4ComparisonReport>::failure(failureStatus(
            referenceCase.caseId, exception.code(), exception.what()));
    } catch (const std::exception& exception) {
        return Result<Mitc4ComparisonReport>::failure(failureStatus(
            referenceCase.caseId, "comparison-failure", exception.what()));
    }
}

Status Mitc4ReferenceComparison::writeDeterministicJson(
    const Mitc4ComparisonReport& report,
    const std::filesystem::path& outputJson) {
    try {
        std::ofstream stream{outputJson, std::ios::binary | std::ios::trunc};
        if (!stream) {
            return failureStatus(
                report.caseId,
                "comparison-report-write-failed",
                "The deterministic MITC4 JSON report cannot be opened.");
        }
        stream.imbue(std::locale::classic());
        stream << std::setprecision(std::numeric_limits<double>::max_digits10);
        stream << "{\"case_id\":\"" << jsonEscape(report.caseId)
               << "\",\"source_element_type\":\""
               << jsonEscape(report.sourceElementType)
               << "\",\"internal_formulation\":\""
               << jsonEscape(report.internalFormulation)
               << "\",\"integration_rule\":\""
               << jsonEscape(report.integrationRule) << "\",\"rows\":[";
        for (std::size_t index = 0U; index < report.rows.size(); ++index) {
            if (index != 0U) {
                stream << ',';
            }
            const auto& row = report.rows[index];
            stream << "{\"case_id\":\"" << jsonEscape(row.caseId)
                   << "\",\"instance_name\":\""
                   << jsonEscape(row.instanceName)
                   << "\",\"source_node_label\":" << row.sourceNodeLabel
                   << ",\"component\":\"" << jsonEscape(row.component)
                   << "\",\"fesa_value\":" << row.fesaValue
                   << ",\"reference_value\":" << row.referenceValue
                   << ",\"absolute_error\":" << row.absoluteError
                   << ",\"tolerance\":" << row.tolerance
                   << ",\"normalized_error\":" << row.normalizedError
                   << ",\"blocking\":" << (row.blocking ? "true" : "false")
                   << ",\"within_tolerance\":"
                   << (row.withinTolerance ? "true" : "false") << '}';
        }
        stream << "],\"metrics\":[";
        for (std::size_t index = 0U; index < report.metrics.size(); ++index) {
            if (index != 0U) {
                stream << ',';
            }
            const auto& metric = report.metrics[index];
            stream << "{\"component\":\"" << jsonEscape(metric.component)
                   << "\",\"reference_scale\":" << metric.referenceScale
                   << ",\"tolerance\":" << metric.tolerance
                   << ",\"maximum_absolute_error\":"
                   << metric.maximumAbsoluteError
                   << ",\"maximum_normalized_error\":"
                   << metric.maximumNormalizedError
                   << ",\"rms_error\":" << metric.rmsError
                   << ",\"vector_norm_error\":" << metric.vectorNormError
                   << ",\"worst_row\":" << metric.worstRow << '}';
        }
        stream << "],\"vector_metrics\":[";
        for (std::size_t index = 0U; index < report.vectorMetrics.size(); ++index) {
            if (index != 0U) {
                stream << ',';
            }
            const auto& metric = report.vectorMetrics[index];
            stream << "{\"instance_name\":\""
                   << jsonEscape(metric.instanceName)
                   << "\",\"source_node_label\":" << metric.sourceNodeLabel
                   << ",\"displacement_norm_error\":"
                   << metric.displacementNormError
                   << ",\"rotation_norm_error\":"
                   << metric.rotationNormError << '}';
        }
        stream << "],\"warnings\":[";
        for (std::size_t index = 0U; index < report.warnings.size(); ++index) {
            if (index != 0U) {
                stream << ',';
            }
            const auto& warning = report.warnings[index];
            stream << "{\"code\":\"" << jsonEscape(warning.code)
                   << "\",\"row\":" << warning.row
                   << ",\"message\":\"" << jsonEscape(warning.message)
                   << "\"}";
        }
        stream << "],\"worst_row\":" << report.worstRow
               << ",\"passed\":" << (report.passed ? "true" : "false")
               << "}\n";
        stream.flush();
        if (!stream) {
            return failureStatus(
                report.caseId,
                "comparison-report-write-failed",
                "The deterministic MITC4 JSON report could not be completed.");
        }
        return Status::ok();
    } catch (const std::exception& exception) {
        return failureStatus(
            report.caseId, "comparison-report-write-failed", exception.what());
    }
}

} // namespace fesa::test