planetsim/test_logic.cpp

// Headless logic test for Phase 1 tectonics. No display / raylib needed.
//
//   g++ -std=c++17 -O2 -Isrc/sim test_logic.cpp src/sim/IcoSphere.cpp \
//       src/sim/Planet.cpp src/sim/PlanetTectonics.cpp src/sim/PlanetDrift.cpp \
//       src/sim/PlanetErosion.cpp src/sim/PlanetHydrology.cpp src/sim/PlanetIO.cpp \
//       -o /tmp/t && /tmp/t
//
// Verifies the invariants documented in CLAUDE.md so Planet::step() and the
// icosphere can be changed with confidence without launching the window.

#include "Planet.hpp"
#include "Projection.hpp"
#include <cstdio>
#include <map>
#include <cmath>
#include <algorithm>

static int failures = 0;

static double angBetweenDeg(const Vec3& a, const Vec3& b) {
    return std::acos(std::clamp(a.dot(b), -1.0, 1.0)) * 180.0 / M_PI;
}

static void check(bool cond, const char* what) {
    std::printf("  [%s] %s\n", cond ? "PASS" : "FAIL", what);
    if (!cond) ++failures;
}

// Euler characteristic for the icosphere: V - E + F == 2 (sphere topology).
static bool eulerOk(const Planet& p, size_t vertCount) {
    const std::vector<int>& tri = p.triIndices();
    size_t F = tri.size() / 3;
    // Each interior edge is shared by exactly 2 triangles; count unique edges.
    std::map<std::pair<int,int>, int> edges;
    for (size_t k = 0; k + 2 < tri.size(); k += 3) {
        int a = tri[k], b = tri[k + 1], c = tri[k + 2];
        int e[3][2] = {{a, b}, {b, c}, {c, a}};
        for (auto& pr : e) {
            int lo = std::min(pr[0], pr[1]), hi = std::max(pr[0], pr[1]);
            edges[{lo, hi}]++;
        }
    }
    size_t E = edges.size();
    long long euler = (long long)vertCount - (long long)E + (long long)F;
    std::printf("    V=%zu E=%zu F=%zu  V-E+F=%lld\n", vertCount, E, F, euler);
    return euler == 2;
}

int main() {
    // ---- Geometry: icosphere level 5 -------------------------------------
    Planet planet;
    PlanetConfig cfg;
    cfg.subdivisions = 5;
    cfg.seed = 1337;
    planet.generate(cfg);

    std::printf("Geometry (level 5):\n");
    check(planet.cells.size() == 10242, "10242 cells at subdivision 5");
    check(eulerOk(planet, planet.cells.size()), "Euler characteristic V-E+F == 2");

    // Vertex degrees: exactly 12 should have degree 5, the rest degree 6.
    int deg5 = 0, deg6 = 0, other = 0;
    for (auto& c : planet.cells) {
        if (c.neighbors.size() == 5) ++deg5;
        else if (c.neighbors.size() == 6) ++deg6;
        else ++other;
    }
    std::printf("    degree5=%d degree6=%d other=%d\n", deg5, deg6, other);
    check(deg5 == 12, "exactly 12 degree-5 vertices");
    check(other == 0, "all remaining vertices are degree 6");

    double cw = planet.cellWidthMeters() / 1000.0;
    std::printf("    cell width ~ %.0f km\n", cw);
    check(cw > 200.0 && cw < 250.0, "cell width ~223 km");

    // ---- Plates ----------------------------------------------------------
    bool allAssigned = true;
    for (auto& c : planet.cells)
        if (c.plateId < 0 || c.plateId >= cfg.plateCount) allAssigned = false;
    check(allAssigned, "every cell assigned to a valid plate");
    check((int)planet.plates.size() == cfg.plateCount, "plateCount plates created");

    // ---- Tectonics: run ~40 ticks ----------------------------------------
    for (int i = 0; i < 40; ++i) planet.step();

    double lo = planet.minElevation(), hi = planet.maxElevation();
    std::printf("Tectonics after 40 ticks:\n");
    std::printf("    elevation range %.0f .. %.0f m\n", lo, hi);
    check(hi > 2000.0, "clear mountains form (max > 2000 m)");
    check(lo < -6000.0, "deep trenches form (min < -6000 m)");

    // Relief must be GRADED, not saturated to the clamp rails (the old bug:
    // uplift so strong every boundary cell railed to +/- the clamp in 1 tick).
    int pinned = 0, midLand = 0, midSea = 0;
    for (auto& c : planet.cells) {
        if (c.elevation >= 8999.0 || c.elevation <= -10999.0) ++pinned;
        if (c.elevation > 800.0 && c.elevation < 2000.0)      ++midLand; // belt flanks
        if (c.elevation < -4500.0 && c.elevation > -6000.0)   ++midSea;  // trench flanks
    }
    double pinnedPct = 100.0 * pinned / planet.cells.size();
    std::printf("    pinned to clamp: %d (%.2f%%)  flank cells: land=%d sea=%d\n",
                pinned, pinnedPct, midLand, midSea);
    check(pinnedPct < 2.0, "not saturated: <2% of cells pinned to clamp rails");
    check(midLand > 0 && midSea > 0, "graded relief: mountains/trenches have flanks");

    bool finite = true;
    for (auto& c : planet.cells)
        if (!std::isfinite(c.elevation)) finite = false;
    check(finite, "no NaN/Inf elevations (numerically stable)");

    // ---- Determinism: same seed -> identical result ----------------------
    Planet p2;
    p2.generate(cfg);
    for (int i = 0; i < 40; ++i) p2.step();
    bool identical = (p2.cells.size() == planet.cells.size());
    for (size_t i = 0; identical && i < p2.cells.size(); ++i)
        if (p2.cells[i].elevation != planet.cells[i].elevation) identical = false;
    check(identical, "same seed reproduces identical world (deterministic)");

    // ---- Equal Earth projection round-trip (used by the 2D map + hover) ---
    {
        double maxErr = 0.0; int rejected = 0;
        for (const auto& c : planet.cells) {
            double lon, lat; dirToLonLat(c.unit, lon, lat);
            double x, y; EqualEarth::forward(lon, lat, x, y);
            double lo2, la2;
            if (!EqualEarth::inverse(x, y, lo2, la2)) { ++rejected; continue; }
            Vec3 d2 = lonLatToDir(lo2, la2);
            maxErr = std::max(maxErr, angBetweenDeg(c.unit, d2));
        }
        std::printf("    Equal Earth round-trip: max err %.2e deg, rejected %d\n", maxErr, rejected);
        check(maxErr < 1e-3 && rejected == 0, "Equal Earth forward/inverse round-trips");
    }

    // ---- Subgrid: continuity at center + neighbor coverage ----------------
    {
        int cell = (int)planet.cells.size() / 2;
        auto sg = planet.makeSubGrid(cell, 16);
        bool sgFinite = true; int ownCell = 0, ownNbr = 0;
        for (auto& s : sg->sub) {
            if (!std::isfinite(s.elevation)) sgFinite = false;
            if (s.nearestMacro == cell) ++ownCell; else ++ownNbr;
        }
        const SubCell& ctr = sg->sub[(16 / 2) * 16 + 16 / 2];
        double centerErr = std::fabs(ctr.elevation - planet.cells[cell].elevation);
        std::printf("    subgrid: %zu subcells, center |diff| %.0f m, own %d / nbr %d\n",
                    sg->sub.size(), centerErr, ownCell, ownNbr);
        check(sgFinite && sg->sub.size() == 256 && centerErr < 800.0 && ownCell > 0 && ownNbr > 0,
              "subgrid is finite, continuous at center, and reaches neighbors");
    }

    std::printf("\n%s (%d failure%s)\n",
                failures == 0 ? "ALL TESTS PASSED" : "TESTS FAILED",
                failures, failures == 1 ? "" : "s");
    return failures == 0 ? 0 : 1;
}