Commit 0ad2d5a
authored
Time-varying demography + demes compliance (#88)
* Design: time-varying demographic parameters (issue #82)
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.
* Plan: foundations (shape primitives + Q1 verification)
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.
* Add Shape type and shapes module scaffold
Stubs for sizeAt, migAt, integratedHazard*, drawWaitingTime*
implementations to be filled in subsequent tasks. No runtime
behavior change yet.
* Add shapes.c to discoal_debug Makefile recipe
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).
* Add test_shapes runner scaffold
Empty Unity test target wired into run_tests. Resets popShape
and migShape arrays in setUp so each test starts from a clean
slate.
* Implement sizeAt for SHAPE_CONSTANT
Returns anchor_value regardless of t. Tests verify per-population
isolation.
* Implement sizeAt for SHAPE_EXPONENTIAL
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.
* Implement sizeAt for SHAPE_LINEAR
N(t) = anchor_value + rate_param * (t - anchor_time).
* Adopt msprime forward-time convention for shape rate_param
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.
* Fix sizeAt SHAPE_LINEAR sign to forward-time convention
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.
* Patch stale residuals from forward-time convention sweep
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.
* Implement migAt for all shapes
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.
* Generalize sign-convention comment to cover migAt; patch plan Task 6
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.
* Implement integratedHazardSize for SHAPE_CONSTANT
H(T) = (k choose 2) * T / N_0. Returns 0 when k < 2.
* Implement integratedHazardSize for SHAPE_EXPONENTIAL
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.
* Use expm1 for integratedHazardSize EXP numerical stability
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.
* Implement integratedHazardSize for SHAPE_LINEAR
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.
* Relax integratedHazardSize LIN quadrature tolerance to 1e-5
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.
* Implement integratedHazardMig for all shapes
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.
* Implement drawWaitingTimeSize for all shapes
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.
* Implement drawWaitingTimeMig for all shapes
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).
* Add KS distribution tests for drawWaitingTimeSize
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.
* Add KS distribution tests for drawWaitingTimeMig
* Add linear-shape zero-crossing stress test
5 (N0, gamma) configurations x 1000 random xi each. T must
always be in (0, T_cross), no NaN, no inf.
* Add detSweepFreqEuler for Q1 verification
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.
* Fix detSweepFreqEuler sign for backward-time integration
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.
* Add --det-sweep-mode runtime flag for Q1 verification
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.
* Add Q1 verification harness for detSweepFreq vs Euler
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.
* Use double precision for x in proposeTrajectory
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).
* Make Q1 harness REPS configurable, default 1000
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.
* Q1 verification result: detSweepFreq vs Euler under constant N
FAIL: 29/108 comparisons reject equality at Bonferroni-corrected
p < 9.26e-05 (alpha = 0.01 / 108 comparisons).
Pattern of failures:
- Recent sweeps (tau=0.01) with smaller alpha (50, 200) reject
heavily (11/12 stats reject in each).
- Larger alpha (1000) does not reject at any tau.
- Older sweeps (tau=0.5) do not reject regardless of alpha.
Conclusion: under the tested grid the closed-form detSweepFreq
and the Euler-step detSweepFreqEuler are NOT statistically
indistinguishable. The conservative dispatch on shape type
specified in design spec 4.5 must be kept: constant-N runs
use the closed form; Euler is only used when N is genuinely
non-constant (Phase 5 of the design).
Run: REPS=1000, n=10, theta=10, rho=10, nsites=10000,
alpha in {50, 200, 1000}, tau in {0.01, 0.1, 0.5}.
Raw .ms and .stats outputs are gitignored (large and
reproducible by re-running the harness).
Full per-comparison detail: test/parity/q1_results/analysis.txt
* Spec addendum: Q1 verification result (FAIL)
* Spec revision: deterministic sweep uses closed-form, not Euler
The ODE dx/dtau = -alpha_eff(tau)*x*(1-x) is separable, so under
arbitrary time-varying N it has a closed-form solution
x(tau) = x_0 * exp(-A(tau)) / (1 - x_0 + x_0 * exp(-A(tau)))
where A(tau) = alpha_0 * integral_0^tau sizeAt(s) ds. For our shape
vocabulary the integral is closed-form per shape (constant, exp,
linear). Phase 5 implements detSweepFreqGeneral plus a sister
integratedSizeRatio accessor; the Euler step and --det-sweep-mode
flag are retired in Phase 5.
Sections 3.3 and 4.5 of the design rewritten to reflect this.
The Q1 addendum gets a "revised conclusion" paragraph noting that
the FAIL outcome motivates the closed-form pivot.
The Q1 harness, --det-sweep-mode flag, and detSweepFreqEuler stay
on the branch through Phase 5 as a record; they will be removed
when the new closed-form path lands.
* Plan: Phase 3 shape accessor wiring (bit-equal foundation)
Wires sizeAt and migAt into the simulation engine with SHAPE_CONSTANT
only. 8 tasks, 57 steps. Bit-equality vs pre-Phase-3 binary across
8 configurations (single-pop, two-pop with size changes / migration /
split, three-pop combos, admixture, ancient samples) is the
verification criterion.
The NHPP per-component draw is deferred to Phase 4 (when EXPONENTIAL
is added and the rate-varies-within-interval problem actually
requires it). Sweep code, importer, and the back-derivation hack
are unchanged in Phase 3.
Reference binary discoal_pre_phase3 is built from a temp git
worktree at the Phase 3 base SHA (currently 18ddde1).
* Add discoal_pre_phase3 Makefile recipe for Phase 3 regression
Builds discoal from the Phase 3 base SHA in a temp git worktree
and copies the binary to build/discoal_pre_phase3. Used as the
bit-equality regression target during Phase 3 wiring tasks.
* Fix discoal_pre_phase3 idempotence with file-target pattern
The original recipe placed the early-exit `if` on its own Make
line, which spawns a separate shell from the build line. `exit 0`
exited only the if-shell and the build ran anyway. Replace with
a phony alias depending on a real-file target -- make's standard
file-target behavior is to skip the recipe if the file exists.
* Add Phase 3 bit-equality regression harness
8 configurations covering single-pop neutral, single-pop with size
changes, two-pop with constant migration, two-pop with split,
two-pop with split + migration, three-pop combo, admixture, and
ancient samples. Each config runs both binaries with identical
seeds and diffs outputs ignoring the line-1 command-line echo.
* Add initializeShapesFromGlobals to seed shape state from globals
Copies currentSize[] -> popShape[] and migMatConst[][] -> migShape[][]
as SHAPE_CONSTANT shapes anchored at t=0. Will be called from
initialize() in the next task to seed per-replicate shape state.
* Call initializeShapesFromGlobals from per-replicate initialize()
Seeds popShape[] and migShape[][] from the current parameter values
at the start of each replicate. Shape state is computed but not
yet read; bit-equality preserved across the regression grid.
* Add shapes.c to test targets that compile discoalFunctions.c
Phase 3 Task 4 wired initializeShapesFromGlobals into initialize()
in discoalFunctions.c, but the unit-test build recipes that compile
discoalFunctions.c (test_node, test_event, test_node_operations,
test_mutations, test_trajectory, etc.) did not include shapes.c
in their source lists. This left make run_tests link-broken with
undefined references to initializeShapesFromGlobals (and would
have done the same for sizeAt as more shapes accessors are wired
into the engine code).
Add $(SRC_CORE)/shapes.h $(SRC_CORE)/shapes.c to the dependency
list and $(SRC_CORE)/shapes.c to the gcc invocation for every
test recipe that compiles discoalFunctions.c.
* Update popShape on 'n' event firing
When -en or other size-change events fire, the new size is now
written to both currentSize[popID] and popShape[popID]. Future
sizeAt() reads at any time after the event get the new value.
Bit-equal regression preserved across the configuration grid.
This commit lands BEFORE the cRate/migration sizeAt/migAt refactor
(originally Tasks 5/6, re-ordered to come after this one) so that
the bit-equality regression stays green at every step: refactoring
the rate computation to read from popShape only makes sense once
popShape is being written by the event handler.
* Read coalescent rate denominator via sizeAt in neutralPhaseGeneralPopNumber
cRate[i] now uses sizeAt(i, currentTime) instead of sizeRatio[i].
Under SHAPE_CONSTANT, sizeAt returns popShape[i].anchor_value
which is kept in sync with currentSize[i] by both
initializeShapesFromGlobals (initial seed) and the 'n' event
handler (mutations on -en events). Output is bit-equal across
all 8 regression configurations including single_pop_size_change.
* Read migration rate via migAt in neutralPhaseGeneralPopNumber
Replaces all migMat[i][j] reads in the function (row-sum mRate[i]
and per-pair event sampling) with migAt(i, j, currentTime). Under
SHAPE_CONSTANT, migAt returns migShape[i][j].anchor_value which
equals migMatConst[i][j] (and hence the migMat[i][j] seed); output
is bit-equal across the regression grid.
* Plan: Phase 4 NHPP sampler and exponential growth shape
6 tasks, ~30 steps. Adds SHAPE_EXPONENTIAL via:
- allShapesConstant() predicate gating inner-loop dispatch
- drawNHPPWaitingTime helper for non-constant path
- 'g' event type + dispatch handler
- CLI flags -eg / -eG mirroring ms
Bit-equality preserved for all-CONST configs by dispatching on
the predicate and keeping the existing Exp(total) sampler unchanged
in that branch. msprime parity validation deferred to Phase 4b.
* Add allShapesConstant predicate for sampler dispatch
Returns 1 iff every active shape (popShape[i] for i<npops and
migShape[i][j] for i,j<npops) has type SHAPE_CONSTANT. Used by
the inner-loop sampler in Phase 4 to keep the Exp(total) fast
path for the common all-constant case while routing non-constant
shapes through the NHPP per-component sampler.
* Add NHPP per-component sampler dispatch for non-constant shapes
neutralPhaseGeneralPopNumber now branches on allShapesConstant():
- All-constant: existing Exp(total_rate) sampler. Bit-equal preserved.
- Otherwise: drawNHPPWaitingTime selects minimum across per-component
closed-form draws (drawWaitingTimeSize / drawWaitingTimeMig for
coalescence and migration; Exp(total) for recomb/gc which are
time-constant). The winning event class and target are returned
to the caller for direct dispatch.
No SHAPE_EXPONENTIAL / SHAPE_LINEAR shapes are wired into the
runtime yet, so this commit exercises only the constant path.
The non-constant path becomes active in Task 3 once the 'g'
event handler starts setting popShape.type = SHAPE_EXPONENTIAL.
* Add 'g' event handler for SHAPE_EXPONENTIAL transitions
When a 'g' event fires, popShape[popID] is set to (EXPONENTIAL,
sizeAt(popID, t), alpha, t). The anchor_value uses sizeAt at the
event time so the size is continuous across the shape-change
boundary (msprime forward-time convention; spec section 4.5
Convention 1). The alpha rate_param comes from events[j].popnSize.
The case body otherwise mirrors 'n' - runs the inter-event
interval through neutralPhase or sweepPhase as appropriate.
No CLI / importer flag emits 'g' events yet; that comes in Task 4.
* Add -eg and -eG CLI flags for exponential-growth events
-eg time popID alpha emits one 'g' event setting popShape[popID]
to SHAPE_EXPONENTIAL with the given growth rate at the given time
(in 2N units, multiplied by 2.0 for internal 4N representation
matching -en convention).
-eG time alpha emits one 'g' event per active population, all
with the same alpha. -p must come before -eG for npops to be set.
alpha is the per-generation forward-time growth rate per the
msprime convention. The CLI matches ms's -eg / -eG flags.
* Add -eg smoke test to confirm EXP shape is exercised
Three configs at the same RNG seed: no demography, -eg with
alpha=50 at t=0.5, and -en with size change at t=0.5. The smoke
test confirms -eg differs from both no-demography and -en --
sanity check that the new event type is actually wired into the
simulation engine and not silently no-op'ing.
* Plan: Phase 4b msprime parity for SHAPE_EXPONENTIAL
6 tasks, ~24 steps. Validates Phase 4 EXP support against msprime
1.4.1 via Bonferroni-corrected KS at alpha=0.01 on segregating
sites, pi, Tajima's D, and SFS for single-pop and two-pop split
configurations.
Includes a parity_utils library with parameter conversions
(discoal CLI 2N -> msprime gen, alpha 4N -> per-gen) and a
runtime sanity check (constant-N pi within 10% tolerance) that
validates conversions before any actual EXP test runs.
* Add Phase 4b msprime parity utility library
parity_utils.py provides:
- ms-format parser for discoal stdout (skips DEBUG lines / cmdline echo)
- msprime tree-sequence to ms-format converter
- Per-replicate stats: segregating sites, pi, Tajima's D, folded SFS
- Parameter conversions (discoal CLI 2N time -> msprime generations,
alpha 4N-scaled -> per-generation, theta/rho -> per-gen-per-site)
- Convention dataclass + three trial conventions (A/B/C) for the
ploidy / popsize choice
- check_constant_parity() trials all three against a no-demography
baseline and reports per-convention rel_err on mean pi
- validate_active_convention() guards EXP-parity tests against
conversion regression
Trial run (n=6, theta=5.0, rho=5.0, L=1000, Ne=1e6, reps=200):
discoal mean pi = 4.9710
A: ploidy=1, popsize=Ne msprime mean pi = 2.4700, rel_err = 0.503
B: ploidy=2, popsize=Ne msprime mean pi = 4.8997, rel_err = 0.014
C: ploidy=1, popsize=2*Ne msprime mean pi = 5.0860, rel_err = 0.023
ACTIVE_CONVENTION locked to CONVENTION_B (ploidy=2, popsize=Ne)
with rel_err 0.014, comfortably within the 10% tolerance.
validate_active_convention() PASSes.
* Add single-pop EXP growth msprime parity test
Runs discoal -eg 0.5 0 50 and msprime with matched
add_population_parameters_change at the converted generation time
and per-generation growth rate. Computes segregating sites, pi,
Tajima's D, and folded SFS across replicates, then runs Bonferroni-
corrected Kolmogorov-Smirnov / chi-squared tests at alpha = 0.01.
Uses ACTIVE_CONVENTION from parity_utils to set msprime's ploidy
and population_size; calls validate_active_convention at startup
to guard against conversion regression.
Usage: python3 test_single_pop_exp.py [reps]
Default reps=1000; use 100 for a quick smoke run.
* Fix parity_utils sample-count and mutation-model bugs
_msp_constant_parity_pi was producing 2x more samples than
discoal under Convention B (ploidy=2) because msprime's
samples=n means n diploid individuals = 2n haploid samples,
but discoal's n is already haploid. Now passes n // ploidy.
Also switched the mutation model to BinaryMutationModel with
discrete_genome=False to match discoal's infinite-sites binary
semantics; the msprime default JC69 + discrete_genome=True
produces multi-allelic sites that break SFS computation downstream.
The blind spot was masked by pi being roughly sample-size-
invariant, so validate_active_convention PASSed at 1.4% rel_err
even with the 2x sample mismatch. The fixes are now applied
consistently across parity_utils and the EXP-parity tests.
* Add two-pop split + EXP growth msprime parity test
Two populations with EXP growth in pop0 (rate alpha=50, growth
onset at t=0.3 in 2N units) and a split at t=1.0. Compares
discoal -p 2 n n -eg t pop alpha -ed t_split 0 1 against
msprime Demography with population_split + parameters_change.
Uses ACTIVE_CONVENTION (B: ploidy=2, popsize=Ne) and the
n // ploidy sample-count + BinaryMutationModel fixes from the
single-pop test. n_per_pop = 6 (divisible by ploidy=2).
Bonferroni-corrected KS at alpha = 0.01 across ss, pi, Tajima's D.
* Add Phase 4b orchestrator and run results (PASS)
run_parity.sh runs both single-pop and two-pop EXP parity tests
at REPS=1000, logs to results/analysis.txt, and exits with
overall PASS/FAIL.
Result at REPS=1000: OVERALL PASS.
- Single-pop EXP: 4/4 comparisons within Bonferroni alpha=2.5e-3.
Smallest p = 6.7e-3 (SFS chi-squared).
- Two-pop split + EXP: 3/3 comparisons within Bonferroni alpha=3.3e-3.
Smallest p = 0.31.
The committed analysis.txt records the outcome.
* Spec addendum: Phase 4b msprime parity result (PASS)
* Fix discoal->msprime growth-rate conversion (off by factor of 2)
discoal_alpha_to_msp_growth was dividing by 4*Ne, but the right
factor is 2*Ne. discoal stores alpha as the rate_param of an
exponential shape evaluated in *internal* time units
(size(t)=anchor*exp(-alpha*(t-t0)), see shapes.c). The pair
coalescent rate in neutralPhase is n*(n-1)/2 per internal unit,
so 1 internal unit = 2*Ne_diploid generations, and the per-gen
growth rate satisfying alpha*t_internal = g*t_gen is g=alpha/(2*Ne).
The CLI-time conversion (t_cli * 4 * Ne) is unaffected and remains
correct: discoal multiplies user-supplied event times by 2.0 to enter
internal units, contributing an extra factor of 2 vs the alpha case.
The bug didn't affect constant-parity (no growth events) and was
masked at the original single-point EXP test (alpha=50, t_cli=0.5,
where deeper growth onset weakens the dependence on the exact
exponential rate). The broadened sweep caught it via the a50_t0.1
cell (small alpha, recent growth onset) where the per-gen growth
rate strongly determines diversity.
Plus broaden single-pop EXP test to a 3x3 (alpha, eg_time) sweep
with 11 statistics per cell (ss, pi, tajD, wtheta, hapdiv, nhap,
plus per-bin SFS), Bonferroni-corrected KS at alpha=0.01 across
99 comparisons. With the alpha fix in place, the sweep PASSes at
REPS=5000.
parity_utils now has watterson_theta, haplotype_diversity, and
num_haplotypes helpers. collect_stats returns these new arrays
in addition to the prior ss/pi/td/sfs.
* Spec addendum: broadened single-pop EXP parity + alpha-conversion fix
* Plan: Phase 5 sweep wiring with closed-form detSweepFreqGeneral
12 tasks. Wires sizeAt into proposeTrajectory,
sweepPhaseEventsConditionalTrajectory, and
sweepPhaseEventsGeneralPopNumber. Adds integratedSizeRatio
accessor + detSweepFreqGeneral that reduces to detSweepFreq
exactly under constant shapes. Bit-equal regression for
SHAPE_CONSTANT sweep configs (6 cells). Smoke test for
sweep + EXP. Removes detSweepFreqEuler and --det-sweep-mode
flag per the spec revision.
* Add integratedSizeRatio for sweep-deterministic-mode integration
Returns int_{t0}^{t0+T} sizeAt(popID, s) ds in closed form for all
three shape types. Used by Phase 5's detSweepFreqGeneral to compute
the integrated selection coefficient A(tau) under time-varying N
without numerical integration.
* Add detSweepFreqGeneral for time-varying alpha_eff sweeps
Generalizes detSweepFreq to take A = int alpha_eff(s) ds rather
than alpha and tau separately. Reduces to detSweepFreq exactly when
A = alpha*tau (constant). Boundary condition anchored to the
user-specified alpha (epsilon = 0.05/alpha, A_ts = -2*log(epsilon)).
Unit test verifies bit-equivalence to detSweepFreq for constant
alpha across a range of tau values.
* Add discoal_pre_phase5 Makefile recipe for sweep regression
* Add Phase 5 sweep bit-equality regression harness
6 sweep configurations covering deterministic / stochastic-forward /
neutral-stochastic modes, single-pop / two-pop with size change,
varying alpha and tau.
* Save/restore popShape across proposeTrajectory; track 'g' events
proposeTrajectory now mutates popShape as it walks events forward
(in addition to the existing currentSizeRatio mutation), saving
and restoring the global state around the function call so the
caller's view of popShape is unchanged. Future tasks read sizeAt
inside the trajectory loop, which depends on this state being
correct during the walk.
The existing currentSizeRatio updates remain as a backward-
compatible cache for the trajectory loop's inner step computations,
to be removed in Phase 5 Task 6 once those reads are migrated to
sizeAt.
* Read size via sizeAt in proposeTrajectory inner loop
Inner trajectory loop now reads sizeAt(0, currentTime + ttau)
each step, replacing the per-epoch currentSizeRatio scalar.
Under SHAPE_CONSTANT (the only shape currently produced by the
events walk), sizeAt returns the anchor_value which equals the
just-set currentSizeRatio, preserving bit-equality.
Under SHAPE_EXPONENTIAL the size now varies per-step within
an epoch — the intended behavior. Phase 5 Task 8 will replace
the deterministic-mode detSweepFreq call with the closed-form
detSweepFreqGeneral for the variable-size case.
* Read sizes via sizeAt in sweepPhaseEventsConditionalTrajectory
Replaces all sizeRatio[i] reads in the function body with
sizeAt(i, cTime + ttau). Under SHAPE_CONSTANT, sizeAt returns
popShape[i].anchor_value which equals the prior sizeRatio[i]
seed; output is bit-equal.
* Dispatch deterministic-mode sweep on shape constancy
For SHAPE_CONSTANT-only configs, keep the existing
detSweepFreq(ttau, alpha*sr_now) path bit-equal with pre-Phase-5.
For non-constant shapes, use detSweepFreqGeneral with
A_now = A_prev + alpha*integratedSizeRatio(0, t-dt, dt)
incrementally tracked.
Mirrors the Phase 4 NHPP dispatch pattern in the neutral phase:
common case stays exact; non-constant case uses the time-varying
closed form with no truncation error.
* Same constant/general dispatch in sweepPhaseEventsGeneralPopNumber
Deterministic-mode sweep in the recurrent-sweep code path now
routes through the same allShapesConstant dispatch as
proposeTrajectory: constant -> detSweepFreq (bit-equal),
non-constant -> detSweepFreqGeneral with incremental A tracking.
Plus replaces sizeRatio[i] reads in the rate-computation block
with sizeAt(i, cTime + ttau) to support time-varying shapes
during a sweep.
The existing detSweepMode (Euler) split is preserved within the
constant branch for Phase 5 Task 11 to remove cleanly.
* Add sweep + EXP growth smoke test
Confirms discoal runs to completion when -wd/-ws and -eg are
combined. Sanity check that the time-varying alpha_eff path
in proposeTrajectory and sweepPhaseEvents* doesn't crash on
typical parameters.
* Remove detSweepFreqEuler and --det-sweep-mode
Per the spec revision, the deterministic-sweep ODE under time-
varying N has a closed-form solution via detSweepFreqGeneral; no
Euler step is needed. The Phase 2 Q1 verification was conclusive
(detSweepFreq != Euler under constant N), and the Euler path was
never used at runtime after Phase 5 wired in the closed-form
general dispatch.
Deleted:
- detSweepFreqEuler function and unit test
- detSweepMode global and the --det-sweep-mode CLI flag
- Q1 verification harness scripts (test/parity/q1_*.sh, *.py)
- Residual if (detSweepMode == 0) splits in discoalFunctions.c
Kept:
- The Q1 spec addendum (historical context for why we chose the
closed-form approach over Euler)
- test/parity/q1_results/analysis.txt (the recorded result)
* Plan: Phase 6 SHAPE_LINEAR engine wiring
5 tasks. Adds CLI flags -el / -eL and 'l' event handler for
linear-growth events. Mirrors Phase 4's 'g' event pattern.
Bit-equality preserved for SHAPE_CONSTANT configs; smoke test
confirms -el produces output distinct from -en, -eg, and
no-demography. msprime parity deferred to Phase 8 if needed
(msprime lacks native linear growth).
* Add -el and -eL CLI flags for linear-growth events
-el time popID gamma emits one 'l' event setting popShape[popID]
to SHAPE_LINEAR with the given linear growth rate at the given
time (in 2N units, multiplied by 2.0 for internal 4N representation
matching -en / -eg convention).
-eL time gamma emits one 'l' event per active population. -p must
come before -eL for npops to be set.
gamma is the per-generation forward-time linear rate of size
change, msprime/forward-time convention; positive gamma means
the population grew forward (was smaller in the past).
The 'l' event handler in the main event-dispatch switch is
added in Task 2; until then 'l' events are parsed and added
to events[] but silently fall through.
* Add 'l' event handler for SHAPE_LINEAR transitions
When an 'l' event fires, popShape[popID] is set to (LINEAR,
sizeAt(popID, t), gamma, t). Continuity preserved across the
shape boundary: anchor_value uses sizeAt at the event time,
matching the 'g'-event pattern from Phase 4.
The case body otherwise mirrors 'g' verbatim — runs the
inter-event interval through neutralPhase or sweepPhase as
appropriate based on activeSweepFlag and recurSweepMode.
* Track 'l' events in proposeTrajectory events walk
Mirrors the 'g'-event tracking added in Phase 5 Task 5. When
proposeTrajectory walks future events to predict the trajectory
shape, 'l' events set popShape to SHAPE_LINEAR. The same
save/restore pattern (memcpy at entry/exit) keeps the global
state unchanged for the caller.
* Add -el smoke test to confirm SHAPE_LINEAR is exercised
Four configs at the same RNG seed: no demography, -el gamma=0.5
at t=0.5, -en size change at t=0.5, and -eg alpha=50 at t=0.5.
The smoke test confirms -el produces distributionally different
output from all three reference configs — sanity check that the
LINEAR shape path is wired into the simulation engine and not
silently treated as constant or exponential.
* Plan: Phase 7 migration shapes + importer rewrite (closes #82)
8 tasks. Adds -em/-eM CLI flags and 'm' event handler. Rewrites
the demes importer: piecewise-constant migration matrices per
interval, emit 'm' events at boundaries, write migMatConst[]
directly for the t=0 interval, lift exp/linear epoch rejection
and emit 'g'/'l' events. Removes the back-derivation hack at
discoalFunctions.c:222-261. Validates via msprime parity on the
issue #82 fixture.
* Add -em and -eM CLI flags for time-varying migration events
-em time srcPop dstPop rate emits one 'm' event setting
migShape[src][dst] to SHAPE_CONSTANT with the given migration
rate at the given time (in 2N units, multiplied by 2.0 for
internal 4N representation matching -en / -eg / -el convention).
-eM time rate emits one 'm' event per off-diagonal pair, all
with the same rate. -p must come before -eM for npops to be set.
Event field convention: events[].popID is the destination
population, events[].popID2 is the source -- matches the legacy
'M' event convention used by the demes importer. The runtime
handler in Task 2 indexes migShape[src][dst] accordingly.
Lowercase 'm' avoids colliding with the legacy 'M' events emitted
by the unrewritten parts of the importer; once Task 4 rewrites
the importer, no 'M' events are emitted any longer.
* Add 'm' event handler for time-varying migration
When an 'm' event fires, migShape[src][dst] is set to
(SHAPE_CONSTANT, new_rate, 0, t_event). src is events[j].popID2,
dst is events[j].popID -- matches the CLI parser convention from
Task 1 and the legacy 'M' event convention.
The case body otherwise mirrors 'g' / 'l' -- runs the inter-event
interval through neutralPhase or sweepPhase as appropriate.
The neutral-phase sampler already reads migAt(src, dst, t) which
indexes migShape, so the new rate takes effect immediately for
subsequent inter-event waiting-time draws.
* Add -em smoke test to confirm migShape mutation is exercised
Three configs at the same RNG seed: no migration, constant
migration via -m, and -em turning off migration mid-simulation.
The smoke test confirms -em produces distributionally different
output from the constant-migration baseline -- sanity check that
the new event type actually mutates migShape and the inner-loop
sampler sees the change.
* Importer: emit interval-based 'm' events for demes migration windows
Replaces the broken paired-'M' event emission (issue #82) with a
proper interval-based scheme. For each unidirectional migration
window in graph->migrations:
- The active rate at t=0 (the present) is summed across windows
covering t=0 and written directly to migMatConst[src][dst].
- For each boundary time (window start or end > 0), if the
active rate after the boundary (going further into the past)
differs from the rate before, emit an 'm' event at that
boundary setting the rate to the post-boundary value.
The runtime now correctly tracks time-varying migration rates,
including the issue #82 fixture (multi-window migration where
some demes don't exist across the full timespan).
Note: this commit replaces ONLY the migration handling, not the
exp/linear epoch rejection (Task 5) or the back-derivation hack
in discoalFunctions.c (Task 6). After Task 6, the runtime no
longer needs to back-derive anything from 'M' events; the
importer writes migMatConst directly.
* Importer: emit 'g'/'l' events for exponential / linear demes epochs
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.
* Remove back-derivation hack from initialize()
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.
* Add msprime parity test for the issue #82 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.
* Spec addendum: Phase 7 issue #82 closure (pi parity PASS, ss residual documented)
* Add branch-mode parity diagnostic for the issue #82 fixture
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).
* Fix currentSize bleed between replicates after 'n' size-change events
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).
* Spec addendum: document currentSize replicate-bleed bug and its fix
* Simplify currentSize restore: trim comments, drop NULL guard, relocate 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.
* Docs: add new -em/-eM/-eg/-eG/-el/-eL flags; correct demes time conversion
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.
* Fix demes time_units=years handling; add byte-equality regression
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.
* Stop tracking docs/superpowers/ (internal planning grist)
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.
* ci: guard Unity clone in Memory Leak job (mirror Basic Functionality)
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.
* Stop tracking test/parity/ (parity scripts kept local for development)
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.
* comment clarity
* Preflight: reject configurations that drive a population size to zero
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.
* Strip plan/phase references from in-code comments
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.
* Abort instead of silently zeroing on unknown shape type
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.
* Document integrated-hazard math; clamp linear migration past zero crossing
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.
* migration event checking
* Note YAML config as a migMatConst source in initialize() comment
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.
* Add end-to-end preflight validation regression suite
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.1 parent e08c756 commit 0ad2d5a
19 files changed
Lines changed: 2466 additions & 269 deletions
File tree
- .github/workflows
- docs
- src/core
- testing
- test/unit
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
128 | 128 | | |
129 | 129 | | |
130 | 130 | | |
131 | | - | |
132 | | - | |
133 | | - | |
134 | | - | |
| 131 | + | |
| 132 | + | |
| 133 | + | |
| 134 | + | |
| 135 | + | |
| 136 | + | |
135 | 137 | | |
136 | 138 | | |
137 | 139 | | |
| |||
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
99 | 99 | | |
100 | 100 | | |
101 | 101 | | |
| 102 | + | |
| 103 | + | |
| 104 | + | |
| 105 | + | |
| 106 | + | |
| 107 | + | |
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
52 | 52 | | |
53 | 53 | | |
54 | 54 | | |
55 | | - | |
| 55 | + | |
56 | 56 | | |
57 | | - | |
| 57 | + | |
58 | 58 | | |
59 | 59 | | |
60 | 60 | | |
61 | | - | |
| 61 | + | |
62 | 62 | | |
63 | | - | |
| 63 | + | |
64 | 64 | | |
65 | 65 | | |
66 | | - | |
| 66 | + | |
67 | 67 | | |
68 | | - | |
| 68 | + | |
| 69 | + | |
| 70 | + | |
| 71 | + | |
| 72 | + | |
| 73 | + | |
| 74 | + | |
| 75 | + | |
| 76 | + | |
| 77 | + | |
| 78 | + | |
| 79 | + | |
| 80 | + | |
| 81 | + | |
| 82 | + | |
| 83 | + | |
| 84 | + | |
| 85 | + | |
| 86 | + | |
| 87 | + | |
| 88 | + | |
| 89 | + | |
| 90 | + | |
| 91 | + | |
| 92 | + | |
| 93 | + | |
| 94 | + | |
| 95 | + | |
| 96 | + | |
| 97 | + | |
| 98 | + | |
| 99 | + | |
| 100 | + | |
| 101 | + | |
| 102 | + | |
| 103 | + | |
69 | 104 | | |
70 | 105 | | |
71 | | - | |
| 106 | + | |
72 | 107 | | |
73 | | - | |
| 108 | + | |
74 | 109 | | |
75 | 110 | | |
76 | 111 | | |
| |||
238 | 273 | | |
239 | 274 | | |
240 | 275 | | |
241 | | - | |
| 276 | + | |
242 | 277 | | |
243 | 278 | | |
244 | 279 | | |
245 | | - | |
| 280 | + | |
246 | 281 | | |
247 | 282 | | |
248 | 283 | | |
249 | 284 | | |
250 | 285 | | |
251 | 286 | | |
252 | | - | |
| 287 | + | |
| 288 | + | |
| 289 | + | |
| 290 | + | |
| 291 | + | |
| 292 | + | |
| 293 | + | |
| 294 | + | |
| 295 | + | |
253 | 296 | | |
254 | 297 | | |
255 | 298 | | |
256 | | - | |
| 299 | + | |
257 | 300 | | |
258 | 301 | | |
259 | 302 | | |
| |||
268 | 311 | | |
269 | 312 | | |
270 | 313 | | |
271 | | - | |
| 314 | + | |
272 | 315 | | |
273 | 316 | | |
274 | 317 | | |
275 | | - | |
| 318 | + | |
276 | 319 | | |
277 | 320 | | |
278 | | - | |
| 321 | + | |
279 | 322 | | |
280 | 323 | | |
281 | 324 | | |
282 | 325 | | |
283 | | - | |
| 326 | + | |
284 | 327 | | |
285 | 328 | | |
286 | 329 | | |
287 | 330 | | |
288 | 331 | | |
289 | 332 | | |
290 | | - | |
| 333 | + | |
291 | 334 | | |
292 | 335 | | |
293 | 336 | | |
| 337 | + | |
| 338 | + | |
294 | 339 | | |
295 | 340 | | |
296 | 341 | | |
| |||
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
39 | 39 | | |
40 | 40 | | |
41 | 41 | | |
42 | | - | |
43 | | - | |
| 42 | + | |
| 43 | + | |
| 44 | + | |
| 45 | + | |
| 46 | + | |
| 47 | + | |
44 | 48 | | |
| 49 | + | |
| 50 | + | |
45 | 51 | | |
46 | 52 | | |
47 | 53 | | |
| |||
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
22 | 22 | | |
23 | 23 | | |
24 | 24 | | |
| 25 | + | |
| 26 | + | |
25 | 27 | | |
26 | | - | |
| 28 | + | |
| 29 | + | |
27 | 30 | | |
28 | 31 | | |
29 | 32 | | |
30 | 33 | | |
31 | | - | |
32 | | - | |
33 | 34 | | |
34 | 35 | | |
35 | 36 | | |
| |||
40 | 41 | | |
41 | 42 | | |
42 | 43 | | |
43 | | - | |
44 | | - | |
45 | | - | |
46 | | - | |
47 | | - | |
| 44 | + | |
| 45 | + | |
| 46 | + | |
| 47 | + | |
| 48 | + | |
| 49 | + | |
| 50 | + | |
| 51 | + | |
| 52 | + | |
| 53 | + | |
| 54 | + | |
| 55 | + | |
| 56 | + | |
48 | 57 | | |
49 | 58 | | |
50 | 59 | | |
| |||
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
120 | 120 | | |
121 | 121 | | |
122 | 122 | | |
123 | | - | |
124 | | - | |
125 | | - | |
126 | | - | |
127 | | - | |
128 | | - | |
129 | | - | |
130 | | - | |
| 123 | + | |
| 124 | + | |
| 125 | + | |
| 126 | + | |
| 127 | + | |
| 128 | + | |
| 129 | + | |
| 130 | + | |
131 | 131 | | |
132 | 132 | | |
133 | 133 | | |
| |||
0 commit comments