// Las pruebas de la semana 15: comprueban lo que dicen los apuntes. Cada grupo es una función.
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <iostream>
#include <numbers>
#include <numeric>
#include <random>
#include <string>
#include <utility>
#include <vector>

#include "../week-10/held_karp.h"
#include "bound.h"
#include "cover.h"
#include "geo.h"
#include "tour.h"

namespace {

using geo::at;
using tour::Matrix;
using tour::Tour;

// Una instancia métrica al azar: puntos en un plano, distancia euclidiana redondeada hacia arriba.
Matrix random_metric(int n, std::mt19937& rng) {
    std::vector<std::pair<int, int>> p;
    for (int i = 0; i < n; i++) p.push_back({static_cast<int>(rng() % 1000), static_cast<int>(rng() % 1000)});
    Matrix d(at(n), std::vector<long long>(at(n), 0));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            const double dx = p[at(i)].first - p[at(j)].first;
            const double dy = p[at(i)].second - p[at(j)].second;
            d[at(i)][at(j)] = static_cast<long long>(std::ceil(std::sqrt(dx * dx + dy * dy)));
        }
    }
    return d;
}

// Pesos simétricos sin ninguna relación entre sí: no cumplen la desigualdad triangular.
Matrix random_symmetric(int n, std::mt19937& rng) {
    Matrix d(at(n), std::vector<long long>(at(n), 0));
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) d[at(i)][at(j)] = d[at(j)][at(i)] = 1 + static_cast<long long>(rng() % 1000);
    }
    return d;
}

// El árbitro del árbol de expansión mínima: Kruskal con union-find.
long long kruskal_weight(const Matrix& d) {
    const int n = static_cast<int>(d.size());
    std::vector<std::pair<long long, std::pair<int, int>>> edges;
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) edges.push_back({d[at(i)][at(j)], {i, j}});
    }
    std::sort(edges.begin(), edges.end());
    std::vector<int> parent(at(n));
    std::iota(parent.begin(), parent.end(), 0);
    auto find = [&](int v) {
        while (parent[at(v)] != v) v = parent[at(v)] = parent[at(parent[at(v)])];
        return v;
    };
    long long total = 0;
    for (const auto& [w, e] : edges) {
        const int a = find(e.first);
        const int b = find(e.second);
        if (a != b) {
            parent[at(a)] = b;
            total += w;
        }
    }
    return total;
}

bool triangle_holds(const Matrix& d) {
    const std::size_t n = d.size();
    for (std::size_t a = 0; a < n; a++) {
        for (std::size_t b = 0; b < n; b++) {
            for (std::size_t c = 0; c < n; c++) {
                if (d[a][c] > d[a][b] + d[b][c]) return false;
            }
        }
    }
    return true;
}

void test_distances() {
    const auto all = geo::read_places("capitales.csv");
    assert(all.size() == 26);
    const auto d = geo::distance_matrix(all);
    assert(triangle_holds(d));
    for (std::size_t i = 0; i < all.size(); i++) {
        assert(d[i][i] == 0);
        for (std::size_t j = 0; j < all.size(); j++) {
            assert(d[i][j] == d[j][i]);
            assert(i == j || d[i][j] > 0);
        }
    }
    // Lima y Callao son el par más cercano; Tacna y Tumbes, el más lejano.
    const auto index_of = [&](const std::string& name) {
        for (int i = 0; i < 26; i++) {
            if (all[at(i)].name == name) return i;
        }
        return -1;
    };
    const int lima = index_of("Lima");
    const int callao = index_of("Callao");
    assert(d[at(lima)][at(callao)] == 11596);
    assert(d[at(index_of("Tacna"))][at(index_of("Tumbes"))] == 1953472);
    for (std::size_t i = 0; i < all.size(); i++) {
        for (std::size_t j = i + 1; j < all.size(); j++) assert(d[i][j] >= 11596 && d[i][j] <= 1953472);
    }
    // Redondear hacia arriba mantiene la desigualdad en puntos cualesquiera; al más cercano, no siempre.
    const double degrees = 1000.4 / (geo::EARTH_RADIUS_M * std::numbers::pi / 180.0);
    const geo::Place a{"a", 0, 0};
    const geo::Place b{"b", 0, degrees};
    const geo::Place c{"c", 0, 2 * degrees};
    assert(std::llround(geo::great_circle_m(a, b)) + std::llround(geo::great_circle_m(b, c)) <
           std::llround(geo::great_circle_m(a, c)));
    assert(geo::distance_m(a, b) + geo::distance_m(b, c) >= geo::distance_m(a, c));
    std::mt19937 rng(1);
    for (int round = 0; round < 3000; round++) {
        std::vector<geo::Place> p;
        for (int i = 0; i < 3; i++) {
            p.push_back({"", -90.0 + rng() % 18001 / 100.0, -180.0 + rng() % 36001 / 100.0});
        }
        assert(geo::distance_m(p[0], p[2]) <= geo::distance_m(p[0], p[1]) + geo::distance_m(p[1], p[2]));
    }
}

void test_spanning_tree() {
    std::mt19937 rng(2);
    for (int round = 0; round < 300; round++) {
        const int n = 1 + static_cast<int>(rng() % 25);
        const auto d = round % 2 == 0 ? random_metric(n, rng) : random_symmetric(n, rng);
        assert(tour::minimum_spanning_tree(d).weight == kruskal_weight(d));
    }
}

// Lo que dice la demostración del factor 2, paso por paso.
void check_double_tree(const Matrix& d, long long optimum) {
    const int n = static_cast<int>(d.size());
    const auto tree = tour::minimum_spanning_tree(d);
    const auto walk = tour::euler_walk(tree);
    assert(walk.size() == at(2 * n - 1));
    assert(tour::walk_length(d, walk) == 2 * tree.weight);
    const auto t = tour::shortcut(walk, n);
    assert(tour::is_tour(t, n));
    assert(tree.weight <= optimum);
    assert(tour::tour_length(d, t) <= 2 * tree.weight);
    assert(tour::tour_length(d, t) <= 2 * optimum);
    // 2-opt desde ahí nunca empeora, y mantiene el factor por haber empezado en el árbol.
    auto improved = t;
    const auto moves = tour::two_opt_plain(d, improved);
    assert(tour::is_tour(improved, n));
    assert(tour::is_two_opt_local_optimum(d, improved));
    assert(optimum <= tour::tour_length(d, improved));
    assert(tour::tour_length(d, improved) <= tour::tour_length(d, t));
    long long previous = tour::tour_length(d, t);
    for (const auto& move : moves) {
        assert(move.length_after < previous);
        previous = move.length_after;
    }
}

void test_factor_two() {
    const auto all = geo::read_places("capitales.csv");
    for (int n = 3; n <= 14; n++) {
        const auto d = geo::distance_matrix(geo::first_places(all, n));
        check_double_tree(d, tsp::held_karp(d));
    }
    std::mt19937 rng(3);
    for (int round = 0; round < 300; round++) {
        const int n = 3 + static_cast<int>(rng() % 9);
        const auto d = random_metric(n, rng);
        assert(triangle_holds(d));
        check_double_tree(d, tsp::held_karp(d));
    }
}

void test_no_guarantee_without_triangle() {
    // El ciclo desordenado del ejemplo: el árbol duplicado da 30 y el óptimo es 6, cinco veces más.
    const int n = 6;
    const long long heavy = 2 * n + 1;
    Matrix w(at(n), std::vector<long long>(at(n), heavy));
    for (int v = 0; v < n; v++) w[at(v)][at(v)] = 0;
    const std::vector<std::pair<int, int>> edges = {{0, 3}, {3, 1}, {1, 5}, {5, 2}, {2, 4}, {4, 0}};
    for (const auto& [u, v] : edges) w[at(u)][at(v)] = w[at(v)][at(u)] = 1;
    assert(!triangle_holds(w));
    assert(tsp::held_karp(w) == 6);
    assert(tour::tour_length(w, tour::double_tree(w)) == 30);
}

// ¿Hay ciclo hamiltoniano? Árbitro: prueba todos los órdenes.
bool has_hamiltonian_cycle(int n, const std::vector<std::vector<char>>& adjacent) {
    std::vector<int> order(at(n));
    std::iota(order.begin(), order.end(), 0);
    do {
        bool ok = true;
        for (int i = 0; i < n && ok; i++) ok = adjacent[at(order[at(i)])][at(order[at((i + 1) % n)])] != 0;
        if (ok) return true;
    } while (std::next_permutation(order.begin() + 1, order.end()));
    return false;
}

void test_hamilton_reduction() {
    std::mt19937 rng(4);
    int with_cycle = 0;
    int without_cycle = 0;
    for (int round = 0; round < 600; round++) {
        const int n = 3 + static_cast<int>(rng() % 5);
        const long long alpha = 1 + static_cast<long long>(rng() % 4);
        const long long heavy = alpha * n + 1;
        std::vector<std::vector<char>> adjacent(at(n), std::vector<char>(at(n), 0));
        Matrix w(at(n), std::vector<long long>(at(n), heavy));
        for (int u = 0; u < n; u++) {
            w[at(u)][at(u)] = 0;
            for (int v = u + 1; v < n; v++) {
                if (rng() % 2 == 0) {
                    adjacent[at(u)][at(v)] = adjacent[at(v)][at(u)] = 1;
                    w[at(u)][at(v)] = w[at(v)][at(u)] = 1;
                }
            }
        }
        const long long optimum = tsp::held_karp(w);
        if (has_hamiltonian_cycle(n, adjacent)) {
            with_cycle++;
            assert(optimum == n);
        } else {
            without_cycle++;
            assert(optimum >= heavy + n - 1);
            assert(optimum > alpha * n);
        }
    }
    assert(with_cycle > 50 && without_cycle > 50);
}

void test_two_opt_depends_on_start() {
    // Sin el árbol no hay garantía, y dos inicios pueden terminar en óptimos locales distintos.
    const auto d = geo::distance_matrix(geo::first_places(geo::read_places("capitales.csv"), 20));
    const long long optimum = tsp::held_karp(d);
    auto from_tree = tour::double_tree(d);
    tour::two_opt_plain(d, from_tree);
    assert(tour::tour_length(d, from_tree) == 4867292);
    assert(optimum == 4640216);
    auto from_neighbor = tour::nearest_neighbor(d);
    tour::two_opt_plain(d, from_neighbor);
    assert(tour::tour_length(d, from_neighbor) == 4827795);
}

void test_tour_helpers() {
    assert(tour::is_tour({2, 0, 1}, 3));
    assert(!tour::is_tour({2, 0, 0}, 3));
    assert(!tour::is_tour({2, 0}, 3));
    assert(!tour::is_tour({2, 0, 3}, 3));
}

// Lo que completion_bound promete: nunca pasa de la mejor forma de terminar.
void test_bound_is_sound() {
    std::mt19937 rng(5);
    for (int round = 0; round < 200; round++) {
        const int n = 3 + static_cast<int>(rng() % 6);
        const auto d = round % 2 == 0 ? random_metric(n, rng) : random_symmetric(n, rng);
        const int last = 1 + static_cast<int>(rng() % static_cast<unsigned>(n - 1));
        std::vector<int> rest;
        for (int v = 1; v < n; v++) {
            if (v != last && rng() % 2 == 0) rest.push_back(v);
        }
        std::sort(rest.begin(), rest.end());
        long long best = bound::NO_TOUR;
        auto order = rest;
        do {
            long long total = 0;
            int current = last;
            for (int v : order) {
                total += d[at(current)][at(v)];
                current = v;
            }
            best = std::min(best, total + d[at(current)][0]);
        } while (std::next_permutation(order.begin(), order.end()));
        assert(bound::completion_bound(d, last, rest) <= best);
    }
}

void test_search_is_exact() {
    std::mt19937 rng(6);
    for (int round = 0; round < 200; round++) {
        const int n = 1 + static_cast<int>(rng() % 9);
        const auto d = round % 2 == 0 ? random_metric(n, rng) : random_symmetric(n, rng);
        const long long optimum = tsp::held_karp(d);
        const auto none = bound::solve(d, bound::Pruning::none);
        const auto partial = bound::solve(d, bound::Pruning::partial);
        const auto spanning = bound::solve(d, bound::Pruning::spanning);
        assert(none.cost == optimum && partial.cost == optimum && spanning.cost == optimum);
        assert(partial.nodes <= none.nodes && spanning.nodes <= none.nodes);
    }
    const auto all = geo::read_places("capitales.csv");
    // Las cuentas de la tabla de los apuntes. El árbol sin podar tiene una fórmula.
    const auto d10 = geo::distance_matrix(geo::first_places(all, 10));
    long long expected = 0;
    long long level = 1;
    for (int k = 0; k <= 9; k++) {
        expected += level;
        level *= 9 - k;
    }
    assert(bound::solve(d10, bound::Pruning::none).nodes == expected);
    assert(expected == 986410);
    assert(bound::solve(d10, bound::Pruning::partial).nodes == 146810);
    assert(bound::solve(d10, bound::Pruning::spanning).nodes == 3157);
    assert(bound::solve(d10, bound::Pruning::spanning).cost == 3167199);
    // El corte por cantidad de nodos avisa que el costo no vale.
    assert(!bound::solve(d10, bound::Pruning::none, 100).finished);
    const auto d15 = geo::distance_matrix(geo::first_places(all, 15));
    assert(bound::solve(d15, bound::Pruning::spanning).cost == tsp::held_karp(d15));
}

void test_cover() {
    std::mt19937 rng(7);
    for (int round = 0; round < 1500; round++) {
        const int n = 1 + static_cast<int>(rng() % 11);
        cover::Edges edges;
        for (int u = 0; u < n; u++) {
            for (int v = u + 1; v < n; v++) {
                if (rng() % 3 == 0) edges.push_back({u, v});
            }
        }
        std::shuffle(edges.begin(), edges.end(), rng);
        const int optimum = cover::minimum_cover_size(n, edges);
        const auto matched = cover::matching_cover(n, edges);
        const int picked = static_cast<int>(matched.cover.size());
        assert(cover::is_cover(matched.cover, n, edges));
        assert(picked == 2 * static_cast<int>(matched.matching.size()));
        // Aristas sin nodos en común: el óptimo necesita un nodo distinto por cada una.
        std::vector<char> used(at(n), 0);
        for (const auto& [u, v] : matched.matching) {
            assert(used[at(u)] == 0 && used[at(v)] == 0);
            used[at(u)] = used[at(v)] = 1;
        }
        // Maximal: toda arista toca al matching.
        for (const auto& [u, v] : edges) assert(used[at(u)] != 0 || used[at(v)] != 0);
        assert(static_cast<int>(matched.matching.size()) <= optimum);
        assert(picked <= 2 * optimum);
        assert(cover::is_cover(cover::greedy_degree_cover(n, edges), n, edges));
    }
    // La estrella: una sola punta elige todas las hojas. Una arista: el factor 2 se alcanza.
    const cover::Edges star = {{0, 1}, {0, 2}, {0, 3}, {0, 4}, {0, 5}, {0, 6}};
    assert(cover::one_endpoint_cover(7, star).size() == 6);
    assert(cover::minimum_cover_size(7, star) == 1);
    assert(cover::matching_cover(2, {{0, 1}}).cover.size() == 2);
    assert(cover::minimum_cover_size(2, {{0, 1}}) == 1);
    // Maximal sin ser máximo: en el camino 0-1-2-3, la arista 1-2 ya no deja meter otra.
    assert(cover::matching_cover(4, {{1, 2}, {0, 1}, {2, 3}}).matching.size() == 1);
    // La trampa: el mayor grado elige toda la derecha, y el matching demuestra que el óptimo es k.
    for (int k : {6, 12, 24, 48, 96, 192}) {
        const auto trap = cover::greedy_trap(k);
        int count = 0;
        for (int i = 2; i <= k; i++) count += k / i;
        assert(trap.right_nodes == count);
        const auto chosen = cover::greedy_degree_cover(trap.nodes, trap.edges);
        assert(static_cast<int>(chosen.size()) == count);
        for (int v : chosen) assert(v >= k);
        assert(static_cast<int>(cover::matching_cover(trap.nodes, trap.edges).matching.size()) == k);
    }
    assert(cover::greedy_degree_cover(84, cover::greedy_trap(24).edges).size() == 60);
}

// Programación lineal, solo con cuentas a mano: no hay un solucionador en el entorno.
void test_linear_programming_examples() {
    // máx 3x + 2y con x + y ≤ 4, x ≤ 3, y ≤ 3, x, y ≥ 0. Se mira cada esquina y una grilla fina.
    const auto feasible = [](double x, double y) { return x >= 0 && y >= 0 && x + y <= 4 && x <= 3 && y <= 3; };
    const std::vector<std::pair<double, double>> corners = {{0, 0}, {3, 0}, {3, 1}, {1, 3}, {0, 3}};
    double best = 0;
    for (const auto& [x, y] : corners) {
        assert(feasible(x, y));
        best = std::max(best, 3 * x + 2 * y);
    }
    assert(best == 11);
    for (int i = 0; i <= 400; i++) {
        for (int j = 0; j <= 400; j++) {
            const double x = i / 100.0;
            const double y = j / 100.0;
            if (feasible(x, y)) assert(3 * x + 2 * y <= best);
        }
    }
    // El triángulo: con ½ en cada nodo se cumple toda arista y la suma es 3/2. Con enteros, 2.
    const cover::Edges triangle = {{0, 1}, {1, 2}, {0, 2}};
    assert(0.5 + 0.5 >= 1 && 3 * 0.5 == 1.5);
    assert(cover::minimum_cover_size(3, triangle) == 2);
}

// Las cuentas de la pregunta 5 de repaso: 12 clientes.
void test_review_numbers() {
    long long transitions = 0;
    std::mt19937 rng(8);
    const auto d = random_metric(12, rng);
    static_cast<void>(tsp::held_karp(d, &transitions));
    assert(transitions == 56320);
    long long factorial = 1;
    for (int i = 2; i <= 11; i++) factorial *= i;
    assert(factorial == 39916800);
}

}  // namespace

int main() {
    test_distances();
    test_spanning_tree();
    test_factor_two();
    test_no_guarantee_without_triangle();
    test_hamilton_reduction();
    test_two_opt_depends_on_start();
    test_tour_helpers();
    test_bound_is_sound();
    test_search_is_exact();
    test_cover();
    test_linear_programming_examples();
    test_review_numbers();
    std::cout << "Todas las pruebas de la semana 15 pasaron.\n";
}