settled-reach/server/tests/hydrology_equilibrium_bench.rs

//! Equilibrium hydrology solver benchmarks (T-1177, body-map-viewer workshop
//! measurement ①).
//!
//! Measures `hydrology_equilibrium::solve` at today's Layer-1 working grid
//! (512×256) and at two 4K-class synthetic grids (~768×432 ≈ 330K cells,
//! matching the workshop's per-gridunit derive measurement ②'s canvas size
//! for direct comparison; 3840×2160 ≈ 8.3M cells, the "computer catches
//! fire" ceiling case). Real GJ1c heightmap data is used at 512×256 (the
//! actual production working-grid size — no upsampling needed there); the
//! two larger grids use synthetic elevation (documented in
//! `synthetic_elevation` below) since no committed heightmap PNG is stored
//! at those resolutions and generating/committing new fixture PNGs is out of
//! scope for a measurement prototype.
//!
//! Run: `cargo test --release --test hydrology_equilibrium_bench -- --ignored --nocapture`
//! (debug numbers are not representative — this crate's other benches use
//! the same release-only convention).
//!
//! Hardware: 16 cores, Rayon default thread pool (14 workers observed
//! elsewhere in this repo's benches on the same machine).

use std::time::Instant;

use settled_reach_server::atlas::heightmap::load_heightmap_png;
use settled_reach_server::atlas::hydrology_equilibrium::{solve, ClimateInputs};

/// Deterministic synthetic elevation for grids larger than any committed
/// heightmap PNG. NOT a real body — a smooth multi-octave ridged surface
/// (a few sine terms at different frequencies/phases, summed and
/// normalized) chosen to produce a realistic MIX of basins or the solver
/// would have nothing to fill: a plain gradient (as `zoom_ladder_bench.rs`'s
/// `bench_hm` uses for its unrelated per-cell derive cost) has almost no
/// interior depressions, which would make this bench measure an
/// unrepresentative best case (priority-flood on a monotonic slope is
/// nearly free — the expensive part is basin interiors + overflow search).
/// Purely a function of `(row, col, width, height)` — the same call always
/// produces the same bytes, so the resulting elevation grid is itself
/// deterministic (D-010), even though it is synthetic rather than sourced
/// from a real body.
fn synthetic_elevation(width: u32, height: u32) -> Vec<f32> {
    let w = width as f64;
    let h = height as f64;
    let n = (width * height) as usize;
    (0..n)
        .map(|i| {
            let row = (i / width as usize) as f64;
            let col = (i % width as usize) as f64;
            let x = col / w;
            let y = row / h;
            // Several sine octaves at different frequencies/phases — enough
            // basins (local minima not at the grid boundary) that the
            // priority-flood + overflow-search work is representative, not
            // a degenerate monotonic slope.
            let v = 0.5
                + 0.25 * (x * std::f64::consts::TAU * 3.0).sin() * (y * std::f64::consts::TAU * 2.0).cos()
                + 0.15 * (x * std::f64::consts::TAU * 7.3 + 1.7).sin() * (y * std::f64::consts::TAU * 5.1).sin()
                + 0.10 * (x * std::f64::consts::TAU * 13.0).cos() * (y * std::f64::consts::TAU * 11.0 + 0.4).sin();
            v.clamp(0.0, 1.0) as f32
        })
        .collect()
}

fn gj1c_512x256() -> (Vec<f32>, f32) {
    let src = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"))
        .join("../wiki/star-systems/GJ-1/bodies/GJ1c/heightmap.png");
    let heightmap =
        load_heightmap_png(&src, "GJ1c", 0.3).expect("decode committed GJ1c heightmap");
    let small = heightmap.downsample(512, 256); // GRID_W x GRID_H, the real production working grid
    (small.data, small.sea_level)
}

fn default_climate() -> ClimateInputs {
    ClimateInputs { moisture_q: 55 }
}

fn run_and_report(label: &str, width: u32, height: u32, elevation: &[f32], sea_level: f32) {
    let n_cells = (width as u64) * (height as u64);

    // Cold run.
    let t0 = Instant::now();
    let result_cold = solve(elevation, width, height, sea_level, default_climate());
    let cold = t0.elapsed();

    // Warm run (same process, allocator/cache warm — same input).
    let t1 = Instant::now();
    let result_warm = solve(elevation, width, height, sea_level, default_climate());
    let warm = t1.elapsed();

    let lake_cells: usize = result_cold
        .basins
        .iter()
        .map(|b| b.cells.len())
        .sum();
    let carved_cells = result_cold.cliff_edge.iter().filter(|&&c| c).count();
    let endorheic_count = result_cold
        .basins
        .iter()
        .filter(|b| {
            !b.cells.is_empty()
                && matches!(
                    b.outcome,
                    settled_reach_server::atlas::hydrology_equilibrium::BasinOutcome::Endorheic { .. }
                )
        })
        .count();
    let overflow_count = result_cold
        .basins
        .iter()
        .filter(|b| {
            !b.cells.is_empty()
                && matches!(
                    b.outcome,
                    settled_reach_server::atlas::hydrology_equilibrium::BasinOutcome::Overflow { .. }
                )
        })
        .count();

    println!("\n=== {label} ({width}x{height} = {n_cells} cells) ===");
    println!(
        "  cold: {:>9.2} ms total, {:>8.1} ns/cell",
        cold.as_secs_f64() * 1000.0,
        cold.as_secs_f64() * 1e9 / n_cells as f64
    );
    println!(
        "  warm: {:>9.2} ms total, {:>8.1} ns/cell",
        warm.as_secs_f64() * 1000.0,
        warm.as_secs_f64() * 1e9 / n_cells as f64
    );
    println!(
        "  basins: {} total ({} overflow, {} endorheic, {} empty/no-depression), \
         lake cells: {lake_cells}, carved gorge cells: {carved_cells}",
        result_cold.basins.len(),
        overflow_count,
        endorheic_count,
        result_cold.basins.len() - overflow_count - endorheic_count,
    );
    std::hint::black_box(&result_warm);
}

#[test]
#[ignore]
fn bench_512x256_real_gj1c() {
    let (elev, sea_level) = gj1c_512x256();
    run_and_report("512x256 (real GJ1c, production working-grid size)", 512, 256, &elev, sea_level);
}

#[test]
#[ignore]
fn bench_768x432_synthetic() {
    let (w, h) = (768u32, 432u32);
    let elev = synthetic_elevation(w, h);
    run_and_report(
        "768x432 (~330K cells, 4K-class synthetic — see synthetic_elevation doc)",
        w,
        h,
        &elev,
        0.35,
    );
}

#[test]
#[ignore]
fn bench_3840x2160_synthetic() {
    let (w, h) = (3840u32, 2160u32);
    let elev = synthetic_elevation(w, h);
    run_and_report(
        "3840x2160 (~8.3M cells, 4K synthetic — see synthetic_elevation doc)",
        w,
        h,
        &elev,
        0.35,
    );
}

/// Determinism proof at bench scale (T-1177 mandatory deliverable): same
/// seed + input → byte-identical solver output, twice, on a non-trivial
/// grid (not just the small fixtures already covered by the module's own
/// unit tests).
#[test]
#[ignore]
fn determinism_at_330k_cells() {
    let (w, h) = (768u32, 432u32);
    let elev = synthetic_elevation(w, h);
    let r1 = solve(&elev, w, h, 0.35, default_climate());
    let r2 = solve(&elev, w, h, 0.35, default_climate());
    assert_eq!(r1.filled_scaled, r2.filled_scaled, "filled surface must be byte-identical");
    assert_eq!(
        r1.channel_depth_scaled, r2.channel_depth_scaled,
        "carved channel depth must be byte-identical"
    );
    assert_eq!(r1.cliff_edge, r2.cliff_edge, "cliff-edge flags must be byte-identical");
    assert_eq!(r1.basins.len(), r2.basins.len(), "basin count must be identical");
    for (a, b) in r1.basins.iter().zip(r2.basins.iter()) {
        assert_eq!(a.basin_id, b.basin_id);
        assert_eq!(a.cells, b.cells);
        assert_eq!(a.spill_level_scaled, b.spill_level_scaled);
        assert_eq!(a.spill_cell, b.spill_cell);
        assert_eq!(format!("{:?}", a.outcome), format!("{:?}", b.outcome));
    }
    println!(
        "\n=== determinism proof (768x432, {} basins) — byte-identical across two solves ===",
        r1.basins.len()
    );
}

/// Rayon-parallel throughput: the REAL production shape is N independent
/// bodies, each solved once (not one body's solve parallelized internally —
/// priority-flood's heap and the overflow Dijkstra search are both globally
/// sequential by nature, same as `road_graph.rs`'s A*). This measures what
/// "always keep hydrology for ~273 bodies" would cost in wall-clock if
/// solved across the Rayon pool, at the 512×256 production grid size —
/// directly answering the workshop's red-flag-2-adjacent question of
/// whether per-body-open hydrology is affordable at scale.
#[test]
#[ignore]
fn bench_parallel_273_bodies_at_512x256() {
    use rayon::prelude::*;

    let (elev, sea_level) = gj1c_512x256();
    let body_count = 273usize;

    let t0 = Instant::now();
    let total_basins: usize = (0..body_count)
        .into_par_iter()
        .map(|_| {
            let result = solve(&elev, 512, 256, sea_level, default_climate());
            result.basins.len()
        })
        .sum();
    let elapsed = t0.elapsed();

    println!(
        "\n=== {body_count} bodies x 512x256, Rayon par_iter ({} threads available) ===",
        std::thread::available_parallelism().map(|n| n.get()).unwrap_or(0)
    );
    println!(
        "  {:>9.2} ms total, {:>7.2} ms/body average, {total_basins} basins summed",
        elapsed.as_secs_f64() * 1000.0,
        elapsed.as_secs_f64() * 1000.0 / body_count as f64
    );
}