Index perf

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.

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):
 

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

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__"]

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()

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]
 

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()
 
 #############################################################################

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__