tide-maxflow 0.1.0

//! Generate a pool of test graphs in DIMACS max-flow format.
//!
//! DIMACS format:
//!   c comment
//!   p max NODES ARCS
//!   n SOURCE s
//!   n SINK t
//!   a FROM TO CAPACITY
//!
//! Vertices are 1-indexed in DIMACS.

use std::fs;
use std::io::{BufWriter, Write};
use std::path::Path;

struct DimacsGraph {
    name: String,
    n: usize,
    source: usize,
    sink: usize,
    edges: Vec<(usize, usize, i64)>,
}

impl DimacsGraph {
    fn write_to(&self, dir: &Path) {
        let path = dir.join(format!("{}.max", self.name));
        let file = fs::File::create(&path).unwrap();
        let mut f = BufWriter::with_capacity(1 << 20, file); // 1 MB buffer
        writeln!(f, "c {}", self.name).unwrap();
        writeln!(f, "p max {} {}", self.n, self.edges.len()).unwrap();
        writeln!(f, "n {} s", self.source).unwrap();
        writeln!(f, "n {} t", self.sink).unwrap();
        for &(u, v, cap) in &self.edges {
            writeln!(f, "a {} {} {}", u, v, cap).unwrap();
        }
        f.flush().unwrap();
    }

    fn estimated_bytes(&self) -> usize {
        // header ~60 bytes + ~20 bytes per edge
        60 + self.edges.len() * 20
    }
}

// All generators produce 1-indexed vertices for DIMACS.

fn layered(name: &str, layers: usize, width: usize) -> DimacsGraph {
    let s = 1;
    let t = layers * width + 2;
    let n = t;
    let vid = |l: usize, p: usize| l * width + p + 2;
    let m_est = width + (layers - 1) * width * 3 + width;
    let mut edges = Vec::with_capacity(m_est);
    for p in 0..width {
        edges.push((s, vid(0, p), (50 + p % 1000) as i64));
    }
    for l in 0..layers - 1 {
        for p in 0..width {
            for d in 0..=2usize {
                edges.push((
                    vid(l, p),
                    vid(l + 1, (p + d) % width),
                    (10 + (p + d) % 1000) as i64,
                ));
            }
        }
    }
    for p in 0..width {
        edges.push((vid(layers - 1, p), t, (50 + p % 1000) as i64));
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

fn grid(name: &str, rows: usize, cols: usize) -> DimacsGraph {
    let s = 1;
    let t = rows * cols + 2;
    let n = t;
    let vid = |r: usize, c: usize| r * cols + c + 2;
    let m_est = cols + rows * (cols - 1) + (rows - 1) * cols + cols;
    let mut edges = Vec::with_capacity(m_est);
    for c in 0..cols {
        edges.push((s, vid(0, c), (50 + c % 1000) as i64));
    }
    for r in 0..rows {
        for c in 0..cols - 1 {
            edges.push((vid(r, c), vid(r, c + 1), (5 + (r + c) % 1000) as i64));
        }
    }
    for r in 0..rows - 1 {
        for c in 0..cols {
            edges.push((vid(r, c), vid(r + 1, c), (5 + (r + c) % 1000) as i64));
        }
    }
    for c in 0..cols {
        edges.push((vid(rows - 1, c), t, (50 + c % 1000) as i64));
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

fn bipartite(name: &str, left: usize, right: usize) -> DimacsGraph {
    let s = 1;
    let t = left + right + 2;
    let n = t;
    let m_est = left + left * right + right;
    let mut edges = Vec::with_capacity(m_est);
    for i in 0..left {
        edges.push((s, 2 + i, (100 + (i * 7) % 50) as i64));
    }
    for i in 0..left {
        for j in 0..right {
            edges.push((2 + i, 2 + left + j, (10 + (i * 13 + j * 7) % 40) as i64));
        }
    }
    for j in 0..right {
        edges.push((2 + left + j, t, (100 + (j * 11) % 50) as i64));
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

fn random_er(name: &str, n_inner: usize, edge_prob_pct: usize, max_cap: i64) -> DimacsGraph {
    let s = 1;
    let t = n_inner + 2;
    let n = t;
    let m_est = n_inner + (n_inner * n_inner * edge_prob_pct / 100) + n_inner;
    let mut edges = Vec::with_capacity(m_est);
    let mut lcg: u64 = 12345;

    for v in 2..=n_inner + 1 {
        lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
        let cap = (lcg % max_cap as u64) as i64 + 1;
        edges.push((s, v, cap));
    }
    for u in 2..=n_inner + 1 {
        for v in 2..=n_inner + 1 {
            if u != v {
                lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                if (lcg % 100) < edge_prob_pct as u64 {
                    lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                    let cap = (lcg % max_cap as u64) as i64 + 1;
                    edges.push((u, v, cap));
                }
            }
        }
    }
    for v in 2..=n_inner + 1 {
        lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
        let cap = (lcg % max_cap as u64) as i64 + 1;
        edges.push((v, t, cap));
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

fn chains(name: &str, k: usize, len: usize) -> DimacsGraph {
    let s = 1;
    let t = k * len + 2;
    let n = t;
    let m_est = k * (1 + (len - 1) + 1);
    let mut edges = Vec::with_capacity(m_est);
    for chain in 0..k {
        let first = 2 + chain * len;
        let last = first + len - 1;
        edges.push((s, first, (50 + (chain * 3) % 1000) as i64));
        for i in 0..len - 1 {
            edges.push((first + i, first + i + 1, (20 + (chain + i) % 1000) as i64));
        }
        edges.push((last, t, (50 + (chain * 3) % 1000) as i64));
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

/// Washington / DIMACS-style "ak" acyclic graph:
/// k layers of n_per_layer vertices, all-to-all between consecutive layers.
/// Source connects to layer 0, layer k-1 connects to sink.
fn washington(name: &str, k: usize, n_per_layer: usize, max_cap: i64) -> DimacsGraph {
    let s = 1;
    let t = k * n_per_layer + 2;
    let n = t;
    let vid = |layer: usize, idx: usize| 2 + layer * n_per_layer + idx;
    let mut edges = Vec::new();
    let mut lcg: u64 = 99991;

    // source -> layer 0
    for i in 0..n_per_layer {
        lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
        let cap = (lcg % max_cap as u64) as i64 + 1;
        edges.push((s, vid(0, i), cap));
    }
    // inter-layer: all-to-all
    for l in 0..k - 1 {
        for i in 0..n_per_layer {
            for j in 0..n_per_layer {
                lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                let cap = (lcg % max_cap as u64) as i64 + 1;
                edges.push((vid(l, i), vid(l + 1, j), cap));
            }
        }
    }
    // layer k-1 -> sink
    for i in 0..n_per_layer {
        lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
        let cap = (lcg % max_cap as u64) as i64 + 1;
        edges.push((vid(k - 1, i), t, cap));
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

/// Vision-style 2D grid with t-links at every vertex.
///
/// Mimics image segmentation graphs:
/// - Source `s` and sink `t`
/// - Every pixel has `s→v` (source t-link) and `v→t` (sink t-link) edges
/// - 4-connected bidirectional neighbor edges (n-links)
/// - Capacities: n-links ~ [1, 50], t-links ~ [1, 500]
fn vision_grid_2d(name: &str, rows: usize, cols: usize) -> DimacsGraph {
    let s = 1;
    let t = rows * cols + 2;
    let n = t;
    let vid = |r: usize, c: usize| r * cols + c + 2;
    // t-links: 2 * rows * cols (s→v, v→t for each pixel)
    // n-links: 2 * (rows*(cols-1) + (rows-1)*cols) bidirectional 4-connected
    let n_tlinks = 2 * rows * cols;
    let n_nlinks = 2 * (rows * (cols - 1) + (rows - 1) * cols);
    let mut edges = Vec::with_capacity(n_tlinks + n_nlinks);
    let mut lcg: u64 = 77777;

    // t-links: s→v and v→t for every pixel
    for r in 0..rows {
        for c in 0..cols {
            let v = vid(r, c);
            lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
            let cap_s = (lcg % 500) as i64 + 1;
            edges.push((s, v, cap_s));
            lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
            let cap_t = (lcg % 500) as i64 + 1;
            edges.push((v, t, cap_t));
        }
    }
    // n-links: bidirectional 4-connected
    for r in 0..rows {
        for c in 0..cols {
            let u = vid(r, c);
            // right neighbor
            if c + 1 < cols {
                let v = vid(r, c + 1);
                lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                let cap = (lcg % 50) as i64 + 1;
                edges.push((u, v, cap));
                edges.push((v, u, cap));
            }
            // down neighbor
            if r + 1 < rows {
                let v = vid(r + 1, c);
                lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                let cap = (lcg % 50) as i64 + 1;
                edges.push((u, v, cap));
                edges.push((v, u, cap));
            }
        }
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

/// Vision-style 3D grid (voxel segmentation) with t-links at every vertex.
///
/// Mimics 3D medical image segmentation:
/// - Source `s` and sink `t`
/// - Every voxel has `s→v` and `v→t` edges (t-links)
/// - 6-connected bidirectional neighbor edges (n-links)
/// - Capacities: n-links ~ [1, 50], t-links ~ [1, 500]
fn vision_grid_3d(name: &str, lx: usize, ly: usize, lz: usize) -> DimacsGraph {
    let total_voxels = lx * ly * lz;
    let s = 1;
    let t = total_voxels + 2;
    let n = t;
    let vid = |x: usize, y: usize, z: usize| x * ly * lz + y * lz + z + 2;
    // t-links: 2 * total_voxels
    // n-links: 2 * (3 forward directions) per voxel (approx)
    let n_tlinks = 2 * total_voxels;
    let n_nlinks_est = 2 * (lx * ly * (lz - 1) + lx * (ly - 1) * lz + (lx - 1) * ly * lz);
    let mut edges = Vec::with_capacity(n_tlinks + n_nlinks_est);
    let mut lcg: u64 = 88888;

    // t-links
    for x in 0..lx {
        for y in 0..ly {
            for z in 0..lz {
                let v = vid(x, y, z);
                lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                let cap_s = (lcg % 500) as i64 + 1;
                edges.push((s, v, cap_s));
                lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                let cap_t = (lcg % 500) as i64 + 1;
                edges.push((v, t, cap_t));
            }
        }
    }
    // n-links: bidirectional 6-connected
    for x in 0..lx {
        for y in 0..ly {
            for z in 0..lz {
                let u = vid(x, y, z);
                // +x
                if x + 1 < lx {
                    let v = vid(x + 1, y, z);
                    lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                    let cap = (lcg % 50) as i64 + 1;
                    edges.push((u, v, cap));
                    edges.push((v, u, cap));
                }
                // +y
                if y + 1 < ly {
                    let v = vid(x, y + 1, z);
                    lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                    let cap = (lcg % 50) as i64 + 1;
                    edges.push((u, v, cap));
                    edges.push((v, u, cap));
                }
                // +z
                if z + 1 < lz {
                    let v = vid(x, y, z + 1);
                    lcg = lcg.wrapping_mul(6364136223846793005).wrapping_add(1);
                    let cap = (lcg % 50) as i64 + 1;
                    edges.push((u, v, cap));
                    edges.push((v, u, cap));
                }
            }
        }
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

/// RFIM (Random Field Ising Model) 2D lattice with periodic boundary conditions.
///
/// - d=2 lattice with PBC (toroidal): every vertex has exactly 4 neighbors
/// - Every vertex has either s→v or v→t edge based on random field h_i
/// - Neighbor edges bidirectional with coupling J
/// - Gaussian disorder approximated via Box-Muller from LCG
///
/// Parameters: L (lattice side length), J (coupling constant), sigma (disorder width)
fn rfim_2d(name: &str, l: usize, j_coupling: i64, sigma: f64) -> DimacsGraph {
    let total = l * l;
    let s = 1;
    let t = total + 2;
    let n = t;
    let vid = |x: usize, y: usize| x * l + y + 2;
    // t-links: total (one per vertex, either s→v or v→t)
    // n-links: 2 * 2 * total (4-connected PBC, bidirectional, but each pair counted once = 2*total forward, doubled = 4*total)
    let mut edges = Vec::with_capacity(total + 4 * total);
    let mut lcg: u64 = 31415;

    // Helper: generate approximate Gaussian from LCG using Box-Muller
    let next_gaussian = |lcg_state: &mut u64| -> f64 {
        *lcg_state = lcg_state.wrapping_mul(6364136223846793005).wrapping_add(1);
        let u1 = (*lcg_state as f64) / (u64::MAX as f64);
        let u1 = u1.max(1e-10); // avoid log(0)
        *lcg_state = lcg_state.wrapping_mul(6364136223846793005).wrapping_add(1);
        let u2 = (*lcg_state as f64) / (u64::MAX as f64);
        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
    };

    // t-links: field-dependent
    for x in 0..l {
        for y in 0..l {
            let v = vid(x, y);
            let h = next_gaussian(&mut lcg) * sigma;
            if h > 0.0 {
                let cap = (h * 1000.0) as i64 + 1; // scale to integer
                edges.push((s, v, cap));
            } else {
                let cap = (-h * 1000.0) as i64 + 1;
                edges.push((v, t, cap));
            }
        }
    }
    // n-links: bidirectional 4-connected with PBC
    // For each forward direction (+x, +y), add both (u→v) and (v→u)
    for x in 0..l {
        for y in 0..l {
            let u = vid(x, y);
            // +x (with PBC wrap)
            let v = vid((x + 1) % l, y);
            edges.push((u, v, j_coupling));
            edges.push((v, u, j_coupling));
            // +y (with PBC wrap)
            let v = vid(x, (y + 1) % l);
            edges.push((u, v, j_coupling));
            edges.push((v, u, j_coupling));
        }
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

/// RFIM (Random Field Ising Model) 3D lattice with periodic boundary conditions.
///
/// - d=3 lattice with PBC (toroidal): every vertex has exactly 6 neighbors
/// - Every vertex has either s→v or v→t edge based on random field h_i
/// - Neighbor edges bidirectional with coupling J
/// - Gaussian disorder approximated via Box-Muller from LCG
fn rfim_3d(name: &str, l: usize, j_coupling: i64, sigma: f64) -> DimacsGraph {
    let total = l * l * l;
    let s = 1;
    let t = total + 2;
    let n = t;
    let vid = |x: usize, y: usize, z: usize| x * l * l + y * l + z + 2;
    // t-links: total
    // n-links: 6*total (3 directions * 2 bidirectional, PBC means every vertex has full degree)
    let mut edges = Vec::with_capacity(total + 6 * total);
    let mut lcg: u64 = 27182;

    let next_gaussian = |lcg_state: &mut u64| -> f64 {
        *lcg_state = lcg_state.wrapping_mul(6364136223846793005).wrapping_add(1);
        let u1 = (*lcg_state as f64) / (u64::MAX as f64);
        let u1 = u1.max(1e-10);
        *lcg_state = lcg_state.wrapping_mul(6364136223846793005).wrapping_add(1);
        let u2 = (*lcg_state as f64) / (u64::MAX as f64);
        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
    };

    // t-links
    for x in 0..l {
        for y in 0..l {
            for z in 0..l {
                let v = vid(x, y, z);
                let h = next_gaussian(&mut lcg) * sigma;
                if h > 0.0 {
                    let cap = (h * 1000.0) as i64 + 1;
                    edges.push((s, v, cap));
                } else {
                    let cap = (-h * 1000.0) as i64 + 1;
                    edges.push((v, t, cap));
                }
            }
        }
    }
    // n-links: bidirectional 6-connected with PBC
    // For each forward direction (+x, +y, +z), add both (u→v) and (v→u)
    for x in 0..l {
        for y in 0..l {
            for z in 0..l {
                let u = vid(x, y, z);
                // +x
                let v = vid((x + 1) % l, y, z);
                edges.push((u, v, j_coupling));
                edges.push((v, u, j_coupling));
                // +y
                let v = vid(x, (y + 1) % l, z);
                edges.push((u, v, j_coupling));
                edges.push((v, u, j_coupling));
                // +z
                let v = vid(x, y, (z + 1) % l);
                edges.push((u, v, j_coupling));
                edges.push((v, u, j_coupling));
            }
        }
    }
    DimacsGraph {
        name: name.to_string(),
        n,
        source: s,
        sink: t,
        edges,
    }
}

fn main() {
    let dir = Path::new("../graph-pool");
    fs::create_dir_all(dir).unwrap();

    // We build the graph list in stages so we can track total size.
    let mut all_graphs: Vec<DimacsGraph> = Vec::new();

    // ---- Small (existing, for correctness validation) ----
    all_graphs.push(layered("layered_5x5", 5, 5));
    all_graphs.push(layered("layered_10x10", 10, 10));
    all_graphs.push(layered("layered_20x10", 20, 10));
    all_graphs.push(layered("layered_50x50", 50, 50));
    all_graphs.push(layered("layered_100x100", 100, 100));

    all_graphs.push(grid("grid_5x5", 5, 5));
    all_graphs.push(grid("grid_10x10", 10, 10));
    all_graphs.push(grid("grid_20x20", 20, 20));
    all_graphs.push(grid("grid_50x50", 50, 50));
    all_graphs.push(grid("grid_100x100", 100, 100));

    all_graphs.push(bipartite("bipartite_5x5", 5, 5));
    all_graphs.push(bipartite("bipartite_10x10", 10, 10));
    all_graphs.push(bipartite("bipartite_20x20", 20, 20));
    all_graphs.push(bipartite("bipartite_50x50", 50, 50));
    all_graphs.push(bipartite("bipartite_100x100", 100, 100));

    all_graphs.push(random_er("random_20_10pct", 20, 10, 100));
    all_graphs.push(random_er("random_20_30pct", 20, 30, 100));
    all_graphs.push(random_er("random_50_10pct", 50, 10, 100));
    all_graphs.push(random_er("random_50_30pct", 50, 30, 100));
    all_graphs.push(random_er("random_200_10pct", 200, 10, 100));

    all_graphs.push(chains("chains_5x10", 5, 10));
    all_graphs.push(chains("chains_10x20", 10, 20));
    all_graphs.push(chains("chains_20x10", 20, 10));
    all_graphs.push(chains("chains_50x50", 50, 50));

    // CLRS textbook example
    all_graphs.push(DimacsGraph {
        name: "clrs".to_string(),
        n: 6,
        source: 1,
        sink: 6,
        edges: vec![
            (1, 2, 16),
            (1, 3, 13),
            (2, 3, 10),
            (2, 4, 12),
            (3, 2, 4),
            (3, 5, 14),
            (4, 3, 9),
            (4, 6, 20),
            (5, 4, 7),
            (5, 6, 4),
        ],
    });

    // ---- Medium ----
    all_graphs.push(layered("layered_200x50", 200, 50));
    all_graphs.push(layered("layered_500x100", 500, 100));
    all_graphs.push(grid("grid_200x200", 200, 200));
    all_graphs.push(grid("grid_500x500", 500, 500));
    all_graphs.push(bipartite("bipartite_200x200", 200, 200));
    all_graphs.push(bipartite("bipartite_500x500", 500, 500));
    all_graphs.push(random_er("random_500_5pct", 500, 5, 1000));
    all_graphs.push(random_er("random_1000_5pct", 1000, 5, 1000));
    all_graphs.push(random_er("random_1000_10pct", 1000, 10, 1000));
    all_graphs.push(chains("chains_100x100", 100, 100));
    all_graphs.push(chains("chains_500x100", 500, 100));
    all_graphs.push(washington("washington_10x100", 10, 100, 1000));
    all_graphs.push(washington("washington_5x200", 5, 200, 1000));

    // ---- Large (bulk of the ~1 GB) ----
    all_graphs.push(layered("layered_500x500", 500, 500)); // ~750k edges, ~15 MB
    all_graphs.push(layered("layered_1000x100", 1000, 100)); // ~300k edges, ~7 MB
    all_graphs.push(layered("layered_1000x1000", 1000, 1000)); // ~3M edges, ~70 MB
    all_graphs.push(grid("grid_1000x1000", 1000, 1000)); // ~2M edges, ~50 MB
    all_graphs.push(grid("grid_2000x2000", 2000, 2000)); // ~8M edges, ~200 MB
    all_graphs.push(bipartite("bipartite_1000x1000", 1000, 1000)); // ~1M edges, ~25 MB
    all_graphs.push(bipartite("bipartite_2000x2000", 2000, 2000)); // ~4M edges, ~100 MB
    all_graphs.push(random_er("random_2000_5pct", 2000, 5, 1000)); // ~200k edges, ~5 MB
    all_graphs.push(random_er("random_2000_10pct", 2000, 10, 1000)); // ~400k edges, ~10 MB
    all_graphs.push(random_er("random_5000_3pct", 5000, 3, 1000)); // ~750k edges, ~20 MB
    all_graphs.push(random_er("random_5000_5pct", 5000, 5, 1000)); // ~1.25M edges, ~30 MB
    all_graphs.push(random_er("random_5000_10pct", 5000, 10, 1000)); // ~2.5M edges, ~65 MB
    all_graphs.push(random_er("random_10000_3pct", 10000, 3, 1000)); // ~3M edges, ~80 MB
    all_graphs.push(random_er("random_10000_5pct", 10000, 5, 1000)); // ~5M edges, ~130 MB
    all_graphs.push(chains("chains_1000x1000", 1000, 1000)); // ~1M edges, ~25 MB
    all_graphs.push(washington("washington_10x500", 10, 500, 1000)); // ~2.25M edges, ~55 MB
    all_graphs.push(washington("washington_20x200", 20, 200, 1000)); // ~760k edges, ~18 MB

    // ---- Extra-large: scaling series for structured families ----
    // Layered scaling: 3M -> 6M -> 12M -> 27M -> 30M edges
    all_graphs.push(layered("layered_2000x1000", 2000, 1000)); // ~6M edges
    all_graphs.push(layered("layered_2000x2000", 2000, 2000)); // ~12M edges
    all_graphs.push(layered("layered_5000x1000", 5000, 1000)); // ~15M edges
    all_graphs.push(layered("layered_3000x3000", 3000, 3000)); // ~27M edges
    all_graphs.push(layered("layered_10000x1000", 10000, 1000)); // ~30M edges
    all_graphs.push(layered("layered_10000x2000", 10000, 2000)); // ~60M edges

    // Grid scaling: 2M -> 8M -> 18M edges
    all_graphs.push(grid("grid_3000x3000", 3000, 3000)); // ~18M edges

    // Chain scaling: 1M -> 4M -> 5M edges
    all_graphs.push(chains("chains_2000x2000", 2000, 2000)); // ~4M edges
    all_graphs.push(chains("chains_5000x1000", 5000, 1000)); // ~5M edges

    // ---- Vision-style grids (2D with t-links) ----
    all_graphs.push(vision_grid_2d("vision2d_100x100", 100, 100)); // ~60K edges
    all_graphs.push(vision_grid_2d("vision2d_500x500", 500, 500)); // ~1.5M edges
    all_graphs.push(vision_grid_2d("vision2d_1000x1000", 1000, 1000)); // ~6M edges
    all_graphs.push(vision_grid_2d("vision2d_2000x2000", 2000, 2000)); // ~24M edges
    all_graphs.push(vision_grid_2d("vision2d_5000x5000", 5000, 5000)); // ~150M edges

    // ---- Vision-style grids (3D with t-links, 6-connected) ----
    all_graphs.push(vision_grid_3d("vision3d_50x50x50", 50, 50, 50)); // ~1M edges
    all_graphs.push(vision_grid_3d("vision3d_100x100x100", 100, 100, 100)); // ~8M edges
    all_graphs.push(vision_grid_3d("vision3d_200x200x200", 200, 200, 200)); // ~64M edges
    all_graphs.push(vision_grid_3d("vision3d_300x300x300", 300, 300, 300)); // ~216M edges

    // ---- RFIM 2D (PBC, Gaussian disorder, J=1000, sigma=1.0) ----
    all_graphs.push(rfim_2d("rfim2d_100", 100, 1000, 1.0)); // ~60K edges
    all_graphs.push(rfim_2d("rfim2d_500", 500, 1000, 1.0)); // ~1.5M edges
    all_graphs.push(rfim_2d("rfim2d_1000", 1000, 1000, 1.0)); // ~6M edges
    all_graphs.push(rfim_2d("rfim2d_2000", 2000, 1000, 1.0)); // ~24M edges

    // ---- RFIM 3D (PBC, Gaussian disorder, J=1000, sigma=2.3 ~ near critical) ----
    all_graphs.push(rfim_3d("rfim3d_32", 32, 1000, 2.3)); // ~230K edges
    all_graphs.push(rfim_3d("rfim3d_64", 64, 1000, 2.3)); // ~1.8M edges
    all_graphs.push(rfim_3d("rfim3d_100", 100, 1000, 2.3)); // ~8M edges
    all_graphs.push(rfim_3d("rfim3d_128", 128, 1000, 2.3)); // ~15M edges
    all_graphs.push(rfim_3d("rfim3d_200", 200, 1000, 2.3)); // ~56M edges
    all_graphs.push(rfim_3d("rfim3d_256", 256, 1000, 2.3)); // ~118M edges

    // Write all and report
    let mut total_bytes: usize = 0;
    for g in &all_graphs {
        let est = g.estimated_bytes();
        total_bytes += est;
        print!(
            "writing {}.max  (V={}, E={}, ~{} MB)...",
            g.name,
            g.n,
            g.edges.len(),
            est / (1024 * 1024)
        );
        std::io::stdout().flush().unwrap();
        g.write_to(dir);
        println!(" done");
    }

    println!(
        "\nTotal: {} graphs, estimated ~{:.1} GB",
        all_graphs.len(),
        total_bytes as f64 / (1024.0 * 1024.0 * 1024.0)
    );
}