#pragma once

#include <algorithm>
#include <cmath>    // std::sqrt
#include <numeric>
#include <string>
#include <utility>  // std::swap, std::pair
#include <vector>

/**
 * Reads a symmetric TSP instance from a TSPLIB-format file.
 *
 * Pairwise distances are stored in a compact flat lower-triangular array
 * (n*(n-1)/2 entries) of doubles.  No integer rounding is applied — raw
 * floating-point values are kept throughout.
 *
 * Nodes are 0-indexed (0 … size()-1), regardless of the 1-based numbering
 * used inside the file.
 *
 * Supported EDGE_WEIGHT_TYPE:
 *   EUC_2D, EUC_3D, MAX_2D, MAX_3D, MAN_2D, MAN_3D, CEIL_2D, GEO, ATT,
 *   EXPLICIT
 *
 * Supported EDGE_WEIGHT_FORMAT (EXPLICIT only):
 *   FULL_MATRIX,
 *   UPPER_ROW,      LOWER_ROW,
 *   UPPER_DIAG_ROW, LOWER_DIAG_ROW,
 *   UPPER_COL,      LOWER_COL,
 *   UPPER_DIAG_COL, LOWER_DIAG_COL
 */
class TSPLIBGraph {
public:
    /// Load a TSPLIB file.  Throws std::runtime_error on any failure.
    static TSPLIBGraph load(const std::string& filepath);

    /// Load a plain-text edge-list file.
    /// Format: first line "n m", then m lines "u v cost" (0-indexed nodes).
    /// Missing edges get infinite cost.
    static TSPLIBGraph load_txt(const std::string& filepath);

    /// Number of nodes.
    int size() const noexcept { return n_; }

    /// Number of edges: n*(n-1)/2.
    size_t num_edges() const noexcept { return data_.size(); }

    // ── node-pair ↔ edge id ────────────────────────────────────────────────

    /// Unique edge id for the pair (u, v).  Same result for (v, u).
    /// Precondition: u != v, both in [0, size()).
    size_t edge_id(int u, int v) const noexcept {
        if (u < v) std::swap(u, v);          // ensure u > v
        return tri_index(u, v);
    }

    /// Endpoints {u, v} of edge e, with u > v.
    /// Precondition: e < num_edges().
    std::pair<int,int> edge_nodes(size_t e) const noexcept {
        // Closed-form row: i*(i-1)/2 <= e < i*(i+1)/2  →  i = floor((1+sqrt(1+8e))/2)
        int i = static_cast<int>((1.0 + std::sqrt(1.0 + 8.0 * static_cast<double>(e))) / 2.0);
        // Clamp for floating-point rounding (off by at most 1).
        while (tri_index(i, 0) > e) --i;
        while (tri_index(i + 1, 0) <= e) ++i;
        return {i, static_cast<int>(e - tri_index(i, 0))};
    }

    // ── distance accessors ─────────────────────────────────────────────────

    /// Distance between nodes u and v (0-based).  Returns 0.0 when u == v.
    /// Symmetric: dist(u,v) == dist(v,u).
    double dist(int u, int v) const noexcept {
        if (u == v) return 0.0;
        if (u < v) std::swap(u, v);          // ensure u > v
        return data_[tri_index(u, v)];
    }

    /// Distance (cost) of edge e.
    /// Precondition: e < num_edges().
    double dist(size_t e) const noexcept { return data_[e]; }

    /// Edge ids sorted by non-decreasing distance. Computed once on first call.
    const std::vector<size_t>& sorted_edges() const {
        if (sorted_edges_.empty() && !data_.empty()) {
            sorted_edges_.resize(data_.size());
            std::iota(sorted_edges_.begin(), sorted_edges_.end(), 0);
            std::sort(sorted_edges_.begin(), sorted_edges_.end(),
                      [this](size_t a, size_t b) { return data_[a] < data_[b]; });
        }
        return sorted_edges_;
    }

private:
    mutable std::vector<size_t> sorted_edges_;
    int n_ = 0;
    std::vector<double> data_;  // lower-triangular, n*(n-1)/2 elements

    /// Flat index for element (i, j) with i > j >= 0.
    static size_t tri_index(int i, int j) noexcept {
        return static_cast<size_t>(i) * (i - 1) / 2 + static_cast<size_t>(j);
    }

    void store(int u, int v, double d) {
        if (u == v) return;
        if (u < v) std::swap(u, v);
        data_[tri_index(u, v)] = d;
    }
};