[evol assembly] Fast toy models with invasion fitness for dev & validation

[dfalster]

Problem statement

Developing and validating the assembly / attractor machinery in regnans
(community_demography, community_selection_gradient,
community_solve_singularity_1D, community_fitness_landscape, the
births_*/deaths/assembler loop) is slow and hard to verify because every
fitness or equilibrium evaluation runs the full plant SCM (run_scm), which
dominates test wall-clock (see test-solve-attractors.R) and has no
closed-form answer to check the algorithms against.

We want one or more fast toy models with analytic (or near-analytic) invasion
fitness to plug into the same community interface, so we can:

• iterate on algorithms in milliseconds rather than SCM-seconds;
• validate numerics against known answers — singular strategies, branching
  points, and viable bounds that are analytically derivable;
• exercise paths the SCM rarely reaches cheaply — especially the nD
  attractor solver (an unmet acceptance criterion in [evol assembly] Enable evolutionary dynamics & community assembly #22) and the
  non-branching / convergence and non-point-attractor cases.

This is the detailed version of the placeholder #21 (which has an empty body) and
feeds the epic #22. Suggest closing #21 as superseded by this.

Integration point

A toy model only needs to satisfy the abstraction the community object already
assumes. community$model_support + the plant_community_* aliases in
R/community_plant.R are explicitly designed to be swappable (AGENTS.md: "This
indirection is so the plant-specific layer could in principle be swapped"). A toy
backend supplies:

• fitness_function(x_new, x, y) → invasion fitness of mutant(s) x_new in
  resident community (x, y);
• a demographic-equilibrium solver for y (analytic where possible, else the
  existing equilibrium_iteration on the difference equation
  y_i(t+1) = y_i(t)(1 + f(x_i, x, y))).

The job is to add a non-plant model_support backend (e.g.
toy_model_support(make_dieckmann_1999())) and confirm the whole community
pipeline runs against it unchanged. The old traitecoevo/Revolve package
(x_misc/Revolve) already implements several of these as make_* closures
returning exactly fitness / equilibrium / single_equilibrium — they can be
ported with light adaptation.

Review of the Revolve catalogue

Source: x_misc/Revolve/doc/models.md + x_misc/Revolve/R/.

Model	In Revolve R/?	Dim	Behaviour	Why useful for dev
Dieckmann & Doebeli 1999 — Gaussian competition + carrying-capacity kernel	✅ dieckmann-1999.R	1D (nD variant exists)	Evolutionary branching; analytic singular strategy at x=0	The canonical AD toy model. Fast, analytic singularity → gold standard for the 1D gradient/singularity solver. Top pick.
Kisdi 1999 — asymmetric (logistic) competition	✅ kisdi-1999.R	1D	Branching + directional; analytic singular strategies (linear & gaussian growth)	Second analytic 1D check; asymmetric kernel exercises a different regime than DD99. Top pick.
Fox & Vasseur 2008 — competition for two essential (Tilman) resources	✅ fox-2008.R	1D trait	Character convergence (not branching); ODE/resource-explicit with analytic resource equilibria	A non-branching counterpoint, and a resource-explicit (R*) structure closer to plant's mechanism. Good.
Geritz et al. seed-size in safe sites (1988/1995/1999)	❌ (described only)	1D	Polymorphism in seed size; patch/Poisson safe-site structure	Most plant-relevant: patch-structured invasion fitness conceptually parallel to plant's patch model, on a real plant trait. Higher implementation cost (Poisson sum over patch occupancy). Recommended.
Ito & Dieckmann 2007 — multi-trait, directional + disruptive	❌ (described only)	2D	Branching in one trait, directional selection in another; recurrent radiations	Directly fills #22's nD attractor gap — the cheapest way to test 2D gradients/singularities. Recommended.
R* (Tilman / Huisman–Weissing)	✅ rstar.R	nD resources	Resource competition; can give oscillations/chaos with ≥3 resources	Useful later as a non-point-attractor stress test; heavier (RefClass, ODE).
Multi-resource DD99 (Ito & Dieckmann 2007)	❌	nD	nD branching	Alternative nD case; parameters never published, so less attractive than Ito 2D directional.

Also present: mutation.R (mvnorm mutation generator — directly reusable by the
stochastic assembler) and utils.R (equilibrium_sys, ODE/iteration helpers).

Recommendation — start with two well-understood "successful" examples

Begin with two models whose outcomes are clean, classic, and analytically
verifiable — one convergence to a stable endpoint, one branching — so the
core algorithms can be validated against known answers before tackling harder
cases:

Starter 1 — Migratory-bird arrival time (Brännström et al. 2013, §4)

A discrete-time, single-trait model (trait x = arrival time) with
frequency dependence via competition for K territories. Fully analytic — every
quantity below can be checked in a unit test. Outcome: convergence to a single
continuously stable strategy (CSS), no branching — the canonical "clean
success" case.

• Competitive ability C(x) = exp(-a·x) (early arrivers win territories).
• Reproductive output R(x) = R₀·exp(-(x - x_opt)² / 2σ²) (peaks at x_opt).
• Year-to-year survival p; resident demography n_{t+1} = K·R(x) + p·n_t,
  so n* = K·R(x)/(1-p).
• Invasion fitness (geometric, dimensionless; w=1 at equilibrium):
  w_x(x') = (1-p)·C(x')R(x') / (C(x)R(x)) + p.
• Selection gradient D(x) = -(1-p)·(a + (x - x_opt)/σ²).
• Singular strategy x* = x_opt - a·σ² — provably convergence-stable and an
  ESS (a CSS).

Why first: trivial to implement, exercises the discrete-time equilibrium path and
community_selection_gradient / community_solve_singularity_1D, and gives an
exact x* to assert against. A nice ecological story too (tragedy of the
commons — birds arrive earlier than the population optimum).

Starter 2 — Geritz seed-size in safe sites (branching)

The Geritz seed-size / asymmetric-competition-for-safe-sites model (Geritz et al.
1988/1995/1999; the canonical AD treatment is Geritz, Kisdi, Meszéna & Metz
1998, Evol. Ecol. Res. 12:35–57). A 1D evolutionary-branching case on a real
plant trait (seed size), with an explicit PIP and a convergence-stable
fitness-minimum singular strategy — the standard demonstration that the
machinery detects branching, complementing the bird model's CSS endpoint. Also
the most plant-relevant toy model (patch-structured, seed size).

An implementation already exists — Daniel's MATLAB version reproducing
Geritz et al. 1999 Fig. 5 PIPs, at
OneDrive-UNSW/research/directions/Offspring-SmithFretwellReview/models/Geritz/.
The task is a MATLAB→R port. The pieces map directly onto the toy-model
interface and match x_misc/Revolve/doc/models.md:

• pre-competitive survival s(m) = max(0, 1 - 2·exp(-β·m)) (f.m);
• competitive ability c(m) = exp(α·m) (c.m);
• establishment g(m', m, N) = Σ_k e^{-N} N^k/k! · c(m')/(c(m') + k·c(m))
  — Poisson sum over seeds arriving per safe site (g.m);
• invasion fitness f̂ = (R/m')·s(m')·g(m', m, N) (Fit.m);
• resident equilibrium density N* by 1D root-find on the self-competition
  balance (res_DE.m + g_res.m).

So invasion fitness is effectively analytic (a truncated Poisson sum) and the
resident solve is a cheap 1D uniroot — no SCM, fast.

Gotcha for assembly: the MATLAB g() handles a single resident only. For
multi-species communities (which regnans assembles) the establishment term needs
the full multi-resident Poisson product-sum from models.md eq. (g). Start with
the 1-resident case (enough for the PIP / singular-strategy / branching tests),
then generalise.

(If the Poisson form proves fiddly to port, the Revolve-implemented DD99 / Kisdi
1999 below are drop-in stand-ins for a branching test.)

Follow-ups (after the two starters work)

3. Dieckmann & Doebeli 1999 + Kisdi 1999 — near-free ports of
   dieckmann-1999.R / kisdi-1999.R; extra analytic branching checks.
4. Ito & Dieckmann 2007 (2D directional) — implement fresh (params in
   models.md); the cheap case to build and validate the nD attractor solver
   ([evol assembly] Enable evolutionary dynamics & community assembly #22) and 2D fitmax births (find_max_fitness_2d, unverified — [evol assembly] Error in find_max_fitness_2d if trait value on the boundary #9).

Other suggestions from the literature (beyond models.md)

• MacArthur–Roughgarden L–V competition kernel (MacArthur 1970; Roughgarden
  1972) — the fully-analytic mechanistic ancestor of DD99; the simplest possible
  limiting-similarity test.
• Rosenzweig–MacArthur predator–prey AD (Dieckmann, Marrow & Law 1995;
  Dercole & Rinaldi 2008) — produces evolutionary limit cycles / Red Queen
  dynamics; a point-attractor solver should fail here, so a valuable
  robustness/negative test (leave for later, not a "successful example").
• Doebeli & Dieckmann 2000 (Nature) / Doebeli 2011 Adaptive
  Diversification — competition along an environmental gradient; many small
  worked examples.
• Brännström, Johansson & von Festenberg 2013, "The Hitchhiker's Guide to
  Adaptive Dynamics", Games 4:304–328 — the source of Starter 1; tutorial
  review with worked models and pseudocode for PIPs, branching, the canonical
  equation, and trait-evolution (dimorphic) plots. The best implementation /
  diagnostic reference to follow throughout.

Acceptance criteria

• A non-plant model_support backend exists; the core community pipeline
  (community_start → community_add → community_demography → community_selection_gradient) runs against a toy model with no SCM call.
• Starter 1 — bird arrival-time model implemented; a test asserts the
  recovered singular strategy matches x* = x_opt - a·σ² (a CSS / clean
  convergence).
• Starter 2 — Geritz seed-size safe-site model ported from MATLAB to R; a
  test confirms a convergence-stable fitness-minimum singular strategy
  (branching detected) and reproduces the Geritz et al. 1999 Fig. 5 PIP.
• (Follow-up) DD99 / Kisdi 1999 ported as extra branching checks.
• (Follow-up) A 2D toy model (Ito & Dieckmann 2007) drives and validates the
  nD attractor path (ties to [evol assembly] Enable evolutionary dynamics & community assembly #22).
• Toy-model tests are fast enough to run in the default suite without the SCM
  wall-clock cost.

Refs: #21 (supersedes), #22 (epic), #9 (find_max_fitness_2d on boundary).
Source material: x_misc/Revolve/doc/models.md, x_misc/Revolve/R/.

[dfalster]

Implemented: model-harness abstraction + three analytic toy models

A swappable harness now wires a community to its demographic backend, so the
assembly/attractor pipeline can run against fast toy models with analytic answers
instead of always running the plant SCM. Landed on master
(449ade8, 025f027, f0ebed7).

Design

The pipeline calls six model-agnostic connectors (R/harness.R), each forwarding
to the harness on the community:
community_parameters, community_make_demography_runner,
community_demography_runner_cleanup, community_viable_bounds,
community_check_for_inviable_strategies, community_update_fitness_function.

• harness_plant() — the default; a thin wrapper over the existing
  plant_community_* code (plant path unchanged; callers passing only
  model_support keep working). version= hook for different plant interfaces.
• harness_analytic() — implements all six connectors generically from two
  primitives: a vectorised log-invasion-fitness function and an equilibrium solve,
  both in C++ (src/toy_models.cpp). This stands up the package's first C++
  layer (Rcpp via LinkingTo, @useDynLib, src/Makevars C++17).

Models (each with analytic test oracles)

harness	x*	branching	source
harness_bird()	x_opt - a·σ² (CSS, no branching)	—	Johansson & Jonzén 2012; Brännström 2013 §4
harness_dd99()	x0	iff σ_C < σ_K	Dieckmann & Doebeli 1999
harness_geritz98()	capacity-weighted mean optimum (0 symmetric)	iff d/σ > √(3/2) ≈ 1.2247	Geritz et al. 1998 soft-selection

Note: harness_geritz98 is the 1998 paper's soft-selection (Levene) worked
example — not the Geritz 1999 seed-size safe-site model (the MATLAB code in
Daniel's OneDrive). The 1999 seed-size model remains a possible future addition.

Validation

tests/testthat/test-harness-{bird,dd99,geritz98}.R assert each singular strategy
is recovered by community_solve_singularity_1D, and that the invasion-fitness
curvature flips sign at the analytic ESS/branching boundary (the Geritz threshold
is recovered to ~1e-3 by uniroot). These run no SCM, so they are fast.
Suite: 177 pass, 0 fail (was 133; plant path unchanged).

Status vs. acceptance criteria

• Non-plant backend; core pipeline runs against a toy model with no SCM call.
• Bird model — analytic x* = x_opt - a·σ² recovered.
• A branching model with analytic oracle — DD99 and Geritz98 (two, both
  with curvature/threshold oracles).
• Fast toy-model tests in the default suite.
• nD attractor path (Ito & Dieckmann 2007 / multi-trait) — not yet; these
  three models are all 1D.
• Multi-resident assembly fully exercised — single-resident oracles done;
  multi-resident equilibria implemented (DD99 linear solve, bird/Geritz
  recursion) but assembly-loop coverage for toy backends is future work.

[dfalster]

Added: Geritz 1999 seed-size safe-site model (harness_geritz99())

The fourth toy model is in (commit 31ee236) — the Geritz, van der Meijden & Metz
1999 seed-size model (Poisson safe sites + size-asymmetric seedling lottery), i.e.
the model in Daniel's MATLAB. Distinct from harness_geritz98 (1998
soft-selection).

C++ in src/toy_models.cpp: s(x)=max(0,1-2e^{-βx}), c(x)=e^{αx}, establishment
averaged over a Poisson number of competitors (stable pmf recurrence; N-dim
odometer for multi-resident), W=(R/x)s(x)g (lifetime reproductive ratio). Only
α·R and β·R matter. No closed-form x*, so the strong oracle is the
equilibrium invariant W_m(m)=1.

Tests reproduce the paper's Fig. 5 behaviour at βR=15: αR=4.5 → CSS
(m*≈0.183, curvature<0), αR=7 → branching (m*≈0.645, curvature>0); plus
g(m,m,N)=(1-e^{-N})/N, non-viable sizes → N*=0, and αR/βR-only dependence.

Four toy models now ship: harness_bird, harness_dd99, harness_geritz98,
harness_geritz99. Suite: 197 pass, 0 fail. Remaining future work unchanged:
an nD model for the nD-attractor path, and fuller multi-resident assembly coverage.

[dfalster]

Consistency pass: naming, per-model C++ files, parameterised plant harness, write-ups

Following review, tidied the API for consistency (commit 62cb440):

• Author/year harness names: harness_dd99, harness_gk98 (was geritz98),
  harness_gm99 (was geritz99), harness_jj12 (was bird).
• harness_analytic → harness_explicit. The family share is that fitness +
  equilibrium are supplied as explicit (C++) functions, not that everything is
  closed-form — gm99's fitness is a Poisson series and its singular strategy is
  numerical. (Honest naming; see the table below.)
• One C++ file per model: src/DD99.cpp, GK98.cpp, GM99.cpp, JJ12.cpp
  (was a single toy_models.cpp); exported as <code>_fitness / <code>_equilibrium.
• Plant harness parameterised like the toy ones: harness_plant(model = "FF16"/"TF24", version = <plant pkg version>).

What's actually analytic vs numerical (the reason for the rename):

	single-resident eq.	fitness	x*	multi-resident eq.
dd99	closed form	closed form	closed form x0	direct linear solve
gk98	closed form	closed form	closed form Σc_jμ_j	iteration
jj12	closed form	closed form	closed form x_opt−aσ²	iteration
gm99	numerical (bisection)	Poisson series	numerical	iteration

Suite still 197 pass, 0 fail.

Per-model write-ups

Added Quarto write-ups (one per model: original-paper summary, fitness equations,
worked regnans example incl. singular-strategy solve + PIP) under
overstorey_staging/ — DD99.qmd, GK98.qmd, GM99.qmd, JJ12.qmd (+ index).
A staging area for the overstorey docs site.

[dfalster]

Community assembly (beyond PIPs) + multi-trait DD99

Commits 0637f5d, 4934f5c. The DD99/GK98/GM99 models handle mixtures, so the
full assembler runs on them — not just pairwise invasibility plots.

Trait-scale option. community_start(trait_scale = "log"/"linear")
(community_trait_transform()) controls how trait space is spaced/searched
during assembly. "log" (default) preserves the plant path for strictly-positive
traits; "linear" lets DD99/GK98 (optima at 0) assemble without shifting them
positive. Threads through the fitness-landscape grid, nearest-resident distance
(births / should_move), and find_max_fitness_2d.

1D assembly demonstrated:

• DD99 (linear) → packed limiting-similarity community (~12 species across [-3,3]).
• GK98 (linear, seeded — soft selection has no empty habitat) → branches to a
  symmetric protected dimorphism (±0.85), matching the paper's Fig. 6.
• GM99 (log, from empty) → seed-size polymorphism (e.g. 0.11 / 0.45 / 0.69).
  Slowest model (multi-resident fitness is a multi-dimensional Poisson sum).

DD99 in 2 trait dimensions — yes. harness_dd99_nd() (cf. Ito & Dieckmann
2007; product-Gaussian kernels), with dd99_nd_fitness/dd99_nd_equilibrium in
src/DD99.cpp and a now dimension-aware explicit harness. Singular point at
x*=x0=(0,0), branching curvature when σ_C<σ_K per dimension. Assembles ~25
coexisting species across the 2D trait plane (limiting similarity in 2D).

Tests: new test-assembly.R (trait scale, DD99 1D packing, GK98 dimorphism, DD99
2D oracles + multi-resident). Suite 214 pass, 0 fail, 0 warn.

Write-up: overstorey_staging/assembly.qmd (DD99 1D/2D, GK98, GM99 demos with plots).

Caveat / ties to #9: the built-in assembler's 2D fitness search
(find_max_fitness_2d) runs but emits warnings from a vestigial 1D
fitness-landscape step; the 2D demo therefore uses a small public-API max-fitness
loop. Hardening the assembler's nD path (drop/skip the 1D landscape for nD, wire
fitness_slopes) is the remaining work, tracked under #9.