Conversation
…ree matcher Building a KRegCCD was dominated by computing the per-clade escape reserve, and within that by the boundary-4 count N_2 (computeReg -> countNj). The generic enumerate is O(m^3) in the clade's observed-subclade count m (choose 3 parts, derive the 4th as the complement), so on the large near-root clades it burned the whole op-budget and, on real 320-taxon data, made construction take days. Replace the boundary-4 path with a meet-in-the-middle count: a 4-part boundary is two disjoint observed pairs whose unions are complementary in C, so hashing pairs by union turns the search into O(m^2). The matcher is allocation-free -- pair unions go into a reused scratch bitset and pairs live in flat int[] arrays with a chained open hash, instead of a HashMap<BitSet>-with-int[]-per-pair. A boundary is then "this pair + the pair whose union is C\union(i,j), with k>j", emitting each boundary once (i<j<k<l). Exact: returns the identical count as the old enumerate for both FLAT (sum of pathcounts) and SHARED (admissible-boundary indicator), validated against the brute-force triple-scan on every clade of the coal-n320 real data. Same enumOps/opsBudget guard, so pathological clades still degrade to the capped N_2=0. precomputeReserves on coal-n320 (parallel): 300 trees 260s->46s, and the change also computes N_2 exactly on the big clades the op-budget previously capped to 0. All KRegCCD/MRegCCD tests pass.
Follow-up to the meet-in-the-middle boundary-4 matcher. Two further exact speedups to computeReg, the reserve hot path: 1. Find a pair's complement by weighted-sum arithmetic instead of materialising the complement bitset and re-hashing it. getCladeInBits()'s hashCode weights word w by (w+1), so weightedSum is additive over disjoint sets: a pair's union sum is the sum of its parts' precomputed sums, and the complement's sum is weightedSum(C) minus it. Pairs are hashed by that long sum; a chain hit is confirmed exactly (immune to sum wraparound) by the cardinality identity |Pi|+|Pj|+...==|C| plus disjointness, which for subclades of C is equivalent to tiling C. Removes the per-pair complement clear/or/andNot + hashCode that dominated the matcher. 2. Compute boundary-3 (N_1) and boundary-4 (N_2) in a single disjoint-pair pass (countN1N2), so the O(m^2) pair enumeration is scanned once, not twice. N_1's complement (a single observed subclade) is looked up in a sum-keyed hash of the subclades by the same arithmetic. This retires the generic O(m^2) enumerate for the default reserve depth k=2; deeper orders and the earlyExit path keep the old code. Exact: N_1 and N_2 match the previous enumerate/meet-in-the-middle counts on every clade of the coal-n320 real data (both FLAT and SHARED); all KRegCCD/MRegCCD tests pass. precomputeReserves on coal-n320: 300 trees 46s->24s (and ~10 min originally).
…ration) getEntropyRecursive's blueContribution still enumerated the boundary-3/4 regions with the old O(m^3) blueWalk, making entropy the workflow bottleneck (144s at just 30 trees on coal-n320). Replace it with blueContributionFast: the same weighted-sum disjoint-pair enumeration used by countN1N2, accumulating the identical per-boundary entropy terms (boundaryMass * sum of parts' forced-red entropy, plus the SHARED log-pathcount self-information). Deeper orders (>=5, only reached when boundaries 3 and 4 are empty) keep the old blueWalk. Exact: entropy value unchanged (148.74 nats at 30 trees) and the KRegCCDEntropyTest brute-force/Monte-Carlo checks pass. Entropy 144s -> 2.9s at 30 trees.
Reserve depth k=1 (eps solved from N_1 alone, tail=0) previously fell through to the slow generic enumerate. Route it through the fast weighted-sum pair pass too: countN1N2 gains a computeN2 flag that, when false, skips the boundary-4 stash and match entirely, and computeReg's fast path now covers reserveBoundary 3 (k=1) as well as 4 (k=2). k=1 gives a coarser but much cheaper reserve -- N_2 only shifts eps by <=7% (its eps^2 term is second order). On coal-n320, 1000 trees: construct 150s->29s, entropy 478s->113s. MAP tree and reserves match k=2 to ~1e-2 nats (all-red MAP unaffected by the dropped order); self-consistent (getMaxLogTreeProbability == logProb(MAPtree)). All KRegCCD/MRegCCD tests pass.
getEntropyRecursive was a single-threaded loop over clades, so entropy was the workflow bottleneck (~478s at 1000 trees on coal-n320) even though it does the same boundary enumeration as the already-parallel precomputeReserves. Both entropyRedForced(c) and entropyFree(c) read only strictly-smaller clades (partition children and blue-region boundary parts), so once each smaller size level is memoised, all clades of a given size are independent. Bucket clades by size and compute each level with a parallelStream into concurrent memo maps -- the same level-by-level parallelism computeReg gets. Per-clade arithmetic is unchanged, so the entropy value is identical to serial (KRegCCDEntropyTest still passes). Entropy at 300 trees: ~29s on 8 cores.
getMAPTree's discounted DP needs every clade's reserve, but solved them lazily inside its serial size-ordered loop, so a getMAPTree() on a cold model paid the whole per-clade O(m^2) reserve cost single-threaded -- roughly 4 h on 1000 trees of coal-n320 thinned over the chain (14,338 clades), which reads as a hang. computeDiscountedMap now calls precomputeReserves() up front, as getEntropyRecursive already did. Verified to give an identical MAP tree and log-probability; this changes cost only. Sampling drew one tree at a time. AbstractCCD.sampleTrees(n) fans the draws out instead: per-draw seeds are taken from the shared RNG before the fan-out, so results stay reproducible under setRandom however the draws are scheduled; runningInnerIndex and the sampling RNG become thread-local; and a prepareForSampling() hook lets KRegCCD warm its reserves so the draws are pure reads. 7.8x on 100 draws (62.0 s -> 8.0 s, 300 thinned trees, 8 cores). That exposed a latent bug: KRegCCD kept its blue-region reservoir-sampling state (boundarySeen / boundaryWeightSeen / boundaryPick) in instance fields, which concurrent draws corrupted (ArrayIndexOutOfBoundsException out of Tree.listNodes). The state is live only within a single sampleBoundary call, so it moves into a call-local BoundaryReservoir.
Compares, at the resolution of small induced subtrees, the topology distributions implied by the empirical posterior sample and by CCD0, CCD1 and (optimised) regCCD built from it. The purpose is to establish that CCD1 is an adequate basis to build on: it conditions each split on its parent clade and so carries the posterior's correlation structure, which regularisation then smooths to give full support. CCD0 is the negative control (it factorises over clades and cannot represent those correlations); regCCD is CCD1 plus smoothing and is expected to sit slightly further from the posterior than CCD1 does. The question is whether CCD1's deviation is small relative to the posterior's own sampling noise, not which model wins on total variation -- so the tool reports a split-half noise floor alongside the model columns. The marginal probability of a full topology is not Monte-Carlo-estimable (almost every sampled tree is unique), but the induced distribution over a fixed small taxon subset is, and it interpolates between single-clade marginals and the full topology. Trees are restricted to each subset by a non-destructive induced-topology walk, O(k*height) with no tree copy. Supports sampled or exact (sum-product DP over the clade DAG) induced probabilities for CCD0/CCD1, and a calibration mode that summarises the posterior alone.
…gical uncertainty The induced-subtree comparison can only measure something on subsets where the POSTERIOR itself spreads over more than one induced shape; the downstream analysis discards the rest. Nothing checked that any such subset existed before doing the expensive part, so a degenerate posterior bought a full run -- regCCD parameter selection plus three 50,000-tree sampling passes -- and a table the analyser then threw away entirely. Found on tornabene-2016: all 1801 post-burn-in trees have the same topology, so 0 of 1000 subsets were informative at k = 4, and no larger k would have helped either. That cost 12 minutes; on a larger dataset it would have been hours. The empirical tally already computes what is needed, so the check is free. Model construction is now deferred until after that tally and gated on the informative-subset count (-minSubsets, default 1; 0 forces the run). The tool also reports the number of distinct posterior topologies, which is the one line that makes the diagnosis obvious. When the gate trips it writes the empirical columns and leaves the model columns at 0, as calibration mode already does. Deferring construction does not change any result: the three model RNG seeds are still drawn at the same point in the sequence, and construction never touches that Random -- so subset selection and model sampling are bit-identical to before. Gated, tornabene now exits in 14 seconds.
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
3 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.
Summary
O(m^2) meet-in-the-middle + sum-arithmetic pass, and parallelise reserve
solving, entropy, and tree sampling across clades/draws. On real 320-taxon
data this takes construction from ~days to ~24s (300 trees) and entropy from
478s to 113s (1000 trees), with identical results on every existing test
(KRegCCD/MRegCCD/entropy suites all pass, brute-force-validated on real data).
state lived in instance fields and got corrupted by concurrent draws.
muto 0.005.ccd.tools.SubtreeMarginal: compares induced-subtree topologydistributions between the empirical posterior and CCD0/CCD1/regCCD, to
validate that CCD1 adequately captures the posterior's correlation
structure. Skips model construction/sampling when the posterior has no
topological uncertainty on any subset (avoids a 12-minute run for a
discarded result on a degenerate real dataset; now exits in 14s).
Maven release with UnrecognizedPropertyException Exception beast3#117.