Time-varying demography + demes compliance#88
Merged
nspope merged 101 commits intoMay 6, 2026
Conversation
Spec for adding constant/exponential/linear shape primitives for population sizes and migration rates. Defines NHPP sampler with analytical-inversion closed forms, sweep-phase sizeAt accessor, new 'g' and 'em' event types, and CLI -eg/-eG/-em/-eM. Includes TDD-driven phased plan, sweep parity tests vs current discoal, and neutral parity tests vs msprime under demes models.
Implementation plan covering Phases 0-2 of the design spec: build
src/core/shapes.{h,c} with sizeAt, migAt, integratedHazardSize,
integratedHazardMig, drawWaitingTimeSize, drawWaitingTimeMig for
constant/exponential/linear shapes; add detSweepFreqEuler alongside
detSweepFreq; run Q1 verification harness comparing summary statistic
distributions under constant N. 23 TDD-driven tasks, 111 steps.
Phases 3-8 deferred to subsequent plans.
Stubs for sizeAt, migAt, integratedHazard*, drawWaitingTime* implementations to be filled in subsequent tasks. No runtime behavior change yet.
Code-review followup. Keeps the debug build linkable when shapes module is referenced from other TUs in later tasks. Also adds libyaml/demes-c/libcyaml prerequisites and replaces the hardcoded include path list with $(INCLUDE)/$(EXTERN_LIB) so the debug recipe stays consistent with the other targets (the debug build was already broken by the yaml integration that predates this branch).
Empty Unity test target wired into run_tests. Resets popShape and migShape arrays in setUp so each test starts from a clean slate.
Returns anchor_value regardless of t. Tests verify per-population isolation.
N(t) = anchor_value * exp(-rate_param * (t - anchor_time)). Forward growth rate alpha > 0 maps to backward decline in coalescent time. Tests cover anchor, forward, backward, and the alpha=0 degenerate case.
N(t) = anchor_value + rate_param * (t - anchor_time).
Spec section 3.2 and plan tasks 5/9/10/11/12/15 updated so all shape types use the msprime forward-time per-generation rate convention: positive rate_param means forward-time growth, which means the past held a smaller value. Affected formulas: - LIN size: N(s) = N0 - gamma s (was N0 + gamma s) - EXP migration: m(s) = m0 exp(-beta s) (was m0 exp(beta s)) - LIN migration: m(s) = m0 - delta s (was m0 + delta s) Plus integrated-hazard and waiting-time closed forms updated consistently. Linear zero-crossing now triggers when gamma > 0 (forward growth = backward decline) instead of gamma < 0. EXP size formulas were already in this convention; only the LIN shapes and the EXP-migration spec text needed flipping. Code fix in src/core/shapes.c follows in the next commit.
N(t) = anchor_value - rate_param * (t - anchor_time). rate_param is the per-generation forward-time growth rate (msprime convention); positive value means past was smaller. Tests renamed to reflect: 'declines_backward_with_positive_gamma' and 'grows_backward_with_ negative_gamma'. Sign-convention comment above sizeAt rewritten to describe the uniform forward-time convention across all shapes.
Three sites in the design spec and the plan still encoded the old LIN sign: - spec line 493 (test description in Phase 1) - spec line 601 (Phase 6 zero-crossing description) - plan Task 9 commit-message template All three flipped to match the canonical math in spec section 3.2 and the implemented code at 821625c.
Mirrors sizeAt structure but indexes migShape[src][dst]. Same forward-time msprime convention: positive rate_param means past held a smaller migration rate. Tested per shape and for src/dst pair isolation.
shapes.c: extend the comment block above sizeAt to also document migAt. Same forward-time msprime convention applies uniformly to both sizes and migration rates; the comment now uses generic 'value(t)' wording so it covers both functions. plan: Task 6 block updated to match the implemented LIN sign and test names (declines_backward_with_positive_delta / grows_backward_with_negative_delta), and step 2 updated from '4 tests FAIL' to '5 tests FAIL'. Followup to code-review of commit 0c161ed.
H(T) = (k choose 2) * T / N_0. Returns 0 when k < 2.
Closed form (k choose 2)*(exp(alpha*T) - 1) / (N(t0)*alpha) verified to 1e-6 against midpoint-rule quadrature with 16k subintervals. Includes alpha=0 degenerate case and offset t0. N(t0) is now obtained via sizeAt(popID, t0) so the function handles t0 != anchor_time correctly across all shape types.
Replace exp(a*T) - 1 with expm1(a*T) to avoid catastrophic cancellation when a*T is small. No change in behavior for the typical regime; strictly more accurate near a*T -> 0.
H = (k choose 2) * log(N0 / (N0 - gamma*T)) / gamma. Returns +infinity when N(t) crosses zero within T (gamma > 0) to signal forced coalescence to the waiting-time draw. Uses log1p form for numerical stability when gamma*T/N0 is small.
The midpoint-rule reference at N=16384 has ~4e-6 truncation error for this test's parameterization (integrand 1/N(s) ranges from 0.54 to 20 because N at the far end of the interval is 0.05). The closed form is exact within float precision; tightening N to hit 1e-6 would 4x the test runtime for no extra coverage. 1e-5 still catches any algebraic mistake (sign/factor errors produce diffs of order 1). Plan updated to match: comment added in the test code and the log1p numerical-stability form of the implementation made explicit in the Step 3 snippet.
Constant: k*m0*T. Exponential: k*m0*(-expm1(-beta*T))/beta with expm1 form for numerical stability when beta*T is small. Linear: k*(m0*T - delta*T^2/2). Forward-time convention (msprime); quadrature match to 1e-6.
Closed-form inversion of integrated hazard. log1p used for numerical stability near alpha*xi -> 0. Linear shape (msprime forward-time convention: gamma > 0 means past was smaller) clips T at the zero-crossing boundary using nextafter to ensure strict less-than.
Constant: T = xi/(k*m0). Exponential: T = -log1p(-b*xi/(k*m0))/b with unreachable-xi detection. Linear: smaller positive root of the quadratic k*m0*T - k*delta*T^2/2 = xi. Forward-time convention (msprime).
10k samples compared to theoretical CDF 1 - exp(-H(t)) for each shape. Threshold D < 0.05 is loose enough to avoid false positives but tight enough to catch any inversion error.
5 (N0, gamma) configurations x 1000 random xi each. T must always be in (0, T_cross), no NaN, no inf.
Parallel implementation that advances sweep frequency by one Euler step on the deterministic logistic ODE. Unit test confirms agreement with closed-form detSweepFreq to 1e-3 over a full sweep with 100k steps.
The deterministic-sweep closed form detSweepFreq parameterizes the trajectory in backward (coalescent) time: ttau=0 maps to x near 1 (sweep just fixed), ttau=ts maps to x near epsilon (sweep starting). The backward-time ODE is dx/dttau = -alpha*x*(1-x), NOT +alpha*x*(1-x). The original Task 18 used positive sign, which would walk x in the opposite direction from detSweepFreq when plugged into proposeTrajectory in Task 19, defeating the Q1 verification. Sign flipped in alleleTraj.c, header comment corrected, test rewritten to verify backward-direction agreement (start near 1, end near epsilon, matches closed form to 1e-3 relative). Plan Task 18 updated to match.
Default 'closed' preserves existing detSweepFreq behavior exactly. 'euler' switches deterministic-mode sweeps to the detSweepFreqEuler integrator. Used by the Q1 parity harness to compare distributions empirically.
Runs discoal under closed-form and Euler deterministic-sweep modes across alpha x tau grid, computes niceStats summary statistics, runs KS comparison with Bonferroni-corrected threshold. Output PASS or FAIL with per-statistic detail.
The Euler deterministic-sweep integrator (Task 18) advances x by dx = -alpha*x*(1-x)*dt per step, which for alpha=200, dt=2.5e-8, x near 1 gives |dx| ~ 2.5e-12. Float ULP near 1.0 is ~6e-8, so the increment was silently rounded away and x stayed pinned at 1.0 forever, hitting the 5e8 step safety limit. Changed 'float x' to 'double x' in proposeTrajectory and cast to float at buffer-write time. The buffer is still float[1024] for storage compactness; only the in-flight integrator state needs double precision. Closed-form mode is byte-identical to the prior behavior (verified via diff on the smoke test).
Full grid at 10k reps takes ~7.5 hours on a single core; 1000 reps finishes in ~45 minutes and still gives n=1000 KS statistics with plenty of power to reject distributional differences past the Bonferroni-corrected threshold for ~45 comparisons. Override with REPS=10000 if a longer run is desired.
Lifts the rejection of size_function: exponential / linear epochs in convertDemesToEvents. Computes the per-generation forward-time growth rate from (start_size, end_size, duration) and converts to the discoal internal alpha or gamma (= rate_per_gen * 2 * N, where N is the discoal effective population size; matches the Phase 4b parity-validated conversion). Emits 'g' events for exponential and 'l' events for linear epochs at the more-recent boundary (end_time converted to discoal internal time). The neutral and sweep phases already handle these shapes (Phases 4 and 6). Demes graphs with growth-mode epochs can now be simulated by discoal without manual approximation. Updates the corresponding unit test, which previously asserted that exponential epochs were rejected, to assert that the importer now accepts them.
The hack at discoalFunctions.c:222-261 scanned 'M' events to reconstruct the t=0 migration matrix, working around the importer's broken paired-'M' emission (issue #82). With the Phase 7 importer rewrite, migMatConst[src][dst] is written directly at parser time for the t=0 active matrix; the runtime no longer needs to back-derive anything. Deletes the entire fprintf-laced scanning block plus the migration matrix debug printout. initialize() now does only the standard migMat = migMatConst copy and proceeds. Closes the structural part of issue #82. The msprime parity test in Task 7 confirms behavioral closure on the issue's fixture.
Loads config_examples/demes_example.demes.yaml in both discoal (via the -Y flag with the wrapper YAML config) and msprime (via Demography.from_demes), runs matched replicates at the YAML's sample size, theta, rho, and seeds, and compares summary statistics via Bonferroni-corrected KS at alpha = 0.01. Uses Convention B (ploidy=2, popsize=Ne) — the validated convention from Phase 4b — with per-population SampleSets declared ploidy=1. This works around the [10, 5, 5] haploid sample-count divisibility issue while preserving Convention B's coalescent semantics. Result at REPS=500 (FAIL): discoal: mean ss=165.09, pi=38.7469, td=-0.6836 msprime: mean ss=152.13, pi=38.2243, td=-0.4386 pi: D=0.10 p=1.3e-2 (within Bonferroni) hapdiv: D=0.06 p=2.9e-1 (within Bonferroni) nhap: D=0.06 p=2.9e-1 (within Bonferroni) ss: D=0.31 p=3.5e-22 (REJECT) wtheta: D=0.31 p=3.5e-22 (REJECT) td: D=0.32 p=1.8e-22 (REJECT) Mean pi matches to 1.4% — coalescent times of the genealogies agree, so the importer reproduces the demography correctly. The discrepancy is in the realized SFS: discoal places ~22% more mutations on terminal branches (singletons) than msprime under the same trees, inflating S, watterson theta, and the magnitude of Tajima's D. The same direction of mismatch shows up in phase4b smoke runs without migration, suggesting the SFS shape difference is not specific to the issue #82 fixture or its multi-window migration. The Phase 7 importer rewrite, the back-derivation removal, and the regressions all pass on their own terms; this test records a separate residual discoal vs msprime divergence in the SFS shape. Diagnosing the source of that divergence is left as follow-up work.
Computes tskit branch-mode statistics (diversity, segregating_sites, Tajima's D) directly on tree sequences from both discoal (-ts output) and msprime, bypassing the mutation process entirely. Branch-mode stats discriminate between two hypotheses for the residual site-mode divergence reported in test_issue82_fixture.py at REPS=500: 1. Mutation-placement convention difference (e.g. theta-to-mu). 2. Coalescent-simulator divergence in tree shape. If branch-mode stats match across simulators, the divergence is in mutations (case 1). If branch-mode stats differ, the trees themselves differ (case 2). Tajima's D in branch mode is scale-invariant and compared directly. diversity and segregating_sites differ by an overall time-scale factor (~4 * Ne_diploid = 20000 between discoal coalescent units and msprime generations); each replicate's value is divided by the across-replicate mean before KS comparison so distribution SHAPE is tested instead of absolute values. Complements but does not replace the existing site-mode test. Both tests run on the same fixture (config_examples/demes_example.demes.yaml).
The 'n' event handler in the main events dispatch mutates the currentSize[popID] global to the new relative size. The runtime calls initializeShapesFromGlobals() at the start of each replicate to reset popShape from currentSize, but currentSize itself was never restored to its initial values between replicates. So the mutated value bled into the next replicate, where popShape would be initialized to the post-mutation size instead of the original sample-time size. Result: any non-reference population that had an 'n' size change in replicate N would start replicate N+1 with the wrong size. For multi-population fixtures that combine -en with a downstream -ed/-ej (population merge), the daughter pop in subsequent reps coalesced at the post-'n' rate from t=0 instead of switching at the size-change time. Branch-mode parity vs msprime on the issue #82 fixture caught this as a Tajima's D divergence (D=0.42, p=4e-39 at REPS=500). Fix: snapshot currentSize into currentSizeConst[MAXPOPS] right after getParameters() returns (mirroring how migMatConst captures the initial migration-rate state). At the start of each replicate in initialize(), restore currentSize from the snapshot before initializeShapesFromGlobals() reads it. The unit-test path that drives initialize() directly doesn't allocate currentSize, so the restore is guarded with a NULL check. Branch-mode parity on the issue #82 fixture at REPS=500: - pi: D=0.078, p=0.10 (was 0.088, p=0.04) - segregating sites: D=0.042, p=0.77 (was 0.046, p=0.67) - Tajima's D: D=0.056, p=0.41 (was 0.42, p=4e-39) Site-mode parity also passes now at all 6 statistics (was 3 of 6 rejecting at Bonferroni p < 1.7e-3). All bit-equality regressions still pass (the bug was deterministic given a fixed seed; runs with the same seed produced the same wrong output, and the new behavior produces new output that is bit-equal across re-runs).
…e restore Cleanup pass on the currentSize replicate-bleed fix: - Move the restore loop in initialize() down next to its consumer (initializeShapesFromGlobals at line 215) instead of at the top of the function. Locality with the only consumer. - Drop the NULL guard on currentSize in production code; instead allocate currentSize in test_node_operations.c's setUp(), matching test_shapes.c's existing pattern. - Add 'extern double *currentSize;' to discoal.h so the snapshot copy in main() and the restore in initialize() compile without redundant local extern declarations. - Trim the narrative comments to one short line each.
…rsion Phases 4-7 introduced time-varying demographic event flags that were not yet reflected in the user-facing docs: - README.md and binary usage now list -eg/-eG (exponential growth), -el/-eL (linear growth), -em/-eM (migration rate change), and note -ej as an alias for -ed. - yaml_configuration.md no longer claims time-varying migration is unsupported; it points users at -em/-eM and the demes loader. - demes_integration.md no longer lists exponential or linear growth epochs as unsupported (the importer translates them to -eg/-el events). - demes_integration.md time-conversion formula corrected from /4N to /(2*N_ref) to match the actual importer code in demesInterface.c.
The importer's demesTimeToCoalTime() multiplied demes_time by generation_time, which is the wrong direction. For time_units=years with generation_time=25, a 50000-year event was converted to 50000 * 25 / (2N) coal units instead of 50000 / 25 / (2N) — off by generation_time^2. The bug was hidden because all existing fixtures used time_units=generations, where demes-c validates generation_time=1 and multiply-by-1 is a no-op. Fix: divide instead of multiply. Same correction applied to the exp/linear epoch alpha and gamma calculations, which used (start_time - end_time) directly as a generation count. Adds test/parity/phase8_demes_units/ with two demes fixtures describing the same demographic model — one in time_units=years, one in time_units=generations — and a regression script that verifies they produce byte-equal output. Discoal now handles every demes feature except selfing and cloning rates, which are not expressible in the coalescent model.
The plans and design spec under docs/superpowers/ were working notes used while developing the time-varying demography feature. They don't belong in the published codebase. Adding the directory to .gitignore and removing the eight files from the index; local working copies are preserved for reference.
The 'Download Unity test framework' step in the Memory Leak Check job ran 'git clone' unconditionally, which failed when Unity was already present on the runner. The Basic Functionality Tests job already guards the clone with 'if [ ! -d "extern/Unity" ]'; mirror that pattern here.
The phase3-phase8 parity smoke and bit-equality scripts plus the msprime-comparison harnesses under test/parity/ were development infrastructure used to validate the time-varying-demography branch. They reference msprime, demes, scipy and other Python deps that the upstream test suite doesn't require, so they don't belong in the published repo. Local working copies are preserved.
Adds validateShapeTrajectories() in shapes.c, called from
getParameters() right after sortEventArray(). It walks the events
array in time order, tracks each population's shape state, and
exits with an explanatory error on any of:
- An -en event whose target popnSize is <= 0.
- A linear-growth ('l') shape whose forward-time slope drives the
population size to zero before the next event involving that
population (e.g., before a merge or a subsequent -en/-eg/-el).
- A computed size that turns out non-positive at the firing time
of an -en/-eg/-el event (e.g., from a prior linear shape that
carried the size below zero between events).
Exponential shapes never reach zero in finite time, so they don't
trigger the zero-crossing path. Merged populations (after a 'p'
event) stop being checked, since their shape is irrelevant once
all their lineages have moved.
Tests against -en sizes of 0 and -1 reject; -el with a slope that
crosses zero at internal t=1.5 between an -el at t=1.0 and a
merge at t=2.0 rejects with the specific crossing time. Existing
bit-equality regressions and unit tests still pass.
nspope
approved these changes
May 3, 2026
There was a problem hiding this comment.
Looks good ... I get the gist here, but am worried I'm not seeing bugs as there's a lot of additions. I think it'd be worth adding some examples to one of the testing/* suites (or a new suite) to ensure that it at least runs through on a handful of different selection+demography combinations ... and gives different answers when the seed is fixed but e.g. the demography changes.
It'd probably be worth including some failing cases too---like where the linear size/migration changes are known to go negative.
Comments referenced the implementation plan ("Phase 4 wires...", "as of
Phase 7 Task 5", "pre-Phase-4 behavior", "implemented in subsequent
tasks"). These describe project history, not the code, and rot as the
plan recedes. Reworded or dropped them; behavior unchanged.
The default branches in shapes.c switches returned 0.0 (or -1.0 in the draw functions) for unknown shape types. Today these are unreachable, but if a new SHAPE_* enum is added without updating one of these switches, the silent fallback would propagate as wrong rates: zero migration, zero coalescence, or a -1.0 that masquerades as the linear-zero-crossing sentinel. Replaced with a fatal unknownShape() helper that prints the function name and exits, so any future omission fails loudly at first exercise rather than silently corrupting a simulation.
…ssing Each switch case in shapes.c now carries a header block that derives the closed-form integral (or its inverse) directly above the code, so a reader can verify the formula without re-deriving it. The header text also flags the unreachable-xi conditions and the divergent / plateau edges of the linear branches. The linear branch of integratedHazardMig had a math bug: H_raw(T) = k*(m0*T - 0.5*d*T^2) is the antiderivative of the unclamped linear rate, but the physical migration rate is non-negative, so for d > 0 and T past T_zero = m0/d the integrand is identically zero and the integral plateaus at k*m0^2/(2*d). Without the clamp H_raw peaks at T_zero, decreases, and goes negative past T = 2*m0/d -- a meaningless integrated hazard. drawWaitingTimeMig's discriminant guard prevents production from exercising the negative range today, but integratedHazardMig is now self-consistent for any caller. Added a regression test that pins the plateau at T_zero and at T = 100*T_zero.
Nate flagged that the comment listed the CLI parser and the demes importer but missed the YAML config loader (configInterface.c:591). All three paths populate migMatConst before the runtime copies it into migMat at the start of each replicate.
Member
Author
|
@nspope -- this is ready for you now |
|
Looks good ... did you want to add some validation examples here (incl. ones that should fail b/c of negative pop sizes/rates), or do that in another PR? |
Drives the discoal binary with 13 CLI cases that exercise validateShapeTrajectories: 7 well-formed inputs that must exit 0 and produce '//' segment output (covering constant single/two-pop, positive -en, exponential growth, short-interval linear growth, zero -em that turns migration off, and -em involving a merged pop), and 6 ill-formed inputs that must exit non-zero with the documented stderr substring (zero/negative -en sizes, linear -el whose zero crossing precedes the next event, negative -m/-M initial migration, negative -em rate). Mirrors the style of testing/yaml_parser_suite.sh.
Member
Author
|
just added |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #82.
Summary
-eg/-eG,-el/-eL,-em/-eM, plus-ejas alias for-ed.time_units: years.discoalFunctions.c:222-261(which approximated the t=0 migration matrix from paired'M'events) with directmigMatConstwrites from the importer.currentSizebleed between replicates that produced wrong trees for any non-reference pop with an-enevent.Scope of demes compliance
Every demes feature except selfing and cloning rates. Time units, exp/linear epochs, multi-window migration, pulses, and admixture all go through. The time units stuff is still subject to a careful parity check! I think we should run this against msprime for all of
stdpopsimValidation
75 unit tests via
make run_tests(added coverage for the shapes module anddetSweepFreqGeneral).