Skip to content

Drop top n - #12

Merged
jeromekelleher merged 6 commits into
tskit-dev:mainfrom
jeromekelleher:drop-top-n
Oct 28, 2026
Merged

jeromekelleher merged 6 commits into
tskit-dev:mainfrom
jeromekelleher:drop-top-n

Conversation

@jeromekelleher

Copy link
Copy Markdown
Member

No description provided.

Add a tgp_samples_map workload to evaluate.py: chr20p sample matching
with per-site r from the HapMap II GRCh38 chr20 map, mu = 1e-20, 1e-8,
1e-4, 4 samples each (about 5 minutes on 2 workers). Its grid is the
strategies kept after the constant-rate runs, floor 1e-8, and thresholds
at c times the median switch odds (prune_qmed, relative_qmed).
--workloads and --results select workloads and the results directory.

Under the map, q spans 7e-18 to ~1e-9 (median 1.9e-13). Results:

- floor 1e-8 (below q with constant r) is catastrophic: Viterbi losses
  1e5-1e6 on every query, forward-backward logL errors 5-29, TV to 0.89.
- lossless_prune at the smallest q (7e-18) does nothing.
- Viterbi tracks 111-173 nodes (2-7x the constant rate), but no safe
  threshold saves time: every optimal strategy is within noise of exact,
  and prune_qmed 1 is suboptimal on half or more of the queries with no
  saving. Time savings (0.36-0.57) come only with wrong paths.
- relative_qmed 1 saves 16-24% of forward-backward values and time with
  mismatches, with logL error <= 1.3e-3 and TV <= 3e-5.

The README records the strategies dropped from further exploration:
floors, top_n, fixed relative fractions, escalation, and lossless_prune
at the smallest q.
The thresholding evaluations (experiments/thresholding/README.md) dropped
top_n: it turns off sparse sites and sorts every site, so it is 50-270x
slower than exact on sc2ts, it is inexact with mismatches, and it has no
guarantee. Remove it from the C library and the Python API.

With top_n gone, relative pruning is the only experimental mode, so the
(mode, parameter) pair becomes a single threshold:
tsls_hmm_set_relative_prune_threshold in C, and
Hmm(relative_prune_threshold=f) in place of pruning=("relative", f).
0 turns it off and keeps sparse sites; before, ("relative", 0) processed
every site in full. relative_q pruning (f = c q) is kept for
forward-backward.

The other dropped strategies (floors above q, fixed relative fractions,
escalation, lossless_prune at the smallest q) are parameter choices for
mechanisms that stay: likelihood_threshold as an underflow guard,
prune_threshold, and relative pruning.

evaluate.py, thresholding.py and verification.py use the new option, and
no longer run top_n. 448 Python tests and the C tests pass.
Hmm(log_value_step=d) rounds each normalised Viterbi value to the nearest
exp(-k d) for integer k, so that values within a factor of about exp(d)
become equal and can be coalesced. With a prune threshold t, values take
at most about log(1 / t) / d distinct levels. Paths are approximate:
each switch or mismatch can change a path's log odds by up to d / 2.

This is to test whether limiting the number of distinct likelihood
values recovers the compression lost under a recombination map, where
the switch odds differ at every site (Viterbi stored 8x the traceback of
the constant rate on chr20p samples). Forward and backward are unchanged.

Tests: step 0 changes nothing, a tiny step keeps the optimal likelihood
on all examples, a step of 2 reduces the mean number of distinct values
per site (6.9 to 3.7 on a test case with random per-site rates), and
negative or non-finite steps are rejected.
experiments/coarsening/coarsening.py runs Viterbi on 6 chr20p samples
(mu = 1e-20, 1e-8) under the HapMap II chr20 map and coarsened versions,
scoring every path under the map. Ratios to the exact map run:

  variant              time       traceback  loss (mean)
  constant 0.01        0.44-0.50  0.13       911-917
  constant map mean    0.46-0.49  0.15-0.21  416-623
  uniform rate per bp  0.94-1.05  0.95       150-236
  smooth 100kb / 1Mb   0.93-1.04  0.95       148-236
  round r to 2^k       0.90-0.98  0.82-1.00  14-21
  round r to 10^k      0.77-0.85  0.59-0.82  70-104
  value step 0.1       1.22-1.34  0.97       0.1-0.2
  value step 0.5       1.14-1.20  0.87-0.88  2.0-2.7
  value step ln 10     0.90-1.02  0.59-0.62  33-55

The map's cost comes from r varying between sites, not its scale: a
constant r at the map's mean is as cheap as tsinfer's 0.01, and even the
620 distinct r of a uniform rate per bp (from site spacing) cost as much
as the map. Spatial smoothing doesn't help. Rounding r to a factor of 2
loses about 0.1 per switch (within the map's precision) but saves little;
rounding values is more accurate than rounding r for the same traceback,
but the extra log and exp make it no faster. Neither recovers more than
about a quarter of the map's extra time.
This reverts commit 69af95b. Rounding Viterbi values to a log grid
didn't pay off on chr20p samples under the HapMap II map
(experiments/coarsening/README.md): it reduced the traceback by up to 40%
but the extra log and exp per value made it no faster than exact, and it
loses accuracy (log odds 33-55 at a step of ln 10).
The log_value_step engine option has been removed, so coarsening.py no
longer runs the value-rounding strategies. The README keeps their results
and says how to rerun them (commit 69af95b added the option, and b3d970f
has the strategies).
@jeromekelleher
jeromekelleher merged commit e9c3803 into tskit-dev:main Oct 28, 2026
14 of 15 checks passed
@jeromekelleher
jeromekelleher deleted the drop-top-n branch October 28, 2026 19:23
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant