diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 47a1b1f..0a2cf89 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -143,6 +143,31 @@ for user-facing installation/usage docs. `benchmark:` reports as a sum across the process tree). So gz costs more time *and* more peak memory than the parallel-eligible compressions for the same `fetch` workload — worth knowing before running it on a memory-constrained node. +- **`mafutils index`'s hot loop: `TextIOWrapper.tell()` was the bottleneck, and + the fix is binary mode + a byte counter, NOT a chunked scanner.** Profiling a + real 1.5GB MAF slice showed `tell()` at ~60% of runtime — it was called once + per line (`iterMafBlocks` needs a position before every line) and cost ~4.6x + more than the `readline()` it accompanied, because on a `TextIOWrapper` + `tell()` must snapshot the incremental UTF-8 decoder state to build a + seekable cookie (visible as ~1.1M `codecs.getstate`/`setstate` calls). Since + `none`/`gz` have plain additive stream positions, `iterMafBlocks(..., + binary=True)` now reads a binary stream and tracks position with `pos += + len(line)`. Measured on 1.5GB of real data, scan-only: **11.89s -> 1.38s + (8.6x)**; end-to-end `mafutils index`: **23.82s -> 11.63s (2.05x)**, output + byte-identical. + Two approaches that sound faster but measured **slower**, so don't retry them + without new evidence: a chunked `bytes.find(b"\n")` scanner (2.17s — loses to + C-implemented `BufferedReader.readline()` while adding Python-level buffer + slicing), and a block-level scanner that finds `b"\na"` boundaries to avoid + per-line iteration entirely (2.71s — still has to walk every line to count + non-blank lines for `num_seqs`, plus buffer-concatenation cost). Remaining + headroom is real but modest: end-to-end sits at ~129 MB/s against a ~640 MB/s + ceiling (raw read + MD5 with no parsing at all), and closing it would require + assuming things about intra-block blank/comment lines that risk silently + corrupting an index. + **bgzip deliberately stays on the text+`tell()` path** — its virtual offsets + are not additive, so a byte counter cannot reproduce them (same root cause as + the straddling-`a`-line bug above). - The vendored `bgzf.BgzfReader` has no `.name` attribute (unlike `gzip.GzipFile`), so `fetch.py` gets the MAF's display filename from the known `maf_file` path (`WORKER_MAF_FILE`), not from the open file handle. diff --git a/README.md b/README.md index 6e925c5..3d9b0e1 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,12 @@ Show top-level help: mafutils --help ``` +Show the installed version: + +```bash +mafutils --version # also accepts -version, -v, -V +``` + Create block and scaffold indexes for a MAF (defaults to `input.maf.block.idx` / `input.maf.scaffold.idx` if output paths are omitted): diff --git a/benchmarks/baseline-manual-timings.txt b/benchmarks/baseline-manual-timings.txt index f02ac0b..49afae6 100644 --- a/benchmarks/baseline-manual-timings.txt +++ b/benchmarks/baseline-manual-timings.txt @@ -3,6 +3,12 @@ real 197.77 user 185.05 sys 11.46 +time -p mafutils index cricetid-15spec.Mus_musculus.nodupes.maf +real 122.69 +user 110.87 +sys 10.13 + + time -p mafutils validate cricetid-15spec.Mus_musculus.nodupes.maf [ 2026-07-17 12:36:33 - INFO ] mafutils validate called as: mafutils validate cricetid-15spec.Mus_musculus.nodupes.maf [ 2026-07-17 12:36:33 - INFO ] Using index at default location: cricetid-15spec.Mus_musculus.nodupes.maf.block.idx @@ -125,6 +131,11 @@ real 685.60 user 680.49 sys 2.38 +time -p mafutils index cricetid-15spec.Mus_musculus.nodupes.maf.gz +real 182.47 +user 178.17 +sys 2.30 + time -p mafutils validate cricetid-15spec.Mus_musculus.nodupes.maf.gz [ 2026-07-17 15:48:42 - INFO ] mafutils validate called as: mafutils validate cricetid-15spec.Mus_musculus.nodupes.maf.gz [ 2026-07-17 15:48:42 - INFO ] Using index at default location: cricetid-15spec.Mus_musculus.nodupes.maf.gz.block.idx @@ -233,6 +244,13 @@ real 255.30 user 251.65 sys 2.53 +time -p mafutils index cricetid-15spec.Mus_musculus.nodupes.maf.gz + +real 224.75 +user 219.95 +sys 2.93 + + time -p mafutils validate cricetid-15spec.Mus_musculus.nodupes.maf.gz [ 2026-07-17 15:53:14 - INFO ] mafutils validate called as: mafutils validate cricetid-15spec.Mus_musculus.nodupes.maf.gz [ 2026-07-17 15:53:14 - INFO ] Using index at default location: cricetid-15spec.Mus_musculus.nodupes.maf.gz.block.idx diff --git a/mafutils/__init__.py b/mafutils/__init__.py index eabd9b5..797ad3c 100644 --- a/mafutils/__init__.py +++ b/mafutils/__init__.py @@ -1,3 +1,13 @@ """Reusable MAF toolkit package for indexing, fetching, and summarizing MAFs.""" -__all__ = ["cli", "fetch", "index", "stats"] +from importlib.metadata import PackageNotFoundError, version as _installed_version + +try: + __version__ = _installed_version("mafutils") +except PackageNotFoundError: + # Running from a source tree that was never installed (no distribution + # metadata to read). setuptools-scm derives the real version from git + # tags at build/install time, so there's nothing to fall back to here. + __version__ = "unknown" + +__all__ = ["cli", "fetch", "index", "stats", "__version__"] diff --git a/mafutils/cli.py b/mafutils/cli.py index 6b99a21..6a9329e 100644 --- a/mafutils/cli.py +++ b/mafutils/cli.py @@ -1,5 +1,8 @@ +from typing import Annotated + import typer +from . import __version__ from . import fetch as fetch_mod from . import gc as gc_mod from . import index as index_mod @@ -20,5 +23,31 @@ app.command("validate", help="Check whether a MAF file's index is still trustworthy.")(validate_mod.validate_command) +def versionCallback(value: bool) -> None: + if value: + typer.echo(f"mafutils {__version__}") + raise typer.Exit() + + +@app.callback() +def mafutils_callback( + version: Annotated[ + bool, + typer.Option( + "--version", + "-version", + "-v", + "-V", + callback=versionCallback, + is_eager=True, + help="Show the installed mafutils version and exit.", + ), + ] = False, +) -> None: + # is_eager so --version is handled before any subcommand parsing, letting + # `mafutils --version` work on its own without a subcommand. + pass + + def main() -> None: app() diff --git a/mafutils/index.py b/mafutils/index.py index e6aaa82..871064e 100644 --- a/mafutils/index.py +++ b/mafutils/index.py @@ -19,11 +19,10 @@ def process_maf_block(block): ref_seq = block[1].split() ref_scaff = ref_seq[1].split(".", 1)[1] - ref_len = block[1].split()[3] line_len = str(len(block[1])) num_seqs = str(len(block) - 1) seq_len = str(len(ref_seq[6])) - return [ref_scaff, ref_seq[2], ref_len, seq_len, line_len, num_seqs] + return [ref_scaff, ref_seq[2], ref_seq[3], seq_len, line_len, num_seqs] def run_index(maf_file, block_index_path, scaffold_index_path): @@ -48,8 +47,12 @@ def run_index(maf_file, block_index_path, scaffold_index_path): region_start_byte = None region_end_byte = None - for block_text, block_start, block_end in COMMON.iterMafBlocks(maf_stream): - block = block_text.split("\n") + # none/gz get the fast binary+byte-counting path; bgzip must use + # text mode + tell(), since its virtual offsets aren't additive + # (see iterMafBlocks/openMafHashing). + use_binary = maf_compression in ("none", "gz") + + for block, block_start, block_end in COMMON.iterMafBlocks(maf_stream, binary=use_binary): block_info = process_maf_block(block) ref_scaffold = block_info[0] diff --git a/mafutils/lib/common.py b/mafutils/lib/common.py index 7c88e24..5742594 100644 --- a/mafutils/lib/common.py +++ b/mafutils/lib/common.py @@ -148,16 +148,31 @@ def openMafHashing(filename, compression, hash_obj): """ Like openMaf(), but threads a _HashingRawIO underneath the reader so every byte read from disk is fed into hash_obj -- for gz/bgzip this - means hashing the raw compressed bytes, before decompression. Returns a - stream with the same .readline()/.tell() interface iterMafBlocks() - expects. Used only by mafutils index. + means hashing the raw compressed bytes, before decompression. Used only + by mafutils index. + + For "none"/"gz" this returns a BINARY stream, to be consumed by + iterMafBlocks(..., binary=True): those two compression types have plain + additive stream positions, so byte offsets can be tracked with a simple + counter instead of calling tell(). That matters enormously at scale -- + on a TextIOWrapper, tell() must snapshot the incremental UTF-8 decoder's + state to build a seekable cookie, making it ~4.6x more expensive than + the readline() it accompanies, and iterMafBlocks needs a position before + every single line. Measured on a 1.5GB real MAF: 11.89s (text+tell) vs + 1.38s (binary+counting) for the same, byte-identical offsets. + + For "bgzip" this still returns a TEXT stream, because BGZF virtual + offsets pack a compressed-block offset and an in-block offset into one + int -- they are NOT additive, so a byte counter cannot reproduce them + and tell() is genuinely required. See iterMafBlocks for the full + explanation. """ hashing_raw = _HashingRawIO(open(filename, "rb"), hash_obj) if compression == "none": - return io.TextIOWrapper(hashing_raw) + return io.BufferedReader(hashing_raw) elif compression == "gz": - return io.TextIOWrapper(gzip.GzipFile(fileobj=hashing_raw, mode="rb")) + return io.BufferedReader(gzip.GzipFile(fileobj=hashing_raw, mode="rb")) elif compression == "bgzip": return bgzf.BgzfReader(fileobj=hashing_raw, mode="rt") else: @@ -391,52 +406,82 @@ def validateIndexHeader(header, maf_file, detected_compression, LOG, strict=Fals ############################################################################# -def iterMafBlocks(stream): +def iterMafBlocks(stream, binary=False): """ - Splits an open MAF stream into raw per-block text, split on 'a'-prefixed + Splits an open MAF stream into raw per-block lines, split on 'a'-prefixed header lines (skipping blank/comment lines). Yields - (block_text, block_start_offset, block_end_offset) tuples, where the - offsets are raw stream.tell() positions bracketing each block, suitable - for later random-access seeking back into the same file. - - Never computes a boundary by arithmetic on a tell() value (e.g. - tell() - len(line)) -- only ever uses tell() values exactly as the - stream reports them, captured at the right moment (before reading the - line that starts a new block). This matters for BGZF: its virtual - offsets pack a compressed-block-offset and an in-block-offset into one - int, and subtracting a plain byte count from that packed value can - "borrow" across the packing whenever an 'a' line straddles a real BGZF - block boundary, landing on a byte position that was never a valid block - start (confirmed against a real ~114k-block file, where this corrupted - the offset by exactly the size of that borrow). Capturing tell() - upfront sidesteps the packed representation entirely, and is exactly as - correct for the simpler none/gz cases too. + (block_lines, block_start_offset, block_end_offset) tuples, where the + offsets bracket each block and are suitable for later random-access + seeking back into the same file. + + block_lines[0] (the 'a' header) and block_lines[1] (the reference/first + sequence line) are always returned as stripped `str` -- they are the + only two lines whose *content* mafutils index reads. The remaining + sequence lines are returned exactly as readline() produced them, and in + binary mode are still `bytes`: only their *count* (len(block_lines)) is + ever used downstream, so decoding/stripping them would be pure overhead. + That matters enormously at scale -- a block can carry hundreds of + sequence lines (one per species) and a whole-genome MAF hundreds of + millions of blocks. + + Two position-tracking modes, because the two cases are fundamentally + different: + + binary=True (for "none"/"gz", via openMafHashing): the stream is binary + and positions are tracked with a plain additive byte counter. Valid + because both have simple additive stream positions (raw file bytes, and + decompressed-stream bytes, respectively). This exists for speed: on a + TextIOWrapper, tell() must snapshot the incremental UTF-8 decoder state + to build a seekable cookie, costing ~4.6x more than the readline() it + accompanies -- and a position is needed before every single line. + Measured on a 1.5GB real MAF: 11.89s (text+tell) vs 1.38s + (binary+counting), producing byte-identical offsets over 269,308 blocks. + + binary=False (for "bgzip"): positions come from stream.tell() and are + NEVER derived by arithmetic (e.g. tell() - len(line)). BGZF virtual + offsets pack a compressed-block offset and an in-block offset into one + int, so they are not additive -- subtracting or adding a plain byte + count can "borrow" across the packing whenever an 'a' line straddles a + real BGZF block boundary, landing on a byte position that was never a + valid block start (confirmed against a real ~114k-block file, where it + corrupted the offset by exactly the size of that borrow). This is also + why the binary counting path above must not be used for bgzip. """ + if binary: + empty, hash_prefix, a_prefix = b"", b"#", b"a" + else: + empty, hash_prefix, a_prefix = "", "#", "a" + block_lines = [] block_start = None - pos = stream.tell() + pos = 0 if binary else stream.tell() + next_line_is_reference = False while True: line = stream.readline() - if line == "": + if line == empty: break - next_pos = stream.tell() + next_pos = pos + len(line) if binary else stream.tell() - if line.startswith("#") or line.strip() == "": + if line.startswith(hash_prefix) or not line.strip(): pos = next_pos continue - if line.startswith("a"): + if line.startswith(a_prefix): if block_lines: - yield "\n".join(block_lines), block_start, pos + yield block_lines, block_start, pos block_start = pos - block_lines = [line.strip()] + block_lines = [line.decode().strip() if binary else line.strip()] + next_line_is_reference = True + elif next_line_is_reference: + block_lines.append(line.decode().strip() if binary else line.strip()) + next_line_is_reference = False else: - block_lines.append(line.strip()) + block_lines.append(line) pos = next_pos if block_lines: - yield "\n".join(block_lines), block_start, stream.tell() + yield block_lines, block_start, pos if binary else stream.tell() ############################################################################# diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 0000000..3e63b30 --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,65 @@ +import os +import re +import subprocess +import sys + +import pytest + +TEST_DIR = os.path.dirname(__file__) +REPO_ROOT = os.path.abspath(os.path.join(TEST_DIR, "..")) + + +def run(args): + return subprocess.run( + [sys.executable, "-m", "mafutils"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True, + ) + + +@pytest.mark.parametrize("flag", ["--version", "-version", "-v", "-V"]) +def test_version_flag_prints_version_and_exits_cleanly(flag): + result = run([flag]) + assert result.returncode == 0, result.stderr + # e.g. "mafutils 0.3.0" / "mafutils 0.3.1.dev3+g1234abc" / "mafutils unknown" + assert re.match(r"^mafutils \S+$", result.stdout.strip()), result.stdout + + +def test_version_flag_works_without_a_subcommand(): + """ + --version is is_eager so it short-circuits before subcommand parsing -- + without that, `mafutils --version` would fail asking for a COMMAND + (the app is built with no_args_is_help=True). + """ + result = run(["--version"]) + assert result.returncode == 0 + assert "Usage:" not in result.stdout + + +def test_single_dash_version_is_not_parsed_as_grouped_short_flags(): + """ + `-version` is declared alongside the short flag `-v`, so a parser could + plausibly read it as grouped short options (-v -e -r -s -i -o -n) and + either error or silently do the wrong thing. Assert it's treated as the + version alias it's declared as. + """ + result = run(["-version"]) + assert result.returncode == 0, result.stderr + assert result.stdout.strip().startswith("mafutils ") + assert "no such option" not in (result.stdout + result.stderr).lower() + + +def test_top_level_callback_does_not_break_subcommands(): + """ + Guard against the @app.callback() added for --version interfering with + normal subcommand dispatch. + """ + result = run(["index", "--help"]) + assert result.returncode == 0, result.stderr + assert "block index" in result.stdout.lower() + + +def test_package_exposes_dunder_version(): + import mafutils + + assert isinstance(mafutils.__version__, str) + assert mafutils.__version__