main.cpp

/*
 * Pachner Move Simulator
 * =====================
 * 
 * Simulates 2-2 Pachner moves (diagonal flips) on a triangulated grid.
 * Connects to MCMC research: each triangulation is a state, each flip
 * is a Markov chain transition, and path length serves as the Hamiltonian.
 *
 * Usage:
 *   ./pachner_sim                    Interactive mode
 *   ./pachner_sim -m <steps>         MCMC sampling mode
 *   ./pachner_sim -s <rows> <cols>   Custom grid size
 */

#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include <cmath>
#include <numeric>
#include <algorithm>
#include <map>
#include <random>
#include <ctime>

#include "PachnerGrid.h"
#include "PathFinder.h"
#include "Renderer.h"

// MCMC with Metropolis-Hastings criterion
struct MCMCResult {
    std::vector<int> pathLengths;
    int acceptedMoves;
    int rejectedMoves;
    double meanPathLength;
    double variance;
    int minPathLength;
    int maxPathLength;
};

MCMCResult runMCMC(PachnerGrid& grid, int steps, double temperature = 1.0) {
    MCMCResult result;
    result.acceptedMoves = 0;
    result.rejectedMoves = 0;

    std::mt19937 rng(static_cast<unsigned>(std::time(nullptr)));
    std::uniform_real_distribution<double> uniform(0.0, 1.0);

    auto currentPath = PathFinder::findShortestPath(grid, grid.startVertex(), grid.endVertex());
    int currentEnergy = static_cast<int>(currentPath.size()) - 1;

    for (int step = 0; step < steps; ++step) {
        std::uniform_int_distribution<int> rowDist(0, grid.rows - 1);
        std::uniform_int_distribution<int> colDist(0, grid.cols - 1);
        int r = rowDist(rng);
        int c = colDist(rng);

        grid.flip(r, c);

        auto newPath = PathFinder::findShortestPath(grid, grid.startVertex(), grid.endVertex());
        int newEnergy = newPath.empty() ? 999 : static_cast<int>(newPath.size()) - 1;

        int deltaE = newEnergy - currentEnergy;
        bool accept = false;

        if (deltaE <= 0) {
            accept = true;
        } else {
            double prob = std::exp(-static_cast<double>(deltaE) / temperature);
            accept = (uniform(rng) < prob);
        }

        if (accept) {
            currentEnergy = newEnergy;
            currentPath = newPath;
            result.acceptedMoves++;
        } else {
            grid.flip(r, c);
            result.rejectedMoves++;
        }

        result.pathLengths.push_back(currentEnergy);
    }

    double sum = std::accumulate(result.pathLengths.begin(), result.pathLengths.end(), 0.0);
    result.meanPathLength = sum / result.pathLengths.size();

    double sqSum = 0;
    for (int pl : result.pathLengths) {
        sqSum += (pl - result.meanPathLength) * (pl - result.meanPathLength);
    }
    result.variance = sqSum / result.pathLengths.size();

    result.minPathLength = *std::min_element(result.pathLengths.begin(), result.pathLengths.end());
    result.maxPathLength = *std::max_element(result.pathLengths.begin(), result.pathLengths.end());

    return result;
}

void printMCMCStats(const MCMCResult& result) {
    std::cout << "\n";
    std::cout << Color::WHITE << Color::BOLD << "  +======================================+\n";
    std::cout << "  |       MCMC Sampling Results          |\n";
    std::cout << "  +======================================+" << Color::RESET << "\n";
    
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Total Steps:      " << std::setw(8) << result.pathLengths.size() 
              << Color::CYAN << "         |\n";
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Accepted Moves:   " << std::setw(8) << result.acceptedMoves 
              << Color::CYAN << "         |\n";
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Rejected Moves:   " << std::setw(8) << result.rejectedMoves 
              << Color::CYAN << "         |\n";
    
    double acceptRate = 100.0 * result.acceptedMoves / 
                        (result.acceptedMoves + result.rejectedMoves);
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Acceptance Rate:  " << std::setw(7) << std::fixed << std::setprecision(1) 
              << acceptRate << "%" 
              << Color::CYAN << "         |\n";
    
    std::cout << Color::CYAN << "  +--------------------------------------+\n";
    std::cout << Color::CYAN << "  | " << Color::BRIGHT_BLUE 
              << "Mean Path Length: " << std::setw(8) << std::fixed << std::setprecision(2)
              << result.meanPathLength 
              << Color::CYAN << "         |\n";
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Variance:         " << std::setw(8) << std::fixed << std::setprecision(2)
              << result.variance 
              << Color::CYAN << "         |\n";
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Min Path Length:  " << std::setw(8) << result.minPathLength 
              << Color::CYAN << "         |\n";
    std::cout << Color::CYAN << "  | " << Color::WHITE 
              << "Max Path Length:  " << std::setw(8) << result.maxPathLength 
              << Color::CYAN << "         |\n";
    
    std::cout << Color::WHITE << Color::BOLD 
              << "  +======================================+" << Color::RESET << "\n";

    // Path length histogram
    std::map<int, int> histogram;
    for (int pl : result.pathLengths) {
        histogram[pl]++;
    }

    int maxCount = 0;
    for (auto& [len, count] : histogram) {
        maxCount = std::max(maxCount, count);
    }

    std::cout << "\n" << Color::WHITE << Color::BOLD 
              << "  Path Length Distribution (Hamiltonian Energy Levels):" 
              << Color::RESET << "\n\n";

    for (auto& [len, count] : histogram) {
        int barLen = (count * 40) / maxCount;
        double pct = 100.0 * count / result.pathLengths.size();
        
        std::cout << Color::WHITE << "  E=" << std::setw(2) << len << " ";
        std::cout << Color::BRIGHT_BLUE;
        for (int i = 0; i < barLen; ++i) std::cout << "#";
        std::cout << Color::DIM << " " << std::fixed << std::setprecision(1) << pct << "%";
        std::cout << Color::RESET << "\n";
    }
    std::cout << "\n";
}

int main(int argc, char* argv[]) {
    int gridRows = 4;
    int gridCols = 4;
    int mcmcSteps = 0;
    double temperature = 1.0;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "-m" && i + 1 < argc) {
            mcmcSteps = std::stoi(argv[++i]);
        } else if (arg == "-s" && i + 2 < argc) {
            gridRows = std::stoi(argv[++i]);
            gridCols = std::stoi(argv[++i]);
        } else if (arg == "-t" && i + 1 < argc) {
            temperature = std::stod(argv[++i]);
        } else if (arg == "-h" || arg == "--help") {
            std::cout << "Pachner Move Simulator\n\n";
            std::cout << "Usage: ./pachner_sim [options]\n\n";
            std::cout << "Options:\n";
            std::cout << "  -m <steps>        Run MCMC sampling with N steps\n";
            std::cout << "  -s <rows> <cols>  Set grid size (default: 4 4)\n";
            std::cout << "  -t <temp>         Set temperature for Metropolis (default: 1.0)\n";
            std::cout << "  -h, --help        Show this help message\n\n";
            std::cout << "Interactive commands:\n";
            std::cout << "  f <r> <c>   Flip diagonal at cell (r, c)\n";
            std::cout << "  r           Random flip\n";
            std::cout << "  m <N>       Run N MCMC steps\n";
            std::cout << "  p           Print current path\n";
            std::cout << "  s           Print grid statistics\n";
            std::cout << "  reset       Reset to initial state\n";
            std::cout << "  q           Quit\n";
            return 0;
        }
    }

    PachnerGrid grid(gridRows, gridCols);

    // MCMC batch mode
    if (mcmcSteps > 0) {
        std::cout << Color::WHITE << Color::BOLD 
                  << "\n  Running MCMC with " << mcmcSteps << " steps, T=" << temperature 
                  << "...\n" << Color::RESET;
        
        auto result = runMCMC(grid, mcmcSteps, temperature);
        
        auto finalPath = PathFinder::findShortestPath(grid, grid.startVertex(), grid.endVertex());
        Renderer::render(grid, finalPath);
        printMCMCStats(result);
        return 0;
    }

    // Interactive mode
    while (true) {
        auto path = PathFinder::findShortestPath(grid, grid.startVertex(), grid.endVertex());
        Renderer::render(grid, path);

        std::string line;
        if (!std::getline(std::cin, line)) break;

        std::istringstream iss(line);
        std::string cmd;
        iss >> cmd;

        if (cmd == "q" || cmd == "quit" || cmd == "exit") {
            std::cout << Color::DIM << "\n  Goodbye!\n" << Color::RESET;
            break;
        } else if (cmd == "f" || cmd == "flip") {
            int r, c;
            if (iss >> r >> c) {
                if (grid.isFlippable(r, c)) {
                    grid.flip(r, c);
                } else {
                    std::cout << Color::CYAN << "  Invalid cell (" << r << "," << c << "). "
                              << "Use row [0-" << grid.rows-1 << "], col [0-" << grid.cols-1 << "]"
                              << Color::RESET << "\n";
                    std::cout << "  Press Enter to continue...";
                    std::getline(std::cin, line);
                }
            } else {
                std::cout << Color::CYAN << "  Usage: f <row> <col>" << Color::RESET << "\n";
                std::cout << "  Press Enter to continue...";
                std::getline(std::cin, line);
            }
        } else if (cmd == "r" || cmd == "random") {
            auto flipped = grid.randomFlip();
            (void)flipped;
        } else if (cmd == "m" || cmd == "mcmc") {
            int steps = 100;
            iss >> steps;
            
            std::cout << Color::WHITE << "\n  Running " << steps << " MCMC steps...\n" << Color::RESET;
            auto result = runMCMC(grid, steps, temperature);
            printMCMCStats(result);
            std::cout << "  Press Enter to continue...";
            std::getline(std::cin, line);
        } else if (cmd == "p" || cmd == "path") {
            auto currentPath = PathFinder::findShortestPath(grid, grid.startVertex(), grid.endVertex());
            std::cout << Color::WHITE << "\n  Current path: " << Color::RESET;
            for (size_t i = 0; i < currentPath.size(); ++i) {
                auto pos = grid.vertexPos(currentPath[i]);
                std::cout << Color::BRIGHT_BLUE << "(" << pos.first << "," << pos.second << ")" << Color::RESET;
                if (i < currentPath.size() - 1) std::cout << " -> ";
            }
            std::cout << "\n  Path length: " << Color::BRIGHT_BLUE << Color::BOLD 
                      << currentPath.size() - 1 << Color::RESET << " edges\n\n";
            std::cout << "  Press Enter to continue...";
            std::getline(std::cin, line);
        } else if (cmd == "s" || cmd == "stats") {
            auto currentPath = PathFinder::findShortestPath(grid, grid.startVertex(), grid.endVertex());
            std::cout << Color::WHITE << "\n  Grid: " << grid.rows << "x" << grid.cols 
                      << " cells (" << grid.vertexCount() << " vertices, "
                      << grid.triangleCount() << " triangles)\n";
            std::cout << "  Total flips performed: " << grid.flipCount << "\n";
            std::cout << "  State space size: 2^" << (grid.rows * grid.cols) 
                      << " = " << (1 << (grid.rows * grid.cols)) << " triangulations\n";
            std::cout << "  Path length: " << currentPath.size() - 1 << " edges\n\n";
            std::cout << Color::RESET;
            std::cout << "  Press Enter to continue...";
            std::getline(std::cin, line);
        } else if (cmd == "reset") {
            grid.reset();
        } else if (cmd == "t" || cmd == "temp") {
            double newTemp;
            if (iss >> newTemp) {
                temperature = newTemp;
                std::cout << Color::WHITE << "  Temperature set to " << temperature << "\n" << Color::RESET;
            } else {
                std::cout << Color::WHITE << "  Current temperature: " << temperature << "\n" << Color::RESET;
            }
            std::cout << "  Press Enter to continue...";
            std::getline(std::cin, line);
        } else if (!cmd.empty()) {
            std::cout << Color::DIM << "  Unknown command: " << cmd << "\n" << Color::RESET;
            std::cout << "  Press Enter to continue...";
            std::getline(std::cin, line);
        }
    }

    return 0;
}