#include <fesa/validation/comparison.hpp>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <limits>
#include <map>
#include <optional>
#include <set>
#include <sstream>
#include <string>
#include <string_view>
#include <tuple>
#include <utility>
#include <vector>

namespace fesa {
namespace {

using PositionKey = std::tuple<
    ReferenceQuantity,
    std::string,
    std::int64_t,
    std::optional<std::int64_t>>;

void add_failure(
    std::vector<Diagnostic>& failures,
    std::string code,
    std::string message) {
    failures.push_back({
        DiagnosticStage::validation,
        Severity::error,
        std::move(code),
        std::move(message),
        std::nullopt,
    });
}

std::string_view quantity_name(const ReferenceQuantity quantity) {
    switch (quantity) {
        case ReferenceQuantity::displacement:
            return "displacement";
        case ReferenceQuantity::reaction:
            return "reaction";
        case ReferenceQuantity::internal_force:
            return "internal_force";
        case ReferenceQuantity::centroid_stress:
            return "centroid_stress";
    }
    return "unknown";
}

std::size_t expected_component_count(
    const ReferenceQuantity quantity) {
    switch (quantity) {
        case ReferenceQuantity::displacement:
        case ReferenceQuantity::reaction:
        case ReferenceQuantity::internal_force:
            return 6;
        case ReferenceQuantity::centroid_stress:
            return 1;
    }
    return 0;
}

std::string component_name(
    const ReferenceQuantity quantity,
    const std::size_t index) {
    switch (quantity) {
        case ReferenceQuantity::displacement:
            return std::string{NodalFrame::displacement_components[index]};
        case ReferenceQuantity::reaction:
            return std::string{NodalFrame::reaction_components[index]};
        case ReferenceQuantity::internal_force:
            return std::string{
                BeamElementFrame::section_force_components[index]};
        case ReferenceQuantity::centroid_stress:
            return std::string{BeamElementFrame::axial_stress_component};
    }
    return "unknown";
}

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

std::string tolerance_text(const Tolerance tolerance) {
    return " relative_tolerance=" + number_text(tolerance.relative) +
           " absolute_scale=" + number_text(tolerance.absolute_scale);
}

std::string position_text(
    const ReferenceQuantity quantity,
    const ResultPosition& position) {
    std::string text = "quantity=" + std::string{quantity_name(quantity)} +
                       " instance=" + position.instance_name +
                       " entity=" + std::to_string(position.entity_label);
    if (position.end_node_label.has_value()) {
        text += " end_node=" +
                std::to_string(*position.end_node_label);
    } else {
        text += " end_node=n/a";
    }
    return text;
}

std::string scalar_failure_text(
    const ComparisonSample& sample,
    const std::size_t component,
    const double normalized_error) {
    return position_text(sample.quantity, sample.position) +
           " component=" + component_name(sample.quantity, component) +
           " reference=" + number_text(sample.reference[component]) +
           " actual=" + number_text(sample.actual[component]) +
           " normalized_error=" + number_text(normalized_error) +
           tolerance_text(sample.tolerance);
}

std::string unevaluable_failure_text(
    const ReferenceQuantity quantity,
    const ResultPosition& position,
    const Tolerance tolerance,
    const std::string_view reason) {
    return position_text(quantity, position) +
           " component=n/a reference=n/a actual=n/a "
           "normalized_error=inf" + tolerance_text(tolerance) +
           " reason=" + std::string{reason};
}

bool origin_matches(
    const EntityOrigin& origin,
    const std::string& instance_name,
    const std::int64_t local_label) {
    return origin.instance_name == instance_name &&
           origin.local_label == local_label;
}

ComparisonSampleMatch matching_failure(
    const ReferenceQuantity quantity,
    const ResultPosition& position,
    const Tolerance tolerance,
    std::string code,
    const std::string_view reason) {
    std::vector<Diagnostic> failures;
    add_failure(
        failures,
        std::move(code),
        unevaluable_failure_text(
            quantity, position, tolerance, reason));
    return {std::nullopt, std::move(failures)};
}

std::vector<double> as_vector(const std::array<double, 6>& values) {
    return {values.begin(), values.end()};
}

}  // namespace

ComparisonSampleMatch make_comparison_sample(
    const ResultFrame& frame,
    const ReferenceQuantity quantity,
    const ResultPosition& position,
    const std::span<const double> reference,
    const Tolerance tolerance) {
    std::vector<double> actual;

    if (
        quantity == ReferenceQuantity::displacement ||
        quantity == ReferenceQuantity::reaction) {
        if (position.end_node_label.has_value()) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.invalid_result_position",
                "nodal_position_has_end_node");
        }

        std::optional<std::size_t> matched_index;
        for (
            std::size_t index = 0;
            index < frame.nodal.origins.size();
            ++index) {
            if (origin_matches(
                    frame.nodal.origins[index],
                    position.instance_name,
                    position.entity_label)) {
                if (matched_index.has_value()) {
                    return matching_failure(
                        quantity,
                        position,
                        tolerance,
                        "validation.unknown_result_origin",
                        "ambiguous_result_origin");
                }
                matched_index = index;
            }
        }
        if (!matched_index.has_value()) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.unknown_result_origin",
                "unknown_result_origin");
        }

        const auto& field =
            quantity == ReferenceQuantity::displacement
                ? frame.nodal.displacement
                : frame.nodal.reaction;
        if (*matched_index >= field.size()) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.component_count_mismatch",
                "missing_actual_components");
        }
        actual = as_vector(field[*matched_index]);
    } else {
        if (!position.end_node_label.has_value()) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.invalid_element_node_pair",
                "missing_end_node");
        }

        const BeamElementFrame* matched_beam = nullptr;
        for (const BeamElementFrame& beam : frame.element.beams) {
            if (origin_matches(
                    beam.origin,
                    position.instance_name,
                    position.entity_label)) {
                if (matched_beam != nullptr) {
                    return matching_failure(
                        quantity,
                        position,
                        tolerance,
                        "validation.unknown_result_origin",
                        "ambiguous_result_origin");
                }
                matched_beam = &beam;
            }
        }
        if (matched_beam == nullptr) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.unknown_result_origin",
                "unknown_result_origin");
        }

        std::optional<NodeId> end_node;
        for (
            std::size_t index = 0;
            index < frame.nodal.origins.size() &&
            index < frame.nodal.node_ids.size();
            ++index) {
            if (origin_matches(
                    frame.nodal.origins[index],
                    position.instance_name,
                    *position.end_node_label)) {
                if (end_node.has_value()) {
                    return matching_failure(
                        quantity,
                        position,
                        tolerance,
                        "validation.unknown_result_origin",
                        "ambiguous_end_node_origin");
                }
                end_node = frame.nodal.node_ids[index];
            }
        }
        if (!end_node.has_value()) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.unknown_result_origin",
                "unknown_end_node_origin");
        }

        const BeamSectionResult* matched_end = nullptr;
        for (const BeamSectionResult& end : matched_beam->end_results) {
            if (end.end_node == *end_node) {
                matched_end = &end;
                break;
            }
        }
        if (matched_end == nullptr) {
            return matching_failure(
                quantity,
                position,
                tolerance,
                "validation.invalid_element_node_pair",
                "node_is_not_element_end");
        }

        if (quantity == ReferenceQuantity::internal_force) {
            actual = as_vector(matched_end->section_force);
        } else {
            actual = {matched_end->centroid_sigma_xx};
        }
    }

    ComparisonSample matched{
        quantity,
        position,
        {reference.begin(), reference.end()},
        std::move(actual),
        tolerance,
    };
    return {std::move(matched), {}};
}

ComparisonReport compare_samples(
    const std::span<const ComparisonSample> samples) {
    ComparisonReport report{true, 0.0, {}};
    std::set<PositionKey> positions;

    for (const ComparisonSample& sample : samples) {
        const PositionKey key{
            sample.quantity,
            sample.position.instance_name,
            sample.position.entity_label,
            sample.position.end_node_label,
        };
        if (!positions.insert(key).second) {
            add_failure(
                report.failures,
                "validation.duplicate_result_position",
                unevaluable_failure_text(
                    sample.quantity,
                    sample.position,
                    sample.tolerance,
                    "duplicate_result_position"));
            report.maximum_normalized_error =
                std::numeric_limits<double>::infinity();
            continue;
        }

        const std::size_t expected =
            expected_component_count(sample.quantity);
        if (
            sample.reference.size() != expected ||
            sample.actual.size() != expected) {
            add_failure(
                report.failures,
                "validation.component_count_mismatch",
                unevaluable_failure_text(
                    sample.quantity,
                    sample.position,
                    sample.tolerance,
                    "component_count_mismatch") +
                    " expected=" + std::to_string(expected) +
                    " reference_count=" +
                    std::to_string(sample.reference.size()) +
                    " actual_count=" +
                    std::to_string(sample.actual.size()));
            report.maximum_normalized_error =
                std::numeric_limits<double>::infinity();
            continue;
        }

        if (
            !std::isfinite(sample.tolerance.relative) ||
            !std::isfinite(sample.tolerance.absolute_scale) ||
            sample.tolerance.relative < 0.0 ||
            sample.tolerance.absolute_scale < 0.0) {
            add_failure(
                report.failures,
                "validation.invalid_tolerance",
                unevaluable_failure_text(
                    sample.quantity,
                    sample.position,
                    sample.tolerance,
                    "invalid_tolerance"));
            report.maximum_normalized_error =
                std::numeric_limits<double>::infinity();
            continue;
        }

        for (std::size_t component = 0; component < expected; ++component) {
            const double reference = sample.reference[component];
            const double actual = sample.actual[component];
            if (!std::isfinite(reference) || !std::isfinite(actual)) {
                add_failure(
                    report.failures,
                    "validation.nonfinite_comparison_value",
                    scalar_failure_text(
                        sample,
                        component,
                        std::numeric_limits<double>::infinity()));
                report.maximum_normalized_error =
                    std::numeric_limits<double>::infinity();
                continue;
            }

            const double denominator =
                sample.tolerance.absolute_scale +
                sample.tolerance.relative * std::abs(reference);
            if (!(denominator > 0.0) || !std::isfinite(denominator)) {
                add_failure(
                    report.failures,
                    "validation.invalid_tolerance",
                    scalar_failure_text(
                        sample,
                        component,
                        std::numeric_limits<double>::infinity()));
                report.maximum_normalized_error =
                    std::numeric_limits<double>::infinity();
                continue;
            }

            const double normalized_error =
                std::abs(actual - reference) / denominator;
            report.maximum_normalized_error = std::max(
                report.maximum_normalized_error,
                normalized_error);
            if (!(normalized_error <= 1.0)) {
                add_failure(
                    report.failures,
                    "validation.tolerance_exceeded",
                    scalar_failure_text(
                        sample, component, normalized_error));
            }
        }
    }

    report.passed = report.failures.empty();
    return report;
}

CorrelationReport correlate_samples(
    const std::span<const ComparisonSample> samples) {
    struct Accumulator final {
        double squared_error{};
        double squared_reference{};
        double squared_absolute_scale{};
        std::size_t value_count{};
    };

    CorrelationReport report{true, {}, {}};
    if (samples.empty()) {
        add_failure(
            report.failures,
            "validation.empty_comparison",
            "No comparison samples were provided for correlation.");
        report.evaluable = false;
        return report;
    }

    std::set<PositionKey> positions;
    std::map<std::pair<ReferenceQuantity, std::size_t>, Accumulator>
        accumulators;

    for (const ComparisonSample& sample : samples) {
        const PositionKey key{
            sample.quantity,
            sample.position.instance_name,
            sample.position.entity_label,
            sample.position.end_node_label,
        };
        if (!positions.insert(key).second) {
            add_failure(
                report.failures,
                "validation.duplicate_result_position",
                unevaluable_failure_text(
                    sample.quantity,
                    sample.position,
                    sample.tolerance,
                    "duplicate_result_position"));
            continue;
        }

        const std::size_t expected =
            expected_component_count(sample.quantity);
        if (
            sample.reference.size() != expected ||
            sample.actual.size() != expected) {
            add_failure(
                report.failures,
                "validation.component_count_mismatch",
                unevaluable_failure_text(
                    sample.quantity,
                    sample.position,
                    sample.tolerance,
                    "component_count_mismatch") +
                    " expected=" + std::to_string(expected) +
                    " reference_count=" +
                    std::to_string(sample.reference.size()) +
                    " actual_count=" +
                    std::to_string(sample.actual.size()));
            continue;
        }

        if (
            !std::isfinite(sample.tolerance.relative) ||
            !std::isfinite(sample.tolerance.absolute_scale) ||
            sample.tolerance.relative < 0.0 ||
            sample.tolerance.absolute_scale < 0.0) {
            add_failure(
                report.failures,
                "validation.invalid_tolerance",
                unevaluable_failure_text(
                    sample.quantity,
                    sample.position,
                    sample.tolerance,
                    "invalid_tolerance"));
            continue;
        }

        for (std::size_t component = 0; component < expected; ++component) {
            const double reference = sample.reference[component];
            const double actual = sample.actual[component];
            if (!std::isfinite(reference) || !std::isfinite(actual)) {
                add_failure(
                    report.failures,
                    "validation.nonfinite_comparison_value",
                    scalar_failure_text(
                        sample,
                        component,
                        std::numeric_limits<double>::infinity()));
                continue;
            }

            const double error = actual - reference;
            Accumulator& accumulator =
                accumulators[{sample.quantity, component}];
            accumulator.squared_error += error * error;
            accumulator.squared_reference += reference * reference;
            accumulator.squared_absolute_scale +=
                sample.tolerance.absolute_scale *
                sample.tolerance.absolute_scale;
            ++accumulator.value_count;
        }
    }

    for (const auto& [key, accumulator] : accumulators) {
        const auto [quantity, component] = key;
        const double error_l2 = std::sqrt(accumulator.squared_error);
        const double root_mean_square_error =
            error_l2 /
            std::sqrt(static_cast<double>(accumulator.value_count));
        const double reference_l2 =
            std::sqrt(accumulator.squared_reference);
        const double absolute_scale_l2 =
            std::sqrt(accumulator.squared_absolute_scale);
        const double denominator =
            std::max(reference_l2, absolute_scale_l2);
        const double relative_l2_error = denominator > 0.0
                                             ? error_l2 / denominator
                                             : error_l2 == 0.0
                                                   ? 0.0
                                                   : std::numeric_limits<
                                                         double>::infinity();

        report.metrics.push_back({
            quantity,
            component,
            accumulator.value_count,
            root_mean_square_error,
            relative_l2_error,
        });
        if (
            !std::isfinite(root_mean_square_error) ||
            !std::isfinite(relative_l2_error)) {
            add_failure(
                report.failures,
                "validation.nonfinite_correlation_metric",
                "quantity=" + std::string{quantity_name(quantity)} +
                    " component=" + component_name(quantity, component) +
                    " rmse=" + number_text(root_mean_square_error) +
                    " relative_l2=" + number_text(relative_l2_error));
        }
    }

    report.evaluable = report.failures.empty();
    return report;
}

}  // namespace fesa