[speed] Optimise GM99 multi-resident establishment (Poisson sum is exponential in species count)

[dfalster]

Problem

The GM99 toy model (harness_gm99(), Geritz et al. 1999 seed-size safe-site
model) is the slowest of the toy harnesses. Single- and two-resident cases are
fine, but cost explodes with the number of coexisting residents, which makes
community assembly to ≥3 species painfully slow:

residents	one gm99_equilibrium() call
1	<1 ms (closed-form bisection)
2	~3.5 ms
3	~1.1 s

Assembly (assembler_run()) to a 3-species seed-size polymorphism takes ~6 s
per step at 3 residents, and is impractical beyond that — whereas DD99/GK98
assemble dozens of species in well under a second.

Root cause

The bottleneck is the establishment expectation in src/GM99.cpp:

// E[ c_mut / (c_mut + sum_j k_j cres_j) ],  k_j ~ Poisson(N_j) independent
static double estab(double c_mut, const std::vector<double>& cres,
                    const std::vector<double>& N)

It is computed with a full Cartesian odometer over all competitor-count
combinations (k_1, …, k_n), each dimension truncated at
Kmax_j ≈ N_j + 10·sqrt(N_j) + 20. Cost per call is therefore
∏_j (Kmax_j + 1) — exponential in the number of residents (≈ 30³ ≈ 27,000
terms for 3 residents, ≈ 30⁴ for 4).

estab is then called:

• once per resident per iteration of the fixed-point population recursion
  N_i <- N_i · W_i inside gm99_equilibrium() (multi-resident branch,
  max_iter = 10000), and
• once per grid point when building an invasion-fitness landscape against a
  multi-resident community.

So the exponential per-call cost is multiplied by many calls.

Suggested directions (pick what works)

1. Replace the Cartesian odometer with a 1-D convolution. The competition
   sum S = Σ_j k_j c_j is a sum of independent (scaled) Poisson variables.
   Discretise S onto a grid and build its distribution by convolving the
   per-resident contributions (resident j contributes Poisson(N_j) mass at the
   points k·c_j). Then E[g] = Σ_s P(S≈s) · c_mut/(c_mut + s). This is
   O(n · G · Kmax) (linear in residents) instead of ∏ Kmax. Choose the grid
   resolution G to hit a target accuracy; this should preserve the smoothness
   the equilibrium/gradient code relies on.
2. Tighter truncation. Kmax = N + 10·sqrt(N) + 20 is very generous; the
   Poisson tail beyond ~N + 6·sqrt(N) + 6 is < 1e-10. Smaller Kmax helps
   even the odometer.
3. Faster equilibrium solve. The multi-resident equilibrium currently uses
   fixed-point iteration (N_i <- N_i·W_i, up to 10000 steps). A damped Newton /
   nleqslv solve on log N (regnans already depends on nleqslv) would cut the
   number of estab evaluations sharply; add a relative-change stopping rule.
4. Avoid recomputation. Within one equilibrium solve the resident c_j are
   fixed; precompute per-resident Poisson pmfs once per iteration (already partly
   done) and reuse across the n focal evaluations.
5. (Optional) Monte-Carlo estab as a fallback for large n, but note the noise
   would hurt root-finding and finite-difference gradients — prefer the
   deterministic convolution.

Where

• src/GM99.cpp — estab(), gm99_fitness(), gm99_equilibrium().
• R/harness.R — harness_gm99() (no change expected; pars are R, alpha, beta).

Acceptance criteria

• 3-resident gm99_equilibrium() is ≥10× faster (target < 100 ms);
  4–5-resident assembly steps complete in seconds, not minutes.
• Results unchanged within tolerance — the existing oracles must still hold:
  • tests/testthat/test-harness-gm99.R: W_m(m)=1 (log fitness 0),
    g(m,m,N)=(1-e^{-N})/N, the αR=4.5 CSS vs αR=7 branching at βR=15, and the
    αR/βR-only dependence.
  • tests/testthat/test-assembly.R continues to pass.
• A small before/after benchmark (e.g. 2- and 3-resident equilibrium timings)
  reported in the PR.
• (Stretch) GM99 assembles a seed-size polymorphism via assembler_run() to
  ~4 species in comparable time to DD99/GK98.

Context

Added in #33 (the toy-model harnesses). The model is correct and matches the
paper (and Daniel's MATLAB); this is purely a performance issue with the
multi-resident Poisson establishment. See overstorey_staging/GM99.qmd and
overstorey_staging/assembly.qmd.