// Comprueba lo que afirman los apuntes. Se compila con sanitizers.
#include <algorithm>
#include <cassert>
#include <functional>
#include <iostream>
#include <random>
#include <vector>

#include "all_pairs.h"
#include "held_karp.h"
#include "rates.h"
#include "shortest_paths.h"

namespace {

std::vector<Edge> random_graph(int n, int percent, int low, int high, std::mt19937& gen) {
    std::vector<Edge> edges;
    for (int a = 0; a < n; a++) {
        for (int b = 0; b < n; b++) {
            if (a != b && static_cast<int>(gen() % 100) < percent) {
                edges.push_back({a, b, low + static_cast<long long>(gen() % static_cast<unsigned>(high - low + 1))});
            }
        }
    }
    std::shuffle(edges.begin(), edges.end(), gen);
    return edges;
}

// Árbitro: el menor largo de un camino de source a cada nodo con a lo sumo limit aristas, probándolos todos.
std::vector<long long> walks_up_to(int n, const std::vector<Edge>& edges, int source, int limit) {
    std::vector<long long> best(idx(n), INF);
    const std::function<void(int, long long, int)> go = [&](int u, long long length, int used) {
        best[idx(u)] = std::min(best[idx(u)], length);
        if (used == limit) return;
        for (const Edge& e : edges) {
            if (e.from == u) go(e.to, length + e.weight, used + 1);
        }
    };
    go(source, 0, 0);
    return best;
}

// Árbitro: recorre los ciclos simples entre los nodos permitidos y dice si alguno pesa menos que cero.
// Cada ciclo se ve una vez, desde su nodo menor.
bool has_negative_cycle(int n, const std::vector<Edge>& edges, const std::vector<bool>& allowed) {
    std::vector<bool> used(idx(n), false);
    bool found = false;
    const std::function<void(int, int, long long)> go = [&](int start, int u, long long length) {
        used[idx(u)] = true;
        for (const Edge& e : edges) {
            if (e.from != u || e.to < start || !allowed[idx(e.to)]) continue;
            if (e.to == start) {
                if (length + e.weight < 0) found = true;
            } else if (!used[idx(e.to)]) {
                go(start, e.to, length + e.weight);
            }
        }
        used[idx(u)] = false;
    };
    for (int s = 0; s < n; s++) {
        if (allowed[idx(s)]) go(s, s, 0);
    }
    return found;
}

std::vector<bool> reachable_from(int n, const std::vector<Edge>& edges, int source) {
    const std::vector<long long> hops = bfs_distances(n, edges, source);
    std::vector<bool> reached(idx(n));
    for (int v = 0; v < n; v++) reached[idx(v)] = hops[idx(v)] != INF;
    return reached;
}

void test_recurrence_and_one_row() {
    std::mt19937 gen(1);
    int with_negative = 0;
    for (int t = 0; t < 1500; t++) {
        const int n = 2 + static_cast<int>(gen() % 4);
        const std::vector<Edge> edges = random_graph(n, 45, -3, 8, gen);
        // La tabla es la definición, con cualquier signo y con ciclos negativos también.
        const auto table = bellman_ford_table(n, edges, 0, 4);
        for (int i = 0; i <= 4; i++) assert(table[idx(i)] == walks_up_to(n, edges, 0, i));

        const std::vector<bool> reached = reachable_from(n, edges, 0);
        const bool negative = has_negative_cycle(n, edges, reached);
        if (negative) with_negative++;
        const auto full = bellman_ford_table(n, edges, 0, n);
        std::vector<long long> previous;
        const BellmanFordResult r = bellman_ford(n, edges, 0, [&](int round, const BellmanFordResult& partial) {
            for (int v = 0; v < n; v++) {
                assert(partial.dist[idx(v)] <= full[idx(round)][idx(v)]);  // nunca por encima de la tabla
                if (!negative) assert(partial.dist[idx(v)] >= full[idx(n)][idx(v)]);  // ni por debajo de la verdad
                if (!previous.empty()) assert(partial.dist[idx(v)] <= previous[idx(v)]);  // solo bajan
            }
            previous = partial.dist;
        });
        if (!negative) {
            // El lema de V - 1 aristas: la tabla ya no cambia, y una fila llega a lo mismo.
            assert(full[idx(n - 1)] == full[idx(n)]);
            assert(r.dist == full[idx(n - 1)]);
            assert(r.witness == -1);
        } else {
            assert(r.witness != -1);
        }
        // d[v] >= d[pred[v]] + w, salvo en el origen sin relajar.
        for (int v = 0; v < n; v++) {
            if (r.pred[idx(v)] != -1) {
                assert(r.dist[idx(v)] >= r.dist[idx(r.pred[idx(v)])] + edge_weight(edges, r.pred[idx(v)], v));
            }
        }
    }
    assert(with_negative > 100);  // la prueba sí ejercita los dos lados
}

void test_cycle_recovery() {
    std::mt19937 gen(2);
    int seen = 0;
    int off_cycle = 0;
    for (int t = 0; t < 4000; t++) {
        const int n = 2 + static_cast<int>(gen() % 7);
        const std::vector<Edge> edges = random_graph(n, 30, -4, 9, gen);
        const BellmanFordResult r = bellman_ford(n, edges, 0);
        const bool expected = has_negative_cycle(n, edges, reachable_from(n, edges, 0));
        assert((r.witness != -1) == expected);
        if (r.witness == -1) continue;
        seen++;
        const std::vector<int> cycle = recover_cycle(n, r.pred, r.witness);
        std::vector<int> sorted = cycle;
        std::sort(sorted.begin(), sorted.end());
        assert(std::adjacent_find(sorted.begin(), sorted.end()) == sorted.end());  // ciclo simple
        for (std::size_t i = 0; i < cycle.size(); i++) {
            assert(edge_weight(edges, cycle[i], cycle[(i + 1) % cycle.size()]) != INF);
        }
        assert(cycle_weight(edges, cycle) < 0);
        if (std::find(cycle.begin(), cycle.end(), r.witness) == cycle.end()) off_cycle++;
    }
    assert(seen > 300);
    assert(off_cycle > 0);  // el testigo no siempre está en el ciclo

    // El ejemplo de los apuntes: el nodo 5 cuelga del ciclo.
    const std::vector<Edge> tail = {{0, 1, 2}, {1, 2, 1}, {2, 3, 1}, {3, 1, -3}, {3, 4, 1}, {4, 5, 1}};
    const BellmanFordResult r = bellman_ford(6, tail, 0);
    assert(r.witness == 5);
    assert(cycle_weight(tail, recover_cycle(6, r.pred, r.witness)) == -1);
}

void test_unreached_nodes() {
    const std::vector<Edge> unreachable = {{0, 1, 7}, {2, 3, -5}, {3, 2, 1}};
    const BellmanFordResult r = bellman_ford(4, unreachable, 0);
    assert(r.witness == -1);
    assert(r.dist[2] == INF && r.dist[3] == INF);
    assert(r.dist[1] == 7);
    // n = 1, sin aristas.
    assert(bellman_ford(1, {}, 0).dist[0] == 0);
}

void test_week5_graph() {
    const std::vector<Edge> edges = {{0, 1, 1}, {0, 2, 4}, {2, 1, -5}, {1, 3, 2}};
    assert((dijkstra(4, edges, 0) == std::vector<long long>{0, -1, 4, 3}));
    assert((bellman_ford(4, edges, 0).dist == std::vector<long long>{0, -1, 4, 1}));
}

void test_decision_table() {
    std::mt19937 gen(3);
    for (int t = 0; t < 1500; t++) {
        const int n = 2 + static_cast<int>(gen() % 8);
        const std::vector<Edge> unit = random_graph(n, 30, 1, 1, gen);
        assert(bfs_distances(n, unit, 0) == bellman_ford(n, unit, 0).dist);
        const std::vector<Edge> positive = random_graph(n, 40, 0, 20, gen);
        assert(dijkstra(n, positive, 0) == bellman_ford(n, positive, 0).dist);
        std::vector<Edge> dag;
        for (const Edge& e : random_graph(n, 50, -9, 9, gen)) {
            if (e.from < e.to) dag.push_back(e);
        }
        assert(dag_distances(n, dag, 0) == bellman_ford(n, dag, 0).dist);
    }
}

void test_floyd_warshall() {
    std::mt19937 gen(4);
    int negative_graphs = 0;
    for (int t = 0; t < 3000; t++) {
        const int n = 2 + static_cast<int>(gen() % 5);
        const std::vector<Edge> edges = random_graph(n, 40, -3, 9, gen);
        const bool negative = has_negative_cycle(n, edges, std::vector<bool>(idx(n), true));
        const FloydResult f = floyd_warshall(n, edges);
        assert(f.negative_cycle == negative);  // la guardia ve todo ciclo negativo, y solo esos
        if (negative) {
            negative_graphs++;
            // Con la guardia, los valores siguen en la escala de un camino simple.
            for (const auto& row : f.dist) {
                for (long long x : row) assert(x == INF || (x >= -6 * n * 3 && x <= 6 * n * 9));
            }
            // Sin ella, también queda negativa la diagonal al final.
            Matrix m = initial_matrix(n, edges);
            for (int k = 0; k < n; k++) floyd_round(m, n, k);
            bool diagonal = false;
            for (int v = 0; v < n; v++) diagonal = diagonal || m[idx(v)][idx(v)] < 0;
            assert(diagonal);
            continue;
        }
        assert(f.dist == all_pairs_bellman_ford(n, edges));

        // d_k: comparar con la definición, y ver que la fila y la columna k no cambian en la vuelta k.
        Matrix m = initial_matrix(n, edges);
        for (int k = 0; k < n; k++) {
            const Matrix before = m;
            floyd_round(m, n, k);
            for (int v = 0; v < n; v++) {
                assert(m[idx(k)][idx(v)] == before[idx(k)][idx(v)]);
                assert(m[idx(v)][idx(k)] == before[idx(v)][idx(k)]);
            }
            // Dos filas, como en la definición: da lo mismo que en el mismo lugar.
            Matrix two_rows = before;
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    if (before[idx(i)][idx(k)] == INF || before[idx(k)][idx(j)] == INF) continue;
                    two_rows[idx(i)][idx(j)] =
                        std::min(two_rows[idx(i)][idx(j)], before[idx(i)][idx(k)] + before[idx(k)][idx(j)]);
                }
            }
            assert(two_rows == m);
        }
    }
    assert(negative_graphs > 100);

    const std::vector<Edge> trap = {{1, 3, 5}, {3, 2, 5}, {2, 0, 9}, {1, 0, 30}};
    assert(floyd_warshall(4, trap).dist[1][0] == 19);
    assert(floyd_warshall_k_inside(4, trap)[1][0] == 30);
}

void test_rates() {
    assert(log_sum({1.25, 0.8}) < 0);  // el producto real es 1, y la suma sale negativa
    assert(log_sum({2, 0.5}) == 0);
    assert(round_n_relaxes(5, consistent_table(), 0));  // falso arbitraje sin tolerancia
    assert(!round_n_relaxes(5, consistent_table(), 1e-9));
    assert(round_n_relaxes(400, long_cycle(400), 1e-9));  // 400 * 1e-9 < ln(1,000001)
    assert(!round_n_relaxes(400, long_cycle(400), 1e-8));  // 400 * 1e-8 > ln(1,000001)
}

void test_held_karp() {
    std::mt19937 gen(5);
    for (int t = 0; t < 400; t++) {
        const int n = 1 + static_cast<int>(gen() % 8);
        tsp::Costs cost(tsp::idx(n), std::vector<long long>(tsp::idx(n), 0));
        for (auto& row : cost) {
            for (long long& c : row) c = gen() % 5 == 0 ? tsp::NONE : static_cast<long long>(gen() % 50) + 1;
        }
        assert(tsp::held_karp(cost) == tsp::brute_force(cost));
    }
    for (int n = 3; n <= 9; n++) {
        tsp::Costs cost(tsp::idx(n), std::vector<long long>(tsp::idx(n), 1));
        long long tours = 0;
        long long transitions = 0;
        assert(tsp::brute_force(cost, &tours) == n);
        assert(tsp::held_karp(cost, &transitions) == n);
        long long factorial = 1;
        for (int i = 2; i < n; i++) factorial *= i;
        assert(tours == factorial);
        assert(transitions == (n - 1) * (n - 2) * (1LL << (n - 3)));
    }
}

}  // namespace

int main() {
    test_recurrence_and_one_row();
    test_cycle_recovery();
    test_unreached_nodes();
    test_week5_graph();
    test_decision_table();
    test_floyd_warshall();
    test_rates();
    test_held_karp();
    std::cout << "Todas las pruebas de la semana 10 pasaron.\n";
}