SwarmSim

Emergent Flocking & Collective Evasion

PHY 351 — Independent Summer Research  |  UNCG

A computational study of force-based flocking agents on a periodic 2D domain, extended with predator-prey dynamics. Based on Charbonneau (2017), Ch. 10.

Where the Law Runs Out — the summer's results site

Three complex systems, one measured claim: typical behavior obeys sharp laws, and the events that matter escape them.

▶ Open Interactive Demo  ·  📋 Research Log

Project status — where we are

Summer complete (2026-07-18): 123 findings across three chapters — flocking F1–F94 (findings.md), sandpiles S1–S22 (findings_sandpile.md), earthquakes E1–E7 (findings_earthquake.md). The argument the summer adds up to, with its central figures, is on the results site: Where the Law Runs Out.

Flocking threads (first chapter), and whether they're closed:

Most recent work — the co-adaptation thread (F87–F94): prey evolve a heritable escape weight under capture/removal predation. The F70 "dangerous valley" is a strong evolutionary brake (F87): escape is stable once present but can't evolve from scratch — it's a barrier to origination, not invasion (F88), and the resulting escape equilibrium is a robust mixed strategy (F89). With a co-evolving predator the arms race is asymmetric (F90), a self-test corrected an over-reading of it (F91), and under an energy-budget fitness model the brake proves model-independent while escape's persistence does not (F92). A second heritable trait, the alignment strength α, is also bistable — high alignment is a protective ESS that can't evolve from low (F93) — but unlike escape it does not invade from a seeded minority either (F94), because alignment is a mutual coupling rather than a free-rideable shared signal. All of it runs on a validated, fast harness (vectorized_predator.py, vectorized_predator_prey.py).

Where to pick up next: see the end of findings.md (Open Questions / Next Directions). The natural next experiments — a third fitness model (proximity-survival) or a two-trait predator (lead + aggression) — each rest on the same harness.

Second chapter — Sandpiles & Self-Organized Criticality

A new topic (Charbonneau Ch. 5), distinct from the flocking work below. Code in sandpile/, findings in findings_sandpile.md (S-series). The arc: validated 1-D core → rigorous critical exponents → the chapter's 2-D "Grand Challenge" → a universality comparison against the canonical abelian (Bak–Tang–Wiesenfeld) sandpile.

Headline result: the phenomenon of self-organized criticality (scale-free avalanches as a dynamical attractor, no parameter tuning) is robust, but the critical exponents are not — they depend on both dimension and toppling rule:

The 2-D slope-vs-BTW split is now confirmed by a second, convention-free measure: the avalanche duration exponent, τ_T ≈ 0.56 (slope) vs ≈ 1.22 (BTW), stable across lattices up to 512 a side (S10). That measurement needed an active-list, numba-compiled engine (Exercise 5, S9) that reproduces the original engines bit-for-bit and runs ~600× faster in 2-D.

And the criticality is exact only because the bulk is conservative: breaking conservation (destroying a fraction of toppled sand) collapses the cutoff scaling from N² to sub-N¹ and truncates the avalanche distribution at a dissipation-set size — singling out conservation among the four SOC ingredients (the OFC earthquake-model sensitivity, in the simplest sandpile).

Validated against the chapter: mass conservation to 3×10⁻¹¹, N-independent power-law avalanche PDFs, initial-condition independence (the SOC attractor), and a clarification of the angle-of-repose claim (the mean slope sits ~16% below Z_c independent of grain size; only the peak pre-avalanche slope approaches Z_c as grains shrink).

sandpile/
  sandpile1d.py    1-D slope model (eqs 5.1–5.10), avalanche measurement
  validate1d.py    conservation, SOC signatures, IC-independence, eps-sweep
  fss1d.py         1-D finite-size scaling: τ_E, D_E, D_T + data collapse
  repose_peak.py   angle-of-repose mean-vs-peak slope reconciliation
  sandpile2d.py    2-D bond-slope generalization (Grand Challenge)
  fss2d.py         2-D finite-size scaling + 1-D-vs-2-D comparison
  btw_compare.py   canonical 2-D abelian BTW, same pipeline (universality)
  dissipation.py   break bulk conservation -> is conservation necessary for SOC?
  dissipation2d.py 2-D version of the conservation test (transfers from 1-D)
  duration_compare.py  duration-exponent universality cross-check (self-test, S6)
  falloff.py       boundary-evacuation avalanches vs bulk avalanches (Exercise 2)
  sandpile_fast.py active-list, numba-compiled 1-D+2-D engine (Exercise 5, S9)
  duration_fss2d.py large-L duration FSS; resolves S6 (S10)
  moments.py       moment-scaling machinery + FSS self-test (S11)
  moment_slope.py  slope-model moment spectrum: filamentary area (S12)
  conditional.py   conditional avalanche exponents; ballistic front (S13)
  geometry2d.py    per-avalanche footprint geometry: the filament (S14)
  stochastic_split.py  tunable split toward Manna; the causal test (S15)
  geometry1d.py    1-D footprint geometry = the dimensional anchor (S16)
  equilibrate2d.py reusable warmup to a verified-stationary L>256 (S16)
  duration_closure.py  single-scale exact in space, residual in time (S17)
  area_multifractality.py  area multifractality is asymptotic, not a corona (S18)
  filament_fattening.py  typical-vs-tail footprint split: the filament does not fatten (S19)
  repose_scaling.py     repose model selection — log(L) beats 1/L by AIC (S20)
  repose_scaling_L1024.py  out-of-sample confirmation at held-out L = 768/1024 (S21)
  manna.py              literal 2-D Manna under the same pipeline — the measured anchor (S22)

Exercises addressed: 2 (falloff avalanches, S8), 3 (initial-condition independence, S1), 4 (eps-dependence of the angle of repose, S2), 5 (fast active-list engine, S9), 6 (the 2-D Grand Challenge, S4) — plus the conservation study (S5/S7) and a self-test that was inconclusive (S6) then resolved (S10).

The scaling-theory arc (S11–S18). Beyond the exercises, the chapter builds one extended argument about what a 2-D slope avalanche actually is. Moment-spectrum analysis (S11–S12) finds the avalanche footprint is filamentary — its area grows only linearly with system size (D_area ≈ 1), against the compact D ≈ 2 of canonical BTW through the same pipeline — while the total toppling activity is space-filling (≈ L²): a thin front that sweeps its own bonds ≈ L times. Conditional exponents (S13) and a direct footprint recording (S14) confirm it geometrically — a constant-width, one-bond-wide, ballistic filament (topple time ≈ distance from the seed), thinner than the exactly-solvable directed sandpile (D = 3/2). A tunable stochastic split toward the Manna class (S15) shows the filament is caused by the deterministic rule (its dimension climbs 1 → 1.87 as randomness is switched on), yet the model never collapses onto simple Manna scaling. A 1-D dimensional anchor (S16) proves D ≈ 1 is intrinsic to a slope avalanche, not a 2-D artifact, and ships an equilibration tool that reaches a verified-stationary L = 512. With it, the last two findings pin down how single-scale the model is: the spatial relations become exactly single-scale as L → ∞ (size ∝ area², area ∝ duration), but duration stays a permanently loose proxy for spatial extent (S17), and the avalanche area distribution stays multifractal asymptotically rather than healing to simple scaling (S18) — the model is single-scale in its means (the typical avalanche) and multifractal in its tails (the largest avalanches). Four closing findings tie the chapter off: the typical filament does not fatten with system size — the apparent fattening at five seeds was tail undersampling (S19); the angle of repose most likely diverges logarithmically with system size rather than saturating, chosen by model selection (S20) and then confirmed out of sample by frozen-fit predictions at held-out L = 768/1024, which land above the saturating model's own ceiling (S21); and the Manna comparison point was measured under the identical pipeline rather than quoted (S22) — closing the placement triangle. It sits in its own place in the SOC landscape: more filamentary than the directed sandpile, outside the stochastic Manna class, a distinct anomalous class from BTW. The whole arc is drawn together in a standalone write-up, report_sandpile.md; full per-finding detail (S-series) is in findings_sandpile.md.

Third chapter — Earthquakes (the Olami–Feder–Christensen model)

A new topic (Charbonneau Ch. 8), the direct continuation of the sandpile conservation thread. Code in earthquake/, findings in findings_earthquake.md (E-series). The OFC model is a 2-D lattice of "force" values driven uniformly until a node hits threshold and topples, passing a fraction α to each neighbor (0 ≤ α ≤ 0.25). At α = 0.25 redistribution is conservative; for α < 0.25 the bulk dissipates (1−4α) per topple — so α is a continuous knob on the exact ingredient (bulk conservation) that finding S5 singled out for SOC.

Headline result: OFC reproduces the Gutenberg–Richter law, and conservation tunes criticality. The avalanche-size cutoff grows with system size when (near‑)conservative and saturates to a dissipation‑set scale when not — the S5 result made quantitative on the canonical earthquake model. And unlike the sandpile, the nonconservative model is deterministic and develops synchronized domains that drive quasi‑periodic recurrent avalanching — temporal structure the sandpile entirely lacks.

Validated against the chapter: the fig.-8.7 PDF slopes (all three match), the fig.-8.6 synchronization domains forming from a random start, the fig.-8.4 recurrent avalanching (raw period ~10,600 at α=0.15, matching the book's ~10,960; forcing-corrected ~3,982 vs the book's 4,002; vanishing at the conservative α=0.25). The size–duration relation steepens with conservation (E ~ T^γ, γ = 1.47 → 1.89 as α: 0.15 → 0.25). A self-test (Exercise 4) confirms the quasi‑periodicity is fragile: a ±0.01 per‑topple jitter in α halves the recurrence peak, because synchronization needs the redistribution to preserve equality exactly.

Grand Challenge (earthquake prediction): a forecaster that phase-locks to the recurrence rhythm beats chance only modestly (1.5–1.7×) and misses the largest events (top‑decile recall 0.09) — reproducing, in a fully-known noise-free model, the book's point that large events sit in an unpredictable tail.

earthquake/
  ofc.py                core OFC model + 6 validation self-tests (open/periodic BC,
                        per-topple stochastic α, Exercise-3 skip-forcing, bit-identical)
  ofc_gr.py             E1: Gutenberg–Richter validation (fig. 8.7)
  ofc_fss.py            E2: finite-size scaling / conservation test (cf. S5/S7)
  ofc_quasiperiodic.py  E3: recurrence period + synchronization domains (figs. 8.4/8.6)
  ofc_et.py             E4: avalanche size–duration relation (Exercise 5)
  ofc_speedup.py        E5: forcing skip-ahead speedup (Exercise 3)
  ofc_stochastic_alpha.py  E6: self-test — does stochastic α break periodicity? (Ex. 4)
  ofc_predict.py        E7: Grand Challenge — earthquake prediction (Exercise 6)

Exercises addressed: 2 (finite-size scaling, E2), 3 (skip-forcing, E5), 4 (stochastic-α self-test, E6), 5 (size–duration relation, E4), 6 (the Grand Challenge prediction, E7) — plus the validation (E1) and the recurrence/synchronization study (E3).

Model

N agents move on a periodic unit square under four forces each timestep:

Key analytical result: equilibrium cruise speed is v_eq = v₀ + α/μ — not v₀. Derived from force balance in an aligned flock; verified experimentally to within 0.002 across four α values.

Default parameters: N=350, r₀=0.005, ε=0.1, r_f=0.1, α=1.0, v₀=1.0, μ=10.0, η=0.5, dt=0.01

Code Architecture

The codebase has two layers.

Core layer (model.py) — object-oriented foundation for all current and future experiments.

Flock          — holds prey positions, velocities, and all physics parameters
                 flock.evolve(predators=[...]) advances one timestep
Predator       — configurable predator agent; strategy='naive'|'encircle',
                 coord_alpha for predator-predator repulsion, enc_radius for encirclement
simulate()     — runs a full experiment, returns Phi and distance timeseries

New experiments import from model.py. To add a new experiment: import Flock, Predator, and simulate from model.py, set parameters, call simulate().

Legacy layer (flocking.py) — procedural implementation that predates the OOP refactor. The earlier experiment scripts (predator.py, multi_predator.py, encirclement.py, etc.) import from here. They are retained intact for reproducibility and are not broken by the model.py refactor.

Repository Layout

root/
  model.py          ← current API. New experiments import from here.
  geometry.py       ← shared helpers (Rg, AR) — used by both old and new scripts.

  flocking.py       ← legacy: procedural baseline.        [keep for reproducibility]
  analysis.py       ← legacy: validation sweeps.          [keep for reproducibility]
  predator.py       ← legacy: single-predator API.        [keep for reproducibility]
  multi_predator.py ← legacy: multi-predator helpers.     [keep for reproducibility]
  encirclement.py   ← legacy: encirclement scaffolding.   [keep for reproducibility]

predator/    predator strategy experiments
contagion/   panic, epidemic, and vaccination experiments
phase/       phase-transition and statistical-mechanics experiments
3d/          three-dimensional extension experiments
figures/     output PNG figures
outputs/     captured text/log output from runs

The five legacy scripts at root are intentionally retained: every numbered finding F5 through F16 was generated by them, and rerunning them must produce identical figures. New experiments should not import from them — they have no if __name__ == "__main__" guard and will side-effect on import. Use model.py instead.

Experiment scripts in subfolders automatically add the project root to sys.path so they can import the core library files above.

Investigation Progression

The experiments follow a logical arc, each motivated by the result before it: