main.cpp

Go to the documentation of this file.

/// @file main.cpp
/// @brief 2D Ising model  -- classic mode + nucleation with umbrella sampling
///
/// Absorbs ~/projects/IsingNucleation into the numerics project.
/// Exact match to IsingNucleation parameters:
///   n=300, J=1, beta=0.58 (T~=1.724), F=0.1 (nucleation), F=0 (classic)
///   Random spin selection per sweep  (n^2 picks with replacement)
///   DeltaE = 2J*s*(Sigma nbrs) - 2F*s  (same sign convention as mc_step_boltz.c)
///
/// Optimizations vs. the original C implementation:
///   1. Boltzmann lookup table    -- DeltaE discrete -> precomputed exp() table, no
///                                  runtime exp() in the sweep hot path
///   2. Precomputed neighbor arrays  -- up/dn/lt/rt[NN] avoid repeated modulo
///                                    arithmetic per spin selection
///   3. Array-based BFS cluster detection  -- iterative, pre-allocated flat
///                                          queue; no recursion stack risk,
///                                          no heap allocations per call
///   4. SparseMatrix adj retained for energy_per_spin() (SIMD sparse_matvec)
///
/// Two modes (TAB to switch):
///   Classic     -- free Metropolis; magnetisation and E/N^2 graphs
///   Nucleation  -- umbrella sampling; tracks largest spin-down cluster (nucleus)
 
#include "numerics.hpp"
#include "stats/stats.hpp"          // num::RunningStats, num::Histogram, num::autocorr_time
#include "analysis/roots.hpp"       // num::newton  (for mean-field solver)
#include "stochastic/markov.hpp"    // num::markov::metropolis_sweep_prob, umbrella_sweep_prob
#include "spatial/pbc_lattice.hpp"  // num::PBCLattice2D
#include "spatial/connected_components.hpp"  // num::ClusterResult, connected_components
#include <raylib.h>
#include <cmath>
#include <random>
#include <vector>
#include <algorithm>
#include <cstring>
 
using namespace num;
 
static constexpr int  N    = 300;             // matches IsingNucleation/main.c
static constexpr int  NN   = N * N;
static constexpr real J    = 1.0;
static const     real TC   = 2.26918531421;   // Onsager exact
 
static constexpr int  WIN_W  = 1080;
static constexpr int  WIN_H  = 720;
static constexpr int  SIM_SZ = 680;
static constexpr int  PAN_X  = SIM_SZ + 20;
static constexpr int  PAN_W  = WIN_W - PAN_X - 10;
 
// BFS cluster detection uses num::ClusterResult + num::connected_components.
 
struct RingBuffer {
    std::vector<float> data;
    int head = 0, count = 0;
 
    explicit RingBuffer(int cap) : data(cap, 0.f) {}
 
    void push(float v) {
        data[head] = v;
        head  = (head + 1) % cap();
        count = std::min(count + 1, cap());
    }
 
    float get(int i) const { return data[(head - count + i + cap() * 2) % cap()]; }
    float back()     const { return count ? get(count - 1) : 0.f; }
    int   size()     const { return count; }
    int   cap()      const { return static_cast<int>(data.size()); }
};
 
struct IsingLattice {
    Vector       spins;       // NN spins, values +/-1.0
    Vector       ones;        // NN ones  (for dot-product magnetisation)
    Vector       nbr_buf;     // scratch buffer for energy_per_spin (sparse_matvec)
    SparseMatrix adj;         // 4-neighbor periodic CSR (used only for energy)
 
    // Precomputed PBC neighbor indices  -- avoids modulo in hot paths
    num::PBCLattice2D nbr_;
 
    // Boltzmann acceptance table.
    // DeltaE = 2J*s*nbr_sum - 2F*s,   sin{+/-1},  nbr_sumin{-4,-2,0,2,4}
    // boltz[si][ni] = min(1, exp(-beta*DeltaE))
    //   si: 0 -> s=-1,  1 -> s=+1
    //   ni: (nbr_sum/2 + 2) in {0..4}
    real boltz[2][5]{};
 
    std::vector<real> save_buf;   // pre-allocated for sweep_umbrella
 
    std::mt19937 rng;
    std::uniform_real_distribution<real>  dist{0.0, 1.0};   // used in random_init()
 
    real beta = 1.0 / (1.0 / 0.58);   // beta=0.58, T~=1.724  (matches IsingNucleation/main.c)
    real F    = 0.0;
 
    IsingLattice()
        : spins(NN, 1.0)
        , ones(NN, 1.0)
        , nbr_buf(NN, 0.0)
        , adj(build_adj())
        , nbr_(N)
        , save_buf(NN)
        , rng(std::random_device{}())   // non-deterministic seed, like srand(time(0))
    {
        rebuild_boltz();
        random_init();
    }
 
    void set_temperature(real T) { beta = (T > 1e-6) ? 1.0/T : 1e6; rebuild_boltz(); }
    void set_field(real Fnew)    { F = Fnew; rebuild_boltz(); }
    real temperature() const     { return (beta > 1e-6) ? 1.0/beta : 1e6; }
 
    void rebuild_boltz() {
        for (int si = 0; si < 2; ++si) {
            real s = (si == 0) ? -1.0 : 1.0;
            for (int ni = 0; ni < 5; ++ni) {
                real ns = -4.0 + 2.0 * ni;            // nbr_sum in {-4,-2,0,2,4}
                real dE = 2.0 * J * s * ns - 2.0 * F * s;
                boltz[si][ni] = num::markov::boltzmann_accept(dE, beta);
            }
        }
    }
 
    void random_init() {
        for (int i = 0; i < NN; ++i)
            spins[i] = (dist(rng) < 0.5) ? 1.0 : -1.0;
    }
 
    void all_up() { for (int i = 0; i < NN; ++i) spins[i] = 1.0; }
 
    // Square nucleus of spin-down spins at centre, area ~= r (as in IsingNucleation).
    void generate_nucleus(real r = 100.0) {
        all_up();
        int half = static_cast<int>(std::sqrt(r) * 0.5);
        int cx   = N / 2;
        for (int row = cx - half; row < cx + half; ++row)
            for (int col = cx - half; col < cx + half; ++col)
                spins[row * N + col] = -1.0;
    }
 
    // Standard Metropolis sweep: n^2 random spin picks (with replacement).
    // Matches mc_step_boltz() in IsingNucleation/src/monte_carlo.c exactly.
    // Delegates to num::markov::metropolis_sweep_prob  -- Boltzmann table
    // supplies the acceptance probability so no exp() runs at runtime.
    void sweep() {
        num::markov::metropolis_sweep_prob(
            static_cast<idx>(NN),
            [&](idx i) {
                real s  = spins[i];
                real ns = spins[nbr_.up[i]] + spins[nbr_.dn[i]]
                        + spins[nbr_.lt[i]] + spins[nbr_.rt[i]];
                int  si = (s < 0.0) ? 0 : 1;
                int  ni = static_cast<int>(ns * 0.5 + 2.0);
                return boltz[si][ni];
            },
            [&](idx i) { spins[i] = -spins[i]; },
            rng);
    }
 
    // Umbrella sampling sweep: revert full state if nucleus size leaves [lo, hi].
    // Delegates to num::markov::umbrella_sweep_prob; BFS cluster detection is
    // passed as the measure_order callable so det.id is filled for rendering.
    int sweep_umbrella(int lo, int hi, num::ClusterResult& det) {
        auto result = num::markov::umbrella_sweep_prob(
            static_cast<idx>(NN),
            [&](idx i) {
                real s  = spins[i];
                real ns = spins[nbr_.up[i]] + spins[nbr_.dn[i]]
                        + spins[nbr_.lt[i]] + spins[nbr_.rt[i]];
                int  si = (s < 0.0) ? 0 : 1;
                int  ni = static_cast<int>(ns * 0.5 + 2.0);
                return boltz[si][ni];
            },
            [&](idx i) { spins[i] = -spins[i]; },
            [&]() { for (int i = 0; i < NN; ++i) save_buf[i] = spins[i]; },
            [&]() { for (int i = 0; i < NN; ++i) spins[i] = save_buf[i]; },
            [&]() -> idx {
                det = num::connected_components(NN,
                    [&](int i) { return spins[i] < 0.0; },
                    [&](int i, auto&& visit) {
                        visit(nbr_.up[i]); visit(nbr_.dn[i]);
                        visit(nbr_.lt[i]); visit(nbr_.rt[i]);
                    });
                return static_cast<idx>(det.largest_size);
            },
            num::markov::UmbrellaWindow{static_cast<idx>(lo), static_cast<idx>(hi)},
            rng);
        return static_cast<int>(result.order_param);
    }
 
    // |<m>| = |Sigma s_i| / N^2  via SIMD dot product
    real magnetization() const {
        return std::abs(dot(spins, ones, Backend::seq)) / static_cast<real>(NN);
    }
 
    // <E>/N^2 using SparseMatrix adj and SIMD sparse_matvec + dot
    real energy_per_spin() {
        sparse_matvec(adj, spins, nbr_buf);
        real coup = -J * dot(spins, nbr_buf, Backend::seq) / (4.0 * NN);
        real fld  = -F * dot(spins, ones,    Backend::seq) / static_cast<real>(NN);
        return coup + fld;
    }
 
    // Mean-field order parameter: solve m = tanh(beta*4J*m + beta*F) via num::newton().
    static real mean_field_m(real beta_val, real F_val) {
        auto f  = [=](real m) {
            real arg = 4.0 * beta_val * J * m + beta_val * F_val;
            return std::tanh(arg) - m;
        };
        auto df = [=](real m) {
            real arg = 4.0 * beta_val * J * m + beta_val * F_val;
            real th  = std::tanh(arg);
            return 4.0 * beta_val * J * (1.0 - th * th) - 1.0;
        };
        auto res = newton(f, df, 0.9);
        return std::max(-1.0, std::min(1.0, res.root));
    }
 
private:
    static SparseMatrix build_adj() {
        std::vector<idx>  rows_t, cols_t;
        std::vector<real> vals_t;
        rows_t.reserve(4 * NN);
        cols_t.reserve(4 * NN);
        vals_t.reserve(4 * NN);
        for (int row = 0; row < N; ++row) {
            for (int col = 0; col < N; ++col) {
                int i = row * N + col;
                int nbrs[4] = {
                    ((row-1+N)%N)*N + col,  ((row+1)%N)*N   + col,
                    row*N + (col-1+N)%N,    row*N + (col+1)%N
                };
                for (int nb : nbrs) {
                    rows_t.push_back(static_cast<idx>(i));
                    cols_t.push_back(static_cast<idx>(nb));
                    vals_t.push_back(1.0);
                }
            }
        }
        return SparseMatrix::from_triplets(NN, NN, rows_t, cols_t, vals_t);
    }
};
 
static void draw_graph(const RingBuffer& buf, int x, int y, int w, int h,
                       float y_lo, float y_hi, Color col, const char* label,
                       float ref_val = -999.f) {
    DrawRectangle(x, y, w, h, {20, 20, 20, 255});
    DrawRectangleLines(x, y, w, h, GRAY);
    DrawText(label, x + 4, y + 4, 11, LIGHTGRAY);
    if (buf.size() < 2) return;
 
    float range = y_hi - y_lo;
    if (range < 1e-6f) range = 1.f;
 
    if (ref_val > -998.f) {
        int ry = y + h - static_cast<int>((ref_val - y_lo) / range * h);
        DrawLine(x, std::clamp(ry,y,y+h), x+w, std::clamp(ry,y,y+h), {255,220,0,160});
    }
 
    for (int i = 1; i < buf.size(); ++i) {
        float x0f = static_cast<float>(i-1) / (buf.cap()-1) * w;
        float x1f = static_cast<float>(i)   / (buf.cap()-1) * w;
        int py0 = y+h - static_cast<int>((buf.get(i-1)-y_lo)/range*h);
        int py1 = y+h - static_cast<int>((buf.get(i)  -y_lo)/range*h);
        DrawLine(x+static_cast<int>(x0f), std::clamp(py0,y,y+h),
                 x+static_cast<int>(x1f), std::clamp(py1,y,y+h), col);
    }
}
 
static real draw_slider(int x, int y, int w, real val, real lo, real hi,
                        float ref_frac, const char* fmt, Color knob_col) {
    DrawRectangle(x, y-3, w, 6, DARKGRAY);
    if (ref_frac >= 0.f) {
        int rx = x + static_cast<int>(ref_frac * w);
        DrawLine(rx, y-9, rx, y+9, YELLOW);
    }
    int kx = x + static_cast<int>(((val-lo)/(hi-lo)) * w);
    DrawCircle(kx, y, 8, knob_col);
    DrawText(TextFormat(fmt, static_cast<float>(val)), x, y-18, 12, LIGHTGRAY);
 
    if (IsMouseButtonDown(MOUSE_BUTTON_LEFT)) {
        Vector2 mp = GetMousePosition();
        if (mp.x >= x && mp.x <= x+w && mp.y >= y-12 && mp.y <= y+12) {
            float frac = (mp.x - x) / static_cast<float>(w);
            val = lo + frac*(hi-lo);
            val = std::max(lo, std::min(hi, val));
        }
    }
    return val;
}
 
int main() {
    InitWindow(WIN_W, WIN_H, "2D Ising  -- Classic & Nucleation  [TAB to switch]");
    SetTargetFPS(60);
 
    IsingLattice lat;
    num::ClusterResult det;
    det.id.assign(NN, -2);   // default: all excluded (safe before first sweep)
    real T = 1.0 / 0.58;    // beta=0.58 from IsingNucleation/main.c
    lat.set_temperature(T);
 
    std::vector<unsigned char> pixels(NN * 4, 255);
    Image     img = {pixels.data(), N, N, 1, PIXELFORMAT_UNCOMPRESSED_R8G8B8A8};
    Texture2D tex = LoadTextureFromImage(img);
 
    RingBuffer mag_hist(300);
    RingBuffer eng_hist(300);
    RingBuffer nuc_hist(300);
 
    int  sweeps_per_frame = 5;
    bool paused = false;
    int  frame  = 0;
 
    enum class Mode { Classic, Nucleation } mode = Mode::Classic;
 
    // Nucleation state  -- matches run_sim_nucleation defaults:
    //   N_window=15, window_size=30 -> first window [0, 30], centre=15
    int nuc_target = 15;     // umbrella window centre (spins)
    int nuc_half   = 15;     // half-width = window_size/2 = 30/2
    int nuc_size   = 0;
 
    // Online statistics for nucleation observables (num::RunningStats, num::Histogram)
    RunningStats nuc_stats;
    RunningStats mag_stats;
    Histogram    nuc_histo(50, 0.0, 500.0);   // nucleus size distribution
 
    while (!WindowShouldClose()) {
        if (IsKeyPressed(KEY_SPACE)) paused = !paused;
        if (IsKeyPressed(KEY_R))     lat.random_init();
        if (IsKeyPressed(KEY_H))     lat.all_up();
        if (IsKeyPressed(KEY_N))     lat.generate_nucleus(static_cast<real>(nuc_target));
        if (IsKeyPressed(KEY_UP)   && sweeps_per_frame < 20) sweeps_per_frame++;
        if (IsKeyPressed(KEY_DOWN) && sweeps_per_frame >  1) sweeps_per_frame--;
        if (IsKeyPressed(KEY_TAB)) {
            if (mode == Mode::Classic) {
                mode = Mode::Nucleation;
                nuc_target = 15; nuc_half = 15;   // reset to run_sim_nucleation defaults
                lat.set_field(0.1);                // F=0.1  -- nucleation field
                lat.generate_nucleus(static_cast<real>(nuc_target));
                nuc_hist = RingBuffer(300);
                mag_hist = RingBuffer(300);
                nuc_stats.reset(); mag_stats.reset(); nuc_histo.reset();
            } else {
                mode = Mode::Classic;
                lat.set_field(0.0);
                lat.random_init();
                mag_hist = RingBuffer(300);
                eng_hist = RingBuffer(300);
            }
        }
 
        if (!paused) {
            if (mode == Mode::Classic) {
                for (int s = 0; s < sweeps_per_frame; ++s) lat.sweep();
                ++frame;
                if (frame % 3 == 0) {
                    mag_hist.push(static_cast<float>(lat.magnetization()));
                    eng_hist.push(static_cast<float>(lat.energy_per_spin()));
                }
            } else {
                int lo = nuc_target - nuc_half;
                int hi = nuc_target + nuc_half;
                for (int s = 0; s < sweeps_per_frame; ++s)
                    nuc_size = lat.sweep_umbrella(lo, hi, det);
                ++frame;
                if (frame % 3 == 0) {
                    nuc_hist.push(static_cast<float>(nuc_size));
                    mag_hist.push(static_cast<float>(lat.magnetization()));
                    nuc_stats.update(static_cast<real>(nuc_size));
                    mag_stats.update(lat.magnetization());
                    nuc_histo.fill(static_cast<real>(nuc_size));
                }
            }
        }
 
        if (mode == Mode::Classic) {
            for (int i = 0; i < NN; ++i) {
                unsigned char c = (lat.spins[i] > 0.0) ? 240 : 15;
                pixels[4*i+0] = pixels[4*i+1] = pixels[4*i+2] = c;
                pixels[4*i+3] = 255;
            }
        } else {
            // White = spin+1,  Red = nucleus (largest cluster),  Dark = other spin-1
            for (int i = 0; i < NN; ++i) {
                if (lat.spins[i] > 0.0) {
                    pixels[4*i+0] = 240; pixels[4*i+1] = 240; pixels[4*i+2] = 240;
                } else if (det.id[i] == det.largest_id) {
                    pixels[4*i+0] = 220; pixels[4*i+1] = 50;  pixels[4*i+2] = 50;
                } else {
                    pixels[4*i+0] = 15;  pixels[4*i+1] = 15;  pixels[4*i+2] = 15;
                }
                pixels[4*i+3] = 255;
            }
        }
        UpdateTexture(tex, pixels.data());
 
        BeginDrawing();
        ClearBackground({10, 10, 10, 255});
 
        DrawTexturePro(tex,
            {0, 0, (float)N, (float)N},
            {0, 0, (float)SIM_SZ, (float)SIM_SZ},
            {0, 0}, 0.f, WHITE);
        DrawRectangleLines(0, 0, SIM_SZ, SIM_SZ, DARKGRAY);
 
        int px = PAN_X, cy = 10;
 
        DrawText("2D Ising Model", px, cy, 17, WHITE);  cy += 22;
 
        // Mode badge
        bool is_nuc = (mode == Mode::Nucleation);
        Color mode_col = is_nuc ? Color{220,110,50,255} : SKYBLUE;
        DrawText(is_nuc ? "NUCLEATION" : "CLASSIC", px, cy, 13, mode_col);
        DrawText("[TAB]", px + 90, cy, 10, DARKGRAY);  cy += 18;
 
        DrawText(TextFormat("N=%d  spd=%d", N, sweeps_per_frame), px, cy, 11, GRAY);
        cy += 20;
 
        if (!is_nuc) {
            cy += 10;
            float tc_frac = static_cast<float>((TC - 0.05) / (5.0 - 0.05));
            T = draw_slider(px, cy, PAN_W, T, 0.05, 5.0, tc_frac, "T = %.3f", SKYBLUE);
            lat.set_temperature(T);
            DrawText("Tc", px + static_cast<int>(tc_frac * PAN_W) - 5, cy+10, 9, YELLOW);
            cy += 30;
 
            cy += 6;
            real Fnew = draw_slider(px, cy, PAN_W, lat.F, 0.0, 0.2, -1.f,
                                    "F = %.4f", PINK);
            if (Fnew != lat.F) lat.set_field(Fnew);
            cy += 30;
 
            float cur_m = static_cast<float>(lat.magnetization());
            float mf_m  = static_cast<float>(IsingLattice::mean_field_m(lat.beta, lat.F));
            cy += 6;
            DrawText(TextFormat("|m| = %.4f", cur_m),             px, cy,      13, WHITE);
            DrawText(TextFormat("MF  = %.4f", mf_m),              px, cy + 16, 13, YELLOW);
            DrawText(TextFormat("E/N = %.4f", eng_hist.back()),   px, cy + 32, 13, LIGHTGRAY);
            DrawText(TextFormat("beta   = %.3f", (float)lat.beta),   px, cy + 48, 13, GRAY);
            cy += 68;
 
            Color pc = (T < TC) ? Color{80,200,80,255} : Color{200,80,80,255};
            DrawText((T < TC) ? "ORDERED" : "DISORDERED", px, cy, 14, pc);
            cy += 20;
 
            int gh = 128;
            draw_graph(mag_hist, px, cy, PAN_W, gh, 0.f, 1.f, GREEN, "|m|", mf_m);
            cy += gh + 10;
            draw_graph(eng_hist, px, cy, PAN_W, gh, -2.2f, 0.f, ORANGE, "E/spin");
            cy += gh + 10;
 
        } else {
            cy += 10;
            // Temperature slider (still useful to explore metastability)
            float tc_frac = static_cast<float>((TC - 0.05) / (5.0 - 0.05));
            T = draw_slider(px, cy, PAN_W, T, 0.05, 5.0, tc_frac, "T = %.3f", SKYBLUE);
            lat.set_temperature(T);
            cy += 30;
 
            // Field slider  -- drives nucleation of spin-down phase
            cy += 6;
            real Fnew = draw_slider(px, cy, PAN_W, lat.F, 0.0, 0.3, -1.f,
                                    "F = %.4f  (nucleation field)", PINK);
            if (Fnew != lat.F) lat.set_field(Fnew);
            cy += 30;
 
            // Umbrella target and half-width sliders
            cy += 6;
            {
                int old_tgt = nuc_target;
                real tgt = draw_slider(px, cy, PAN_W, static_cast<real>(nuc_target),
                                       10, 2000, -1.f, "Nuc target = %.0f", {220,110,50,255});
                nuc_target = static_cast<int>(tgt);
                if (nuc_target != old_tgt) { nuc_stats.reset(); mag_stats.reset(); nuc_histo.reset(); }
            }
            cy += 30;
 
            cy += 4;
            {
                int old_hw = nuc_half;
                real hw = draw_slider(px, cy, PAN_W, static_cast<real>(nuc_half),
                                      5, 200, -1.f, "Window +/- = %.0f", PURPLE);
                nuc_half = static_cast<int>(hw);
                if (nuc_half != old_hw) { nuc_stats.reset(); mag_stats.reset(); nuc_histo.reset(); }
            }
            cy += 30;
 
            // Current nucleus stats (uses num::RunningStats)
            int lo = nuc_target - nuc_half;
            int hi = nuc_target + nuc_half;
            cy += 6;
            DrawText(TextFormat("N_nuc = %d", nuc_size),                     px, cy,      14, {220,80,80,255});
            DrawText(TextFormat("Window [%d, %d]", lo, hi),                  px, cy + 14, 10, GRAY);
            DrawText(TextFormat("<N> = %.1f  sigma = %.1f",
                                (float)nuc_stats.mean,
                                (float)nuc_stats.std_dev()),                  px, cy + 26, 11, {220,150,80,255});
            DrawText(TextFormat("<|m|> = %.4f  sigma = %.4f",
                                (float)mag_stats.mean,
                                (float)mag_stats.std_dev()),                  px, cy + 38, 11, WHITE);
            DrawText(TextFormat("n = %llu sweeps",
                                static_cast<unsigned long long>(nuc_stats.count)),
                                                                              px, cy + 50, 10, GRAY);
            DrawText(TextFormat("beta = %.3f", (float)lat.beta),                px, cy + 62, 10, GRAY);
            cy += 76;
 
            // Nucleus size graph  -- reference line at target
            int gh = 118;
            float graph_hi = static_cast<float>(nuc_target + nuc_half + 30);
            draw_graph(nuc_hist, px, cy, PAN_W, gh, 0.f, graph_hi,
                       {220,80,80,255}, "Nucleus size  (red cluster)",
                       static_cast<float>(nuc_target));
            cy += gh + 8;
 
            // Magnetisation graph
            draw_graph(mag_hist, px, cy, PAN_W, gh, 0.f, 1.f, GREEN, "|m|");
            cy += gh + 8;
        }
 
        // Controls (bottom of panel)
        cy = WIN_H - 44;
        DrawText("SPACE pause  R rand  H up  N nucleus", px, cy,      9, DARKGRAY);
        DrawText("UP/DOWN speed  TAB mode",               px, cy + 12, 9, DARKGRAY);
 
        if (paused)
            DrawText("PAUSED", SIM_SZ/2 - 45, SIM_SZ/2 - 12, 26, {255,100,100,220});
 
        EndDrawing();
    }
 
    UnloadTexture(tex);
    CloseWindow();
    return 0;
}