Repository navigation
[auto-optimize]: Speed up modisco motifs with byte-identical output - #138
Conversation
…er-pattern work) onto A12
…apped-kmer kernel) onto A15
…cked-Spearman follow-up onto A25
…son_correlation overhead) onto A28
…d A37 (batched SeqletSet init) onto A35
…pile placement onto A43
…aptation) onto A45
git diff autoopt/1 autoopt/4 | git apply -3; clean in affinitymat.py and gapped_kmer.py. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
git diff autoopt/2 autoopt/6 | git apply -3; affinitymat.py clean, precompile.py conflicted in the kernel list. Resolved: _fwd_rev_best takes _jaccard_best's child, and _jaccard_best is queued behind B3's density kernels (B6 queued it behind binary_search_perplexity, which B3 folded into _perplexity_betas). Five groups, as before. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
git diff autoopt/8 autoopt/14 | git apply -3; clean in affinitymat.py, core.py and precompile.py. B15 (B12 + B9) does not touch these three beyond B12, so each is byte-identical to B18's (B14 on B12). Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…o B15 + B14 git diff autoopt/12 autoopt/16 | git apply -3; tfmodisco.py clean. The applied +/- lines are B16's own, and the file differs from B16's by exactly B9's lines (B15), both compared as line multisets. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…15 + B14 + B16 git diff autoopt/12 autoopt/17 | git apply -3; clean in cli.py, cluster.py, extract_seqlets.py and tfmodisco.py. Then cluster._leiden is put back as in B15: find_partition gets weights=weights, not weights.tolist(), and the docstring lines that explained the list are dropped. B17's README measured that hunk at ~200 MB more in each main Leiden clustering for 0.01-0.02 s. The rest of B11's cluster.py change is kept; cli.py and extract_seqlets.py are byte-identical to B17's. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…tted pair-map weights) onto B21
git diff autoopt/26 autoopt/27 | git apply -3; clean in affinitymat.py, aggregator.py, core.py, precompile.py and util.py. The tree is identical to autoopt/27's. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
git diff autoopt/21 autoopt/25 | git apply -3; clean in affinitymat.py and precompile.py. The applied +/- lines are B25's own (line multisets), and every name B25 added or changed in affinitymat.py is byte-identical to B25's; jaccard_from_seqlets and _fwd_rev_best are the same on B21 and B26. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…anges PR 136's tests, merged into master, call private functions B28 changed: - conftest.py's session warm-up and test_sparse_mm_dot (15 cases), test_sparse_mm_dot_max and test_sparse_mm_dot_ties call affinitymat._sparse_mm_dot with its old 7 arguments; it now takes 10 and the pipeline calls _sparse_knn. On all 17 of those inputs, _sparse_knn returns similarities and neighbours byte-identical to upstream's _sparse_mm_dot. - test__jaccard_raises expects _jaccard to reject float64 and int64 inputs; B28 dropped its eager numba signature, so they compile instead. The test checks no value. Removed under the user's rule: tests that fail only because of a change in signature or API are removed; a test that returns an incorrect value is reported instead (none did). The _sparse_mm_dot import, used only by these tests, goes with them. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
… one helper) onto B2
…nd B6 (Jaccard prefetch, aligned accumulators) onto B5 (cherry picked from commit 28c0d64ea2bd45b2b6260fadfccf37cdd4d3c816)
…polishing merge_in_seqlets_filledges holds the parent, once it has been expanded, as a "shell": a SeqletSet whose seqlets are coordinate tuples, with the counts and running sums the SeqletSet would have. Its expansions compute the new coordinates and keep them as SeqletSet would (first of each string), and add up only the new flanks when every seqlet was kept and the sums are the kernel's (the interior positions take the same additions as before). Its seqlets are made once, when it is returned. polish_pattern takes its two trims and its first expansion the same way and makes only the last expansion's seqlets. Anything these paths do not cover (a kernel that does not apply, a NaN, a window that does not fit) runs the code as it was. Also: TrackSet._coordinate_sums (the seqlet sums kernel from coordinates, which _add_seqlet_sums now calls), SeqletSet._support_window (trim_to_support's indices), Seqlet.string as an f-string of str() of each coordinate, and no super().__init__() call in Seqlet.__init__. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> (cherry picked from commit 444ef0d4ecfe6888619bbd0b8331bb9984397ba4)
…e kernel SeqletSet(seqlets, track_set=...) at the places that build patterns from seqlets the track set made: each subpattern (_set_subclusters, through compute_subpatterns' new optional track_set), the patterns a forked child built (_patterns_from_clusters' unpack_share) and the negative metacluster's patterns (_unpack_patterns). TrackSet._add_seqlet_sums then adds the sums up from the coordinates in seqlet order, as B23 did for the expansions. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> (cherry picked from commit f023fcee683bb42a92d8b3f1ad123ab35a822d49)
…f numpy's pairwise_sum
|
Summary: Three (potential) minor issues?
[Below: Codex review 🤖] I compared PR head 7e4377e against base 5b37ef7 on Linux/Python 3.10.18. All 1,450 runnable tests passed, with 28 existing bug skips, and the real-input cases I tested produced matching outputs. I found three issues below. The first two examples use the PR checkout with these dependency versions: python -m pip install .
The PR description already acknowledges this possibility. The example below demonstrates it: a normalized, symmetric graph with finite positive weights and diagonal entries produces four communities using upstream’s graph construction and five using the PR’s construction. Both calls use the same PR clustering code, dependency versions, and seeds; only the graph builder changes. Output: The pairwise membership comparison confirms changed communities beyond label renumbering. Preserving upstream’s ordered parallel-edge representation would avoid this compatibility change.
forked.budget() returns numba.config.NUMBA_NUM_THREADS, the configured maximum, even after the caller reduces the effective limit with numba.set_num_threads(). On Linux with the fork-safe workqueue backend, limiting Numba to one thread still starts three compilation processes: Output: This budget also controls Leiden workers and other process task plans. The effective budget should be captured at run entry, before temporary stage-specific reductions, and respected throughout the run.
The new regular dependency leidenalg==0.11.0 requires igraph>=1.0.0,<2.0, whose available releases require Python ≥3.9. The package still declares requires-python = ">=3.7". The following commands ask pip to resolve Python 3.8 wheels from a newer interpreter; --dry-run leaves the environment unchanged: The PR’s pinned requirement fails: The previous requirement resolves successfully: Would install igraph-0.11.8 leidenalg-0.10.2 texttable-1.7.1 If dropping Python 3.8 is intentional, please update requires-python and document the change. Otherwise, the dependency selection needs to remain compatible with Python 3.8. |
|
Suggested the following: The reason for pinning is sound, but it doesn’t explain issue 1. PR #136 reports different clustering with Leiden 0.11.0 and 0.12.0. Our counterexample uses Leiden 0.11.0 and igraph 1.0.0 for both calls in the same process, with identical inputs and seeds. Changing only graph construction changes four communities into five. PR #136 (#136) I suggest these fixes:
NUMBA_NUM_THREADS is the maximum; get_num_threads() reflects the runtime setting. Numba documentation (https://numba.readthedocs.io/en/stable/user/threading-layer.html#setting-the-number-of-threads) For scheduling that needs a stable total budget, capture that effective value at run entry before precompilation changes thread masks, then divide it among compiler children and execution workers. Restore
Document the compatibility change and retain CI coverage for supported versions. If Python 3.8 must remain supported, use interpreter-specific Leiden dependencies: Update the test extra too, and replace the new os.waitstatus_to_exitcode() calls—it was introduced in Python 3.9. For the current zero-exit checks, this compatible predicate suffices: Then test installation and the fork paths on Python 3.8. Python 3.9 release notes (https://docs.python.org/3.9/whatsnew/3.9.html#os) |
This optimization was found with auto-optimize, a tool developed by the Programmable Genomics Laboratory at UMass Chan Medical School: the
modisco motifspipeline is parallelized across processes, its numba kernels are compiled in background processes and rewritten, and its clustering, merging, polishing and start-up stages are reworked, cutting a 50,000-seqlet benchmark on real attribution scores from 1745 s to 119.1 s with output byte-identical to upstream's under leidenalg 0.11.0, whichpyproject.tomlnow pins.Changes, in order
pyproject.tomlpinsleidenalg==0.11.0. It also removes the tests that failed only because of those changes' signatures and APIs: theconftest.pywarm-up call to_sparse_mm_dot, the three_sparse_mm_dottests (the pipeline now calls_sparse_knn, which returns the same values on all their inputs) andtest__jaccard_raises.@njit(nogil=True)kernel over blocks of 16 rows on a Python thread pool instead of aprangekernel.modiscolite.util._pairwise_sum, a numba translation of numpy's pairwise summation; comments only.Translated code
modiscolite.util._pairwise_sumis a translation into numba ofpairwise_sumin numpy'snumpy/_core/src/umath/loops_utils.h.src(numpy 2.0.1), so that sums over rows of the density-adapted affinity matrix are added in numpy's order. It is under numpy's BSD-3-Clause license, whose notice is reproduced above the function. No other code is copied or translated from another project; routines that reproduce scikit-learn's ROC AUC and scipy's sparse operations are written from their documented behaviour.Where this was measured
This branch is the measured code, unchanged:
masterat 5b37ef7 plus the 39 commits of the search's path. The samechange is also proposed to jmschrei/tfmodisco-lite as a port onto its
main(jmschrei/tfmodisco-lite#78), where
modiscolite/andtests/are byte-identical to thisbranch.
Checks outside the search
Clean room. The measured branch merged into a fresh clone of kundajelab/tfmodisco and installed the ordinary way,
then timed against upstream installed the same way, three runs each on the 50,000-seqlet benchmark at 16 threads:
upstream 1726 s (1723, 1726, 1731), this branch 118.7 s (123.4, 118.5, 118.7). Outputs identical to upstream's on
the benchmark and on 20,000-seqlet, 5,000-seqlet and 480-sequence inputs in every run, and unit tests passing in both
environments. The verdict of the check is
differs: upstream ran 19 s (1.1%) faster than the search measured it(1745 s, on another machine with numpy 2.0.1 and numba 0.67.0), so the gain differs from the search's by 17.9 s,
more than the 5.8 s noise band; this branch's time matched the search's (119.1 s).
CI, run locally on this branch's head:
tests.ymlon Python 3.9, 3.10, 3.11, 3.12 and 3.13 (numpy 2.0.2 to 2.5.4,numba 0.60.0 to 0.68.0): 1,450 passed and 28 skipped on each;
publish.yml's build: the sdist and the wheel build andpass
twine check. The suite has 18 tests fewer thanmain's 1,468: the ones A1 removes.Every command and flag. Both versions, installed, ran 4,972 jobs: every subcommand (
motifs,report,report-simple,meme,seqlet-bed,seqlet-fasta,convert,convert-backward) with every flag one at a time,invalid values included;
motifsat default options on every attribution file (60 models; gradients, DeepLIFT/SHAP,integrated gradients, ISM) of a test-data collection's 10-, 1,000- and 10,000-locus tiers, and the twelve downstream
commands on each output; noise and constant attribution tracks; 40 malformed inputs (NaN, inf, wrong shapes, dtypes,
layouts, keys, truncated and empty files, both legacy h5 layouts); 49 malformed MEME files; and
-n 1000,10000and100000on 44 whole-experiment attribution sets of 3,502 to 242,924 peaks. h5 outputs were compared by content,every other output (HTML, logos, MEME, BED, FASTA), stdout and exit codes by bytes, and failures by their error line.
Every job matched, or failed the same way in both versions, except 19:
subpattern_1andsubpattern_2exchanged, with identical contents). Two Leiden seeds find the same partition, numbered differently;their modularities tie exactly in upstream's arithmetic, and this branch's Leiden graph, which adds each
undirected edge once with weight w + w instead of twice, breaks the tie by one unit in the last place.
-n 1fail in both versions with exit 1 and aValueError, from different lines: upstream's"need at least one array to concatenate" while building the Leiden graph, this branch's "max() iterable argument
is empty" one step later.
-g -1(a negative flank) fail in both versions with exit 1 and aValueError, at different stages:upstream's "negative dimensions are not allowed" in a later polish, this branch's broadcast error in its first,
coordinate-based polish.
In the same run (up to 20 jobs side by side at 8 threads each, so a single job's time includes contention), the
median time of
motifsat-n 100000across the 44 whole-experiment sets was 3,165 s for upstream and 197 s for thisbranch; on the four largest sets upstream took 3.2 to 3.4 hours and this branch 9 to 10 minutes.
What the command-line evaluation covered
Both versions were installed in their own environments and each job ran every command once per version, in the same
working directory, on js031 (2 x AMD EPYC 9575F). Counts are from the job results (
coverage_stats.py).Inputs. Attributions from 60 models and model heads in 17 families (BPNet, ChromBPNet, Cherimoya, ProCapNet, Sei,
Beluga, Malinois, DeepSTARR, LegNet, ENIGMA, Puffin, Puffin-D, CLIPNET, seq2PRINT, Decima, EPInformer, gReLU), by four
attribution methods: gradients (410
motifsjobs), DeepLIFT/SHAP (507), integrated gradients (372) and in silicosaturation mutagenesis (200 over the whole window, 331 over its central 50 to 1,000 bp). 410 attribution sets: 151 of
10 loci, 151 of 1,000, 64 of 10,000, and 44 whole experiments (3,502 to 242,924 peaks) each run at
-n1,000,10,000 and 100,000. 31 locus lists over 13 biosamples (cell lines, tissues, a fly embryo and yeast) and synthetic
sequence, on hg38, hg19, dm6 and sacCer3.
Controls.
motifsjobs, with their downstream commandsmotifsjobs.npy, both legacy h5 layouts, transposed, 2-D, mismatched counts, lengths and channels, wrong or extra npz keys, 0 and 1 loci, odd width, Fortran order, truncated, empty and text files, an h5 with neither layout)-s/-atooreport,report -landreport-simple--helpComparisons. 4,972 jobs, 19,234 command pairs (38,468 processes), 439,585 output files compared.
motifsmemereportseqlet-fastaseqlet-bedreport-simpleconvert-backwardconvertOf the 1,873
motifsruns on data that both versions completed, 1,121 output at least one pattern: 13,010 patternswith 3,480,920 seqlets, identical apart from the three subpattern numberings. The other 752 hold none (almost all from the
10-locus sets), so they compare empty results.
Not covered. 750 downstream jobs did not run, because their
motifsrun failed in both versions; 916motifspairs show only that both fail the same way. Flags were swept one at a time on four representative sets per size,
with four combinations. Everything ran on Linux with leidenalg 0.11.0, the version
pyproject.tomlpins; other numpyand numba versions were exercised only by the unit tests (Python 3.9 to 3.13 above). The Leiden tie behind the three
numbering differences comes from summation order, and on other inputs it could in principle pick a seed with a
different partition; that was not observed.
Usage benchmarks
Machine: AMD EPYC 9575F (2 x 64 cores, 256 hardware threads), Ubuntu 24.04.5 (Linux 7.0.0), Python 3.13.5, numpy 2.5.4,
numba 0.68.0, scipy 1.18.1, scikit-learn 1.9.1, leidenalg 0.11.0, igraph 1.0.0; base is
masterat 5b37ef7 andwinner this branch (B17), each installed with
uv pip install . pytest leidenalg==0.11.0from agit archiveexportinto its own environment. Synthetic
inputs: none were generated for this sweep; the inputs are real DeepLIFT/SHAP attributions from public ENCODE-trained
models in a test-data collection, the search's own Acinar data, and one negative control whose attributions are
Gaussian noise (sd 1, drawn with a seed fixed by the collection). Each row is one
modisco motifsrun per version,one configuration at a time on an otherwise idle machine, each a fresh process with an empty numba cache (start-up and
compilation included),
NUMBA_NUM_THREADSand the BLAS thread counts set to the row's threads. Peak is the wholeprocess tree's
memory.peakfrom its own cgroup (forked workers included, page cache included). Gate: the winner'sh5 equals upstream's in every group, dataset and attribute (dtype, shape and bytes).
-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 2000-n 20000 -w 200-n 20000 -w 200-n 20000 -w 200-n 20000 -w 200-n 20000 -w 200-n 20000-n 20000-n 20000-n 20000-n 20000-n 20000-n 20000-n 20000-n 20000-n 20000-n 10000-n 10000-n 10000-n 10000-n 10000The outputs passed the gate in every configuration.
Path from the baseline
Every value is the median over that iteration's runs (one seed per run); Δ is against the row above. Gates: h5_mismatches_vs_base == 0, h5_mismatches_vs_base_20k == 0, h5_mismatches_vs_base_5k == 0, h5_mismatches_vs_base_tiny == 0.
🤖 Generated with Claude Code