// Nodos y tiempo de la búsqueda en los casos grandes. Una corrida de cada una, con -O2.
#include <chrono>
#include <iomanip>
#include <iostream>
#include <string>

#include "bound.h"
#include "geo.h"

namespace {

void row(const bound::Matrix& d, int n, const std::string& name, bound::Pruning pruning) {
    const auto start = std::chrono::steady_clock::now();
    const auto result = bound::solve(d, pruning);
    const std::chrono::duration<double> elapsed = std::chrono::steady_clock::now() - start;
    std::cout << std::setw(2) << n << "  " << std::left << std::setw(20) << name << std::right << std::setw(12)
              << result.nodes << std::setw(10) << std::fixed << std::setprecision(1)
              << static_cast<double>(result.cost) / 1000.0 << " km" << std::setw(9) << std::setprecision(2)
              << elapsed.count() << " s\n";
}

}  // namespace

int main() {
    const auto all = geo::read_places("capitales.csv");
    std::cout << " n  poda                       nodos       costo     tiempo\n";
    const auto d15 = geo::distance_matrix(geo::first_places(all, 15));
    row(d15, 15, "lo recorrido", bound::Pruning::partial);
    row(d15, 15, "con cota de árbol", bound::Pruning::spanning);
    row(geo::distance_matrix(geo::first_places(all, 20)), 20, "con cota de árbol", bound::Pruning::spanning);
    row(geo::distance_matrix(all), 26, "con cota de árbol", bound::Pruning::spanning);
}

// Cuántos nodos del árbol de búsqueda visita cada forma de podar, en las primeras filas de
// capitales.csv. Las tres dan el mismo costo; lo que cambia es el trabajo.
#include <iostream>
#include <string>

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

namespace {

std::string pad(const std::string& text, std::size_t width) {
    return std::string(width > text.size() ? width - text.size() : 0, ' ') + text;
}

std::string cell(const bound::Matrix& d, bound::Pruning pruning, int n, int up_to) {
    if (n > up_to) return "-";
    return std::to_string(bound::solve(d, pruning).nodes);
}

}  // namespace

int main() {
    const auto all = geo::read_places("capitales.csv");
    std::cout << " n  sin poda      lo recorrido  con cota de árbol  costo (m)\n";
    for (int n : {6, 8, 10, 11, 12, 13}) {
        const auto d = geo::distance_matrix(geo::first_places(all, n));
        std::cout << pad(std::to_string(n), 2) << pad(cell(d, bound::Pruning::none, n, 11), 10)
                  << pad(cell(d, bound::Pruning::partial, n, 13), 14)
                  << pad(cell(d, bound::Pruning::spanning, n, 13), 19)
                  << pad(std::to_string(bound::solve(d, bound::Pruning::spanning).cost), 11) << '\n';
    }
}

// Búsqueda exacta del viajero por backtracking, con tres formas de podar. Solo devuelve el costo y la
// cantidad de nodos visitados del árbol de búsqueda. Un recorrido parcial sale de 0 y crece por el final.
#pragma once
#include <algorithm>
#include <cstddef>
#include <limits>
#include <vector>

namespace bound {

using Matrix = std::vector<std::vector<long long>>;

inline constexpr long long NO_TOUR = std::numeric_limits<long long>::max() / 4;

inline std::size_t at(int i) { return static_cast<std::size_t>(i); }

enum class Pruning {
    none,      // recorre todo el árbol
    partial,   // poda si lo ya recorrido alcanza al mejor conocido
    spanning,  // lo ya recorrido, más una cota inferior de lo que falta
};

struct Result {
    long long cost = NO_TOUR;
    long long nodes = 0;
    bool finished = true;
};

// El peso de un árbol de expansión mínima de los nodos marcados, con Prim en O(r²).
inline long long spanning_weight(const Matrix& d, const std::vector<int>& nodes) {
    if (nodes.size() <= 1) return 0;
    std::vector<long long> best(nodes.size(), NO_TOUR);
    std::vector<char> inside(nodes.size(), 0);
    best[0] = 0;
    long long weight = 0;
    for (std::size_t round = 0; round < nodes.size(); round++) {
        std::size_t v = nodes.size();
        for (std::size_t u = 0; u < nodes.size(); u++) {
            if (inside[u] == 0 && (v == nodes.size() || best[u] < best[v])) v = u;
        }
        inside[v] = 1;
        weight += best[v];
        for (std::size_t u = 0; u < nodes.size(); u++) {
            const long long edge = d[at(nodes[v])][at(nodes[u])];
            if (inside[u] == 0 && edge < best[u]) best[u] = edge;
        }
    }
    return weight;
}

// Una cota inferior para terminar un recorrido que está en last, con los nodos de rest sin visitar.
// Cualquier forma de terminar sale de last a un nodo de rest, recorre rest en un camino y vuelve de
// un nodo de rest a 0. Ese camino es un árbol de expansión de rest, así que pesa al menos el mínimo.
inline long long completion_bound(const Matrix& d, int last, const std::vector<int>& rest) {
    if (rest.empty()) return d[at(last)][0];
    long long out = NO_TOUR;
    long long back = NO_TOUR;
    for (int r : rest) {
        out = std::min(out, d[at(last)][at(r)]);
        back = std::min(back, d[at(r)][0]);
    }
    return out + spanning_weight(d, rest) + back;
}

struct Search {
    const Matrix& d;
    Pruning pruning;
    long long node_limit;
    Result result;

    void visit(int last, std::vector<int>& rest, long long so_far) {
        if (!result.finished) return;
        result.nodes++;
        if (result.nodes > node_limit) {
            result.finished = false;
            return;
        }
        if (rest.empty()) {
            result.cost = std::min(result.cost, so_far + d[at(last)][0]);
            return;
        }
        if (pruning == Pruning::partial && so_far >= result.cost) return;
        if (pruning == Pruning::spanning && so_far + completion_bound(d, last, rest) >= result.cost) return;
        // Probar primero lo más cercano encuentra pronto un buen recorrido, y eso hace podar más.
        std::vector<int> order = rest;
        std::sort(order.begin(), order.end(), [&](int a, int b) { return d[at(last)][at(a)] < d[at(last)][at(b)]; });
        for (int next : order) {
            std::vector<int> remaining;
            remaining.reserve(rest.size() - 1);
            for (int r : rest) {
                if (r != next) remaining.push_back(r);
            }
            visit(next, remaining, so_far + d[at(last)][at(next)]);
        }
    }
};

// node_limit corta la búsqueda: finished queda en false y el costo no vale como óptimo.
inline Result solve(const Matrix& d, Pruning pruning, long long node_limit = std::numeric_limits<long long>::max()) {
    const int n = static_cast<int>(d.size());
    Search search{d, pruning, node_limit, Result{}};
    if (n <= 1) {
        search.result.cost = 0;
        return search.result;
    }
    std::vector<int> rest;
    for (int v = 1; v < n; v++) rest.push_back(v);
    search.visit(0, rest, 0);
    return search.result;
}

}  // namespace bound

capital,latitud,longitud
Abancay,-13.63426,-72.8838
Arequipa,-16.39899,-71.53747
Ayacucho,-13.1638,-74.22345
Cajamarca,-7.16378,-78.50027
Callao,-12.05162,-77.13452
Cerro de Pasco,-10.66577,-76.25309
Chachapoyas,-6.23169,-77.86903
Chiclayo,-6.77008,-79.85495
Cusco,-13.53188,-71.96701
Huacho,-11.11718,-77.60584
Huancavelica,-12.78691,-74.97298
Huancayo,-12.06866,-75.21027
Huaraz,-9.52614,-77.52869
Huánuco,-9.92882,-76.23989
Ica,-14.07538,-75.73422
Iquitos,-3.74814,-73.2529
Lima,-12.04318,-77.02824
Moquegua,-17.19747,-70.93587
Moyobamba,-6.03441,-76.97423
Piura,-5.18192,-80.65715
Pucallpa,-8.37915,-74.55387
Puerto Maldonado,-12.58946,-69.19948
Puno,-15.84003,-70.02198
Tacna,-18.01465,-70.25362
Trujillo,-8.11599,-79.02998
Tumbes,-3.55592,-80.44258

# Fuente de `capitales.csv`

Coordenadas de [GeoNames](https://www.geonames.org/), archivo `PE.zip` de
<https://download.geonames.org/export/dump/>, descargado el 2026-11-04. Licencia CC BY 4.0.

Son los 26 lugares de Perú con código `PPLA` (sede de una división de primer nivel) o `PPLC`
(capital del país): las sedes de los gobiernos regionales, incluidas Huacho y Callao, y Lima.
Solo se conservaron nombre, latitud y longitud, sin otros cambios.

// La comparación de la semana: largo y tiempo de cada método con 10, 15, 20 y 26 ciudades. Cada
// tiempo es el menor de tres corridas, con -O2. El 2-opt de aquí es el lento de tour.h.
#include <algorithm>
#include <chrono>
#include <functional>
#include <iomanip>
#include <iostream>
#include <string>

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

namespace {

struct Measured {
    long long length;
    double seconds;
};

Measured measure(int runs, const std::function<long long()>& method) {
    Measured best{0, 1e18};
    for (int r = 0; r < runs; r++) {
        const auto start = std::chrono::steady_clock::now();
        best.length = method();
        const std::chrono::duration<double> elapsed = std::chrono::steady_clock::now() - start;
        best.seconds = std::min(best.seconds, elapsed.count());
    }
    return best;
}

void row(const std::string& name, const Measured& m, long long reference) {
    std::cout << "  " << std::left << std::setw(26) << name << std::right << std::setw(9) << std::fixed
              << std::setprecision(1) << static_cast<double>(m.length) / 1000.0 << " km" << std::setw(9)
              << std::setprecision(3) << static_cast<double>(m.length) / static_cast<double>(reference)
              << std::setw(10) << m.seconds << " s\n";
}

}  // namespace

int main() {
    const auto all = geo::read_places("capitales.csv");
    for (int n : {10, 15, 20, 26}) {
        const auto d = geo::distance_matrix(geo::first_places(all, n));
        // El óptimo viene de la búsqueda con cota de árbol, que es exacta. Es el divisor de la razón.
        const auto exact = measure(n >= 20 ? 1 : 3, [&] { return bound::solve(d, bound::Pruning::spanning).cost; });
        std::cout << "n = " << n << "                       largo     / óptimo    tiempo\n";
        row("árbol de expansión mínima", measure(3, [&] { return tour::minimum_spanning_tree(d).weight; }), exact.length);
        row("árbol duplicado", measure(3, [&] { return tour::tour_length(d, tour::double_tree(d)); }), exact.length);
        row("duplicado + 2-opt", measure(3, [&] {
            auto t = tour::double_tree(d);
            tour::two_opt_plain(d, t);
            return tour::tour_length(d, t);
        }), exact.length);
        row("vecino más cercano", measure(3, [&] { return tour::tour_length(d, tour::nearest_neighbor(d)); }), exact.length);
        row("vecino + 2-opt", measure(3, [&] {
            auto t = tour::nearest_neighbor(d);
            tour::two_opt_plain(d, t);
            return tour::tour_length(d, t);
        }), exact.length);
        if (n <= 20) row("Held-Karp", measure(3, [&] { return tsp::held_karp(d); }), exact.length);
        row("búsqueda con cota", exact, exact.length);
        std::cout << '\n';
    }
}

// Tres reglas para la cobertura de vértices: la que falla, la que tiene factor 2 y la natural.
#include <algorithm>
#include <cstddef>
#include <iostream>
#include <random>
#include <string>

#include "cover.h"

namespace {

std::string list(const cover::Nodes& nodes) {
    std::string text = "{";
    for (std::size_t i = 0; i < nodes.size(); i++) text += (i > 0 ? ", " : "") + std::to_string(nodes[i]);
    return text + "}";
}

void compare_on(const std::string& title, int n, const cover::Edges& edges) {
    std::cout << title << '\n';
    std::cout << "  una sola punta: " << list(cover::one_endpoint_cover(n, edges)) << '\n';
    const auto matched = cover::matching_cover(n, edges);
    std::cout << "  matching maximal, " << matched.matching.size() << " aristas: " << list(matched.cover) << '\n';
    std::cout << "  mayor grado: " << list(cover::greedy_degree_cover(n, edges)) << '\n';
    std::cout << "  óptimo por fuerza bruta: " << cover::minimum_cover_size(n, edges) << "\n\n";
}

}  // namespace

int main() {
    compare_on("Estrella con centro 0 y 6 hojas:", 7, {{0, 1}, {0, 2}, {0, 3}, {0, 4}, {0, 5}, {0, 6}});
    compare_on("Camino 0-1-2-3, primero la arista 0-1:", 4, {{0, 1}, {1, 2}, {2, 3}});
    compare_on("Camino 0-1-2-3, primero la arista 1-2:", 4, {{1, 2}, {0, 1}, {2, 3}});
    compare_on("Una sola arista:", 2, {{0, 1}});

    std::cout << "La trampa del mayor grado:\n";
    std::cout << "   k  nodos  aristas  mayor grado  matching  mayor grado / k\n";
    for (int k : {6, 12, 24, 48, 96, 192}) {
        const auto trap = cover::greedy_trap(k);
        const auto chosen = cover::greedy_degree_cover(trap.nodes, trap.edges);
        const auto matched = cover::matching_cover(trap.nodes, trap.edges);
        std::cout << std::string(k < 10 ? 3 : (k < 100 ? 2 : 1), ' ') << k << "  " << trap.nodes << "  "
                  << trap.edges.size() << "  " << chosen.size() << "  " << matched.cover.size() << "  "
                  << static_cast<double>(chosen.size()) / k << '\n';
    }

    std::cout << "\n300 grafos al azar, 16 nodos, cada arista con probabilidad 1/4:\n";
    std::mt19937 rng(15);
    long long sum_opt = 0;
    long long sum_matching = 0;
    long long sum_greedy = 0;
    double worst_matching = 0;
    double worst_greedy = 0;
    const int graphs = 300;
    for (int g = 0; g < graphs; g++) {
        cover::Edges edges;
        for (int u = 0; u < 16; u++) {
            for (int v = u + 1; v < 16; v++) {
                if (rng() % 4 == 0) edges.push_back({u, v});
            }
        }
        const int opt = cover::minimum_cover_size(16, edges);
        const auto matching = cover::matching_cover(16, edges).cover.size();
        const auto greedy = cover::greedy_degree_cover(16, edges).size();
        sum_opt += opt;
        sum_matching += static_cast<long long>(matching);
        sum_greedy += static_cast<long long>(greedy);
        if (opt > 0) {
            worst_matching = std::max(worst_matching, static_cast<double>(matching) / opt);
            worst_greedy = std::max(worst_greedy, static_cast<double>(greedy) / opt);
        }
    }
    std::cout << "  promedio del óptimo: " << static_cast<double>(sum_opt) / graphs << '\n';
    std::cout << "  promedio del matching: " << static_cast<double>(sum_matching) / graphs
              << ", peor razón " << worst_matching << '\n';
    std::cout << "  promedio del mayor grado: " << static_cast<double>(sum_greedy) / graphs
              << ", peor razón " << worst_greedy << '\n';
}

// Cobertura de vértices: un método con factor 2 demostrado, dos reglas parecidas y un árbitro.
// Un grafo es una cantidad de nodos y una lista de aristas.
#pragma once
#include <cstddef>
#include <utility>
#include <vector>

namespace cover {

using Edge = std::pair<int, int>;
using Edges = std::vector<Edge>;
using Nodes = std::vector<int>;

inline std::size_t at(int i) { return static_cast<std::size_t>(i); }

inline bool is_cover(const Nodes& chosen, int n, const Edges& edges) {
    std::vector<char> in(at(n), 0);
    for (int v : chosen) in[at(v)] = 1;
    for (const auto& [u, v] : edges) {
        if (in[at(u)] == 0 && in[at(v)] == 0) return false;
    }
    return true;
}

struct Matched {
    Edges matching;
    Nodes cover;
};

// Revisa las aristas en orden. Si ninguna punta de la arista está elegida, entra al matching y se
// eligen las dos. Invariante: lo elegido son las puntas del matching, y las aristas ya revisadas
// tienen alguna punta elegida. Al final el matching es maximal: no le cabe ninguna arista más.
inline Matched matching_cover(int n, const Edges& edges) {
    std::vector<char> used(at(n), 0);
    Matched result;
    for (const auto& [u, v] : edges) {
        if (used[at(u)] == 0 && used[at(v)] == 0) {
            used[at(u)] = used[at(v)] = 1;
            result.matching.push_back({u, v});
            result.cover.push_back(u);
            result.cover.push_back(v);
        }
    }
    return result;
}

// Una idea parecida que falla: por cada arista sin cubrir, elegir solo una punta, la segunda.
inline Nodes one_endpoint_cover(int n, const Edges& edges) {
    std::vector<char> used(at(n), 0);
    Nodes chosen;
    for (const auto& [u, v] : edges) {
        if (used[at(u)] == 0 && used[at(v)] == 0) {
            used[at(v)] = 1;
            chosen.push_back(v);
        }
    }
    return chosen;
}

// La regla natural: repetir con el nodo que cubre más aristas sin cubrir. En un empate, el de menor número.
inline Nodes greedy_degree_cover(int n, const Edges& edges) {
    std::vector<char> alive(edges.size(), 1);
    Nodes chosen;
    while (true) {
        std::vector<int> degree(at(n), 0);
        for (std::size_t e = 0; e < edges.size(); e++) {
            if (alive[e] != 0) {
                degree[at(edges[e].first)]++;
                degree[at(edges[e].second)]++;
            }
        }
        int best = 0;
        for (int v = 1; v < n; v++) {
            if (degree[at(v)] > degree[at(best)]) best = v;
        }
        if (n == 0 || degree[at(best)] == 0) return chosen;
        chosen.push_back(best);
        for (std::size_t e = 0; e < edges.size(); e++) {
            if (edges[e].first == best || edges[e].second == best) alive[e] = 0;
        }
    }
}

// El árbitro: prueba los 2ⁿ subconjuntos. Solo sirve para grafos chicos.
inline int minimum_cover_size(int n, const Edges& edges) {
    int best = n;
    for (unsigned mask = 0; mask < (1U << n); mask++) {
        bool covers = true;
        for (const auto& [u, v] : edges) {
            if (((mask >> u) & 1U) == 0 && ((mask >> v) & 1U) == 0) {
                covers = false;
                break;
            }
        }
        if (covers) {
            int size = 0;
            for (int v = 0; v < n; v++) size += static_cast<int>((mask >> v) & 1U);
            if (size < best) best = size;
        }
    }
    return best;
}

// Una trampa para la regla del mayor grado. Hay k nodos a la izquierda, 0, ..., k - 1. Para cada i de
// 2 a k hay ⌊k/i⌋ nodos a la derecha, y cada uno se une a i nodos de la izquierda, sin repetirlos
// dentro del grupo. Los k nodos de la izquierda ya son una cobertura, así que el óptimo es a lo sumo k.
struct Trap {
    int nodes;
    Edges edges;
    int right_nodes;
};

inline Trap greedy_trap(int k) {
    Trap trap{k, {}, 0};
    for (int i = 2; i <= k; i++) {
        for (int group = 0; group < k / i; group++) {
            const int r = trap.nodes++;
            trap.right_nodes++;
            for (int offset = 0; offset < i; offset++) trap.edges.push_back({group * i + offset, r});
        }
    }
    return trap;
}

}  // namespace cover

// El árbol duplicado, paso a paso, con ocho capitales.
#include <cstddef>
#include <iostream>
#include <string>
#include <vector>

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

namespace {

using geo::at;

std::vector<geo::Place> chosen_places() {
    const auto all = geo::read_places("capitales.csv");
    std::vector<geo::Place> places;
    for (const std::string name : {"Lima", "Callao", "Huacho", "Huancayo", "Ica", "Ayacucho", "Cusco", "Arequipa"}) {
        for (const auto& place : all) {
            if (place.name == name) places.push_back(place);
        }
    }
    return places;
}

std::string km(long long metres) {
    const long long tenths = (metres + 50) / 100;
    return std::to_string(tenths / 10) + "," + std::to_string(tenths % 10) + " km";
}

}  // namespace

int main() {
    const auto places = chosen_places();
    const auto d = geo::distance_matrix(places);
    const int n = static_cast<int>(places.size());
    const auto tree = tour::minimum_spanning_tree(d);

    std::cout << "Aristas del árbol de expansión mínima, con raíz en " << places[0].name << ":\n";
    for (int v = 0; v < n; v++) {
        for (int child : tree.children[at(v)]) {
            std::cout << "  " << places[at(v)].name << " - " << places[at(child)].name << ": " << km(d[at(v)][at(child)]) << '\n';
        }
    }
    std::cout << "peso del árbol: " << km(tree.weight) << "\n\n";

    const auto walk = tour::euler_walk(tree);
    std::cout << "Caminata, cada arista dos veces:\n ";
    for (int v : walk) std::cout << ' ' << places[at(v)].name;
    std::cout << "\nlargo de la caminata: " << km(tour::walk_length(d, walk)) << "\n\n";

    const auto shortened = tour::shortcut(walk, n);
    std::cout << "Atajando las repeticiones:\n ";
    for (int v : shortened) std::cout << ' ' << places[at(v)].name;
    std::cout << " y de vuelta a " << places[at(shortened[0])].name << '\n';
    const long long length = tour::tour_length(d, shortened);
    const long long optimum = tsp::held_karp(d);
    std::cout << "largo del recorrido: " << km(length) << '\n';
    std::cout << "óptimo, con Held-Karp: " << km(optimum) << '\n';
    std::cout << "dos veces el árbol: " << km(2 * tree.weight) << "\n\n";
}

// Distancias entre lugares de la Tierra. El gran círculo es el camino más corto por la superficie
// de una esfera de radio 6371 km: no es una carretera ni una ruta aérea real.
// Cada distancia se guarda en metros enteros, redondeada hacia arriba.
#pragma once
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <fstream>
#include <numbers>
#include <stdexcept>
#include <string>
#include <vector>

namespace geo {

inline constexpr double EARTH_RADIUS_M = 6371000.0;

struct Place {
    std::string name;
    double latitude;
    double longitude;
};

using Matrix = std::vector<std::vector<long long>>;

inline std::size_t at(int i) { return static_cast<std::size_t>(i); }

// La fórmula de haversine, en metros. Se porta bien con puntos muy cercanos.
inline double great_circle_m(const Place& a, const Place& b) {
    const double to_radians = std::numbers::pi / 180.0;
    const double phi1 = a.latitude * to_radians;
    const double phi2 = b.latitude * to_radians;
    const double half_dphi = (phi2 - phi1) / 2.0;
    const double half_dlambda = (b.longitude - a.longitude) * to_radians / 2.0;
    const double h = std::sin(half_dphi) * std::sin(half_dphi) +
                     std::cos(phi1) * std::cos(phi2) * std::sin(half_dlambda) * std::sin(half_dlambda);
    return 2.0 * EARTH_RADIUS_M * std::asin(std::sqrt(std::min(1.0, h)));
}

inline long long distance_m(const Place& a, const Place& b) {
    return static_cast<long long>(std::ceil(great_circle_m(a, b)));
}

// Lo que ocurre si se tratan los grados como un plano: una mala distancia, a propósito.
inline long long flat_m(const Place& a, const Place& b) {
    const double dlat = a.latitude - b.latitude;
    const double dlon = a.longitude - b.longitude;
    const double metres_per_degree = EARTH_RADIUS_M * std::numbers::pi / 180.0;
    return static_cast<long long>(std::ceil(std::sqrt(dlat * dlat + dlon * dlon) * metres_per_degree));
}

inline Matrix distance_matrix(const std::vector<Place>& places) {
    const int n = static_cast<int>(places.size());
    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)] = distance_m(places[at(i)], places[at(j)]);
        }
    }
    return d;
}

// Lee capitales.csv: una cabecera y luego nombre,latitud,longitud.
inline std::vector<Place> read_places(const std::string& path) {
    std::ifstream file(path);
    if (!file) throw std::runtime_error("no se pudo abrir " + path);
    std::vector<Place> places;
    std::string line;
    std::getline(file, line);
    while (std::getline(file, line)) {
        const std::size_t first = line.find(',');
        const std::size_t second = line.find(',', first + 1);
        places.push_back({line.substr(0, first), std::stod(line.substr(first + 1, second - first - 1)),
                          std::stod(line.substr(second + 1))});
    }
    return places;
}

// Las primeras n filas del archivo.
inline std::vector<Place> first_places(const std::vector<Place>& places, int n) {
    return {places.begin(), places.begin() + n};
}

}  // namespace geo

// La reducción que descarta una aproximación sin desigualdad triangular. Si un algoritmo tuviera
// factor alfa = 2, comparar su largo con alfa·n decidiría si hay ciclo hamiltoniano.
#include <cstddef>
#include <iostream>
#include <string>
#include <utility>
#include <vector>

#include "../week-10/held_karp.h"
#include "tour.h"

namespace {

using tour::at;
using EdgeList = std::vector<std::pair<int, int>>;

tour::Matrix weights_from_graph(int n, const EdgeList& edges, long long heavy) {
    tour::Matrix w(at(n), std::vector<long long>(at(n), heavy));
    for (int v = 0; v < n; v++) w[at(v)][at(v)] = 0;
    for (const auto& [u, v] : edges) w[at(u)][at(v)] = w[at(v)][at(u)] = 1;
    return w;
}

void show(const std::string& title, int n, const EdgeList& edges) {
    const int alpha = 2;
    const long long heavy = alpha * n + 1;
    const auto w = weights_from_graph(n, edges, heavy);
    const long long optimum = tsp::held_karp(w);
    const auto doubled = tour::double_tree(w);
    std::cout << title << '\n';
    std::cout << "  n = " << n << ", W = " << heavy << ", alfa·n = " << alpha * n << '\n';
    std::cout << "  óptimo: " << optimum << (optimum <= alpha * n ? " <= " : " > ") << alpha * n << '\n';
    std::cout << "  árbol duplicado:";
    for (int v : doubled) std::cout << ' ' << v;
    std::cout << ", largo " << tour::tour_length(w, doubled) << "\n\n";
}

}  // namespace

int main() {
    show("Un ciclo de 6 nodos con la numeración desordenada: sí hay ciclo hamiltoniano", 6,
         {{0, 3}, {3, 1}, {1, 5}, {5, 2}, {2, 4}, {4, 0}});
    show("Dos nodos unidos a otros cuatro, K(2,4): no hay ciclo hamiltoniano", 6,
         {{0, 2}, {0, 3}, {0, 4}, {0, 5}, {1, 2}, {1, 3}, {1, 4}, {1, 5}});
}

// Cuánta memoria toma de verdad una tabla de Held-Karp. Reserva cada tabla, la llena y mira el
// máximo de memoria del proceso (ru_maxrss). Van de menor a mayor porque ese máximo nunca baja.
#include <sys/resource.h>

#include <cstddef>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

double peak_mb() {
    rusage usage{};
    getrusage(RUSAGE_SELF, &usage);
    return static_cast<double>(usage.ru_maxrss) / 1024.0;
}

template <typename Cell>
void flat_table(const std::string& name, std::size_t rows, std::size_t columns) {
    std::vector<Cell> table(rows * columns, Cell{1});
    std::cout << std::left << std::setw(46) << name << std::right << std::setw(9) << std::fixed
              << std::setprecision(1) << static_cast<double>(table.size() * sizeof(Cell)) / (1024.0 * 1024.0)
              << " MiB" << std::setw(9) << peak_mb() << " MiB  " << (table[rows / 2] == Cell{1} ? "" : "?") << '\n';
}

}  // namespace

int main() {
    const std::size_t half = std::size_t{1} << 19;
    const std::size_t full = std::size_t{1} << 20;
    std::cout << std::left << std::setw(46) << "tabla para n = 20" << std::right << std::setw(13) << "la cuenta"
              << std::setw(13) << "pico" << '\n';
    flat_table<std::int32_t>("2^19 x 20, int", half, 20);
    flat_table<std::int64_t>("2^19 x 20, long long", half, 20);
    flat_table<std::int32_t>("2^20 x 20, int", full, 20);
    flat_table<std::int64_t>("2^20 x 20, long long", full, 20);
    {
        std::vector<std::vector<std::int64_t>> table(full, std::vector<std::int64_t>(20, 1));
        std::cout << std::left << std::setw(46) << "2^20 filas de 20, vector de vectores, long long" << std::right
                  << std::setw(9) << static_cast<double>(full * 20 * 8) / (1024.0 * 1024.0) << " MiB" << std::setw(9)
                  << peak_mb() << " MiB  " << (table[7][3] == 1 ? "" : "?") << '\n';
    }
}

// Comprueba la desigualdad triangular en los 26 lugares de capitales.csv, con tres formas de guardar
// la distancia, y muestra por qué redondear al entero más cercano puede romperla.
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <iostream>
#include <numbers>
#include <string>
#include <vector>

#include "geo.h"

namespace {

using geo::at;
using geo::Place;

struct Check {
    long long triples = 0;
    long long broken = 0;
    double worst = 0;  // cuánto le sobra a un lado, en metros, en el peor caso
};

template <typename Distance>
Check check_triples(const std::vector<Place>& places, Distance distance) {
    Check check;
    const int n = static_cast<int>(places.size());
    for (int a = 0; a < n; a++) {
        for (int b = 0; b < n; b++) {
            for (int c = 0; c < n; c++) {
                if (a == b || b == c || a == c) continue;
                check.triples++;
                const double excess = static_cast<double>(distance(places[at(a)], places[at(c)])) -
                                      static_cast<double>(distance(places[at(a)], places[at(b)])) -
                                      static_cast<double>(distance(places[at(b)], places[at(c)]));
                if (excess > 0) {
                    check.broken++;
                    check.worst = std::max(check.worst, excess);
                }
            }
        }
    }
    return check;
}

void show(const std::string& title, const Check& check) {
    std::cout << "  " << title << ": " << check.broken << " de " << check.triples << " ternas la rompen";
    if (check.broken > 0) std::cout << ", peor exceso " << check.worst << " m";
    std::cout << '\n';
}

}  // namespace

int main() {
    const auto places = geo::read_places("capitales.csv");
    const int n = static_cast<int>(places.size());
    std::cout << std::fixed << std::setprecision(3);

    double nearest = 1e12;
    double farthest = 0;
    int near_a = 0;
    int near_b = 0;
    int far_a = 0;
    int far_b = 0;
    double to_integer = 1;
    long double formula_gap = 0;
    for (int a = 0; a < n; a++) {
        for (int b = a + 1; b < n; b++) {
            const double value = geo::great_circle_m(places[at(a)], places[at(b)]);
            if (value < nearest) {
                nearest = value;
                near_a = a;
                near_b = b;
            }
            if (value > farthest) {
                farthest = value;
                far_a = a;
                far_b = b;
            }
            to_integer = std::min({to_integer, value - std::floor(value), std::ceil(value) - value});
            const double flat = static_cast<double>(geo::flat_m(places[at(a)], places[at(b)]));
            formula_gap = std::max<long double>(formula_gap, std::abs(flat - value) / value);
        }
    }
    std::cout << n << " lugares, " << n * (n - 1) / 2 << " pares\n";
    std::cout << "más cercanos: " << places[at(near_a)].name << " y " << places[at(near_b)].name << ", " << nearest << " m\n";
    std::cout << "más lejanos: " << places[at(far_a)].name << " y " << places[at(far_b)].name << ", " << farthest << " m\n";
    std::cout << std::setprecision(6) << "lo más cerca de un entero que queda una distancia: " << to_integer << " m\n"
              << std::setprecision(3);
    std::cout << "peor error relativo de tratar los grados como un plano: " << 100.0L * formula_gap << " %\n\n";

    std::cout << "Desigualdad triangular, d(a,c) <= d(a,b) + d(b,c), en todas las ternas:\n";
    show("double, en metros", check_triples(places, geo::great_circle_m));
    show("entero, redondeado hacia arriba", check_triples(places, geo::distance_m));
    show("entero, al más cercano", check_triples(places, [](const Place& a, const Place& b) {
        return std::llround(geo::great_circle_m(a, b));
    }));

    double tightest = 1e12;
    int ta = 0;
    int tb = 0;
    int tc = 0;
    for (int a = 0; a < n; a++) {
        for (int b = 0; b < n; b++) {
            for (int c = 0; c < n; c++) {
                if (a == b || b == c || a == c) continue;
                const double slack = geo::great_circle_m(places[at(a)], places[at(b)]) +
                                     geo::great_circle_m(places[at(b)], places[at(c)]) -
                                     geo::great_circle_m(places[at(a)], places[at(c)]);
                if (slack < tightest) {
                    tightest = slack;
                    ta = a;
                    tb = b;
                    tc = c;
                }
            }
        }
    }
    std::cout << "la terna más justa: " << places[at(ta)].name << ", " << places[at(tb)].name << ", "
              << places[at(tc)].name << ", le sobran " << tightest << " m\n\n";

    // Tres puntos sobre el ecuador, separados 1000,4 m: el redondeo al más cercano rompe la desigualdad.
    const double degrees = 1000.4 / (geo::EARTH_RADIUS_M * std::numbers::pi / 180.0);
    const Place a{"a", 0, 0};
    const Place b{"b", 0, degrees};
    const Place c{"c", 0, 2 * degrees};
    std::cout << "Tres puntos sobre el ecuador:\n";
    std::cout << "  reales: ab = " << geo::great_circle_m(a, b) << ", bc = " << geo::great_circle_m(b, c)
              << ", ac = " << geo::great_circle_m(a, c) << '\n';
    std::cout << "  al más cercano: " << std::llround(geo::great_circle_m(a, b)) << " + "
              << std::llround(geo::great_circle_m(b, c)) << " < " << std::llround(geo::great_circle_m(a, c)) << '\n';
    std::cout << "  hacia arriba: " << geo::distance_m(a, b) << " + " << geo::distance_m(b, c) << " >= "
              << geo::distance_m(a, c) << '\n';
}