trviz Cython optimization PR — validation harness (1000-case differential test) and microsecond benchmark

bench_us.py

"""Microsecond-precision benchmark for trviz.cy.decompose.decompose_cy.

Times decompose_cy on 6 representative cases (perfect repeats of varying
motif length, multi-motif mixes, mutated long sequences). Reports
min / median / mean over N iterations per case (default N=20).

Usage:
  python3 bench_us.py [n_iter]

To compare two versions, run the script against each (e.g. on different
git branches with the Cython extension rebuilt in between) and compare
the median columns.

Dependencies:
  trviz (with the Cython extension built — provides
         trviz.cy.decompose.decompose_cy)
"""
import time, statistics, sys
from trviz.cy.decompose import decompose_cy

CASES = [
    ("CGG x 2000",      "CGG" * 2000,                                  ["CGG"]),
    ("CGG x 10000",     "CGG" * 10000,                                 ["CGG"]),
    ("AAAAAA x 5000",   "AAAAAA" * 5000,                               ["AAAAAA"]),
    ("multi-3 x 3000",  "ACGT" * 3000 + "ACTT" * 1000 + "ACCT" * 1000, ["ACGT", "ACTT", "ACCT"]),
    ("multi-3 long mut",("ACGT" * 100 + "ACTT" * 50 + "ACCT" * 50) * 5, ["ACGT", "ACTT", "ACCT"]),
    ("AGTTAT/GTCT mix", "AGTTAT" * 2000 + "GTCT" * 2000 + "AGTTAT" * 1000, ["AGTTAT", "GTCT"]),
]

if __name__ == "__main__":
    n_iter = int(sys.argv[1]) if len(sys.argv) > 1 else 20
    print(f"{'case':<25} {'min(us)':>12} {'med(us)':>12} {'mean(us)':>12}")
    for label, seq, motifs in CASES:
        times = []
        for _ in range(n_iter):
            t0 = time.perf_counter()
            decompose_cy(seq, motifs, {})
            times.append((time.perf_counter() - t0) * 1e6)
        print(f"{label:<25} {min(times):12.1f} {statistics.median(times):12.1f} {statistics.mean(times):12.1f}")

harness.py

"""Differential correctness test for trviz.cy.decompose.decompose_cy.

Generates 1000 deterministic test cases (seeded RNG) covering single perfect
repeats, mutated repeats, multi-motif sequences, score-parameter sweeps,
tie-prone tiny-alphabet inputs, and long sequences. Runs each through
decompose_cy and dumps the inputs + outputs as sorted JSON.

This script does NOT check correctness on its own — it is an output
recorder. To verify a code change is output-preserving, run it once
against the unmodified code and once against the modified code, then
diff the two JSON files:

  python3 harness.py baseline.json      # with original code installed/built
  # ... rebuild or reinstall the changed code ...
  python3 harness.py modified.json      # with modified code installed/built
  diff -q baseline.json modified.json   # exit status 0 == byte-identical

Dependencies:
  trviz (with the Cython extension built — provides
         trviz.cy.decompose.decompose_cy)
"""
import json, random, sys
from trviz.cy.decompose import decompose_cy

ALPHABET = "ACGT"


def _ordered_unique(items):
    seen = []
    for x in items:
        if x not in seen:
            seen.append(x)
    return seen


def _mutate(seq, rng, sub_rate=0.05, ins_rate=0.02, del_rate=0.02):
    """Apply substitutions / insertions / deletions independently per position."""
    out = []
    for ch in seq:
        if rng.random() < del_rate:
            continue  # delete
        if rng.random() < sub_rate:
            out.append(rng.choice(ALPHABET))
        else:
            out.append(ch)
        if rng.random() < ins_rate:
            out.append(rng.choice(ALPHABET))
    return "".join(out) or seq[:1]  # never return empty


def _gen_cases():
    """1000 diverse, deterministic test cases."""
    rng = random.Random(0xC0FFEE)
    out = []

    # ---------- Block A: hand-picked from the existing test suite (24) ----------
    suite = [
        ("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA", ["AAAAAA"], {}),
        ("ACTGACTGACTG", ["ACTG"], {}),
        ("AACCTTTTCTAACCTTTTCT", ["AACCTTTTCT"], {}),
        ("CGG" * 20, ["CGG"], {}),
        ("AAAAAC" "AAAAAA" "AAAAAT" "AAAAAA" "TTAAAA", ["AAAAAA"], {}),
        ("ACTG" "ACTT" "ACTG", ["ACTG"], {}),
        ("AACCTTTTCT" "AACCTTGTCT", ["AACCTTTTCT"], {}),
        ("CGCCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGT", ["CGG"], {}),
        ("AAATA" "AAATT" "AAATAA" "AAATA", ["AAATA"], {}),
        ("AAAAAC" "AAAAAA" "AAAAAT" "AAAAAA" "TTAAAA", ["AAAAAA", "TTAAAA"], {}),
        ("ACTG" "ACTT" "ACTG", ["ACTG", "ACTT"], {}),
        ("AACCTTTTCT" "AACCTTGTCT" "AACCTTGTCT", ["AACCTTTTCT", "AACCTTGTCT"], {}),
        ("CGCCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGGCGT", ["CGG", "CGC", "CGT"], {}),
        ("AAATA" "AAATT" "AAATA" "AAAATA", ["AAATA", "AAATAA"], {}),
        ("ACCCA" "ACCC" "ACCCA" "ACCCA", ["ACCCA"], {}),
        ("ACT" "ACT" "ACC" "ACT", ["ACT"], {}),
        ("ACT" "ACT" "ACCG" "ACT", ["ACT"], {}),
        ("ACT" "ACT" "AC" "CG" "ACT", ["ACT", "AC", "CG"], {}),
        ("AATAA" "AATAAA" "AATAA", ["AATAA"], {}),
        ("ACGTTTACGTTTACGTTTACGTTT", ["ACGTTT"],
         {"match_score": 5, "mismatch_score": -2, "insertion_score": -3, "deletion_score": -3}),
        ("ACGTTTACGTTTACGTTTACGTTT", ["ACGTTT"], {}),
        ("ACGTACGTTACGTAACGT", ["ACGT"],
         {"match_score": 2, "mismatch_score": -1, "insertion_score": -1, "deletion_score": -1}),
        ("ACGTACGTTACGTAACGT", ["ACGT"],
         {"match_score": 2, "mismatch_score": -1, "insertion_score": -2, "deletion_score": -2}),
        ("ACCCAACCCACCCAACCCA", ["ACCCA"],
         {"match_score": 2, "mismatch_score": -1, "insertion_score": -1, "deletion_score": -1}),
    ]
    out.extend(suite)

    # ---------- Block B: single perfect repeats, varied motif length & copies (200) ----------
    for _ in range(200):
        motif = "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 10)))
        out.append((motif * rng.randint(2, 25), [motif], {}))

    # ---------- Block C: single motif + low-rate mutations (200) ----------
    for _ in range(200):
        motif = "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 8)))
        seq = motif * rng.randint(3, 15)
        out.append((_mutate(seq, rng, 0.05, 0.02, 0.02), [motif], {}))

    # ---------- Block D: single motif + heavier mutations (100) ----------
    for _ in range(100):
        motif = "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 8)))
        seq = motif * rng.randint(3, 15)
        out.append((_mutate(seq, rng, 0.15, 0.08, 0.08), [motif], {}))

    # ---------- Block E: multi-motif, sequences sampled from motif set (200) ----------
    for _ in range(200):
        n_motifs = rng.randint(2, 5)
        motifs = _ordered_unique(
            "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 7)))
            for _ in range(n_motifs)
        )
        if len(motifs) < 1:
            continue
        n = rng.randint(3, 15)
        seq = "".join(rng.choice(motifs) for _ in range(n))
        out.append((seq, motifs, {}))

    # ---------- Block F: multi-motif + mutations (150) ----------
    for _ in range(150):
        motifs = _ordered_unique(
            "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 6)))
            for _ in range(rng.randint(2, 4))
        )
        n = rng.randint(3, 12)
        seq = "".join(rng.choice(motifs) for _ in range(n))
        out.append((_mutate(seq, rng, 0.08, 0.04, 0.04), motifs, {}))

    # ---------- Block G: score-parameter sweep (80) ----------
    for _ in range(80):
        motif = "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 6)))
        seq = motif * rng.randint(3, 8)
        out.append((seq, [motif], {
            "match_score": rng.choice([1, 2, 3, 5, 10]),
            "mismatch_score": rng.choice([-5, -3, -2, -1, 0]),
            "insertion_score": rng.choice([-5, -3, -2, -1]),
            "deletion_score": rng.choice([-5, -3, -2, -1]),
        }))

    # ---------- Block H: score sweep with mutations + multi-motif (50) ----------
    for _ in range(50):
        motifs = _ordered_unique(
            "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 5)))
            for _ in range(rng.randint(2, 3))
        )
        seq = "".join(rng.choice(motifs) for _ in range(rng.randint(3, 8)))
        seq = _mutate(seq, rng, 0.1, 0.05, 0.05)
        out.append((seq, motifs, {
            "match_score": rng.choice([1, 2, 5]),
            "mismatch_score": rng.choice([-3, -2, -1]),
            "insertion_score": rng.choice([-3, -2, -1]),
            "deletion_score": rng.choice([-3, -2, -1]),
        }))

    # ---------- Block I: tie-prone cases — uniform scores force many ties (60) ----------
    for _ in range(60):
        # Motif over a tiny alphabet to maximize collisions / ties
        small = rng.choice(["AC", "AT", "GC", "AG"])
        motif = "".join(rng.choice(small) for _ in range(rng.randint(2, 6)))
        seq = motif * rng.randint(3, 10)
        seq = _mutate(seq, rng, 0.1, 0.05, 0.05)
        kw = {"match_score": 1, "mismatch_score": -1,
              "insertion_score": -1, "deletion_score": -1}
        out.append((seq, [motif], kw))

    # ---------- Block J: long sequences (heavy j>1 branch) (30) ----------
    for _ in range(30):
        motif = "".join(rng.choice(ALPHABET) for _ in range(rng.randint(3, 8)))
        seq = motif * rng.randint(50, 250)
        out.append((seq, [motif], {}))

    # ---------- Block K: long sequences w/ mutations (30) ----------
    for _ in range(30):
        motif = "".join(rng.choice(ALPHABET) for _ in range(rng.randint(3, 8)))
        seq = motif * rng.randint(50, 200)
        out.append((_mutate(seq, rng, 0.05, 0.02, 0.02), [motif], {}))

    # ---------- Block L: multi-motif long sequences (20) ----------
    for _ in range(20):
        motifs = _ordered_unique(
            "".join(rng.choice(ALPHABET) for _ in range(rng.randint(3, 6)))
            for _ in range(rng.randint(2, 4))
        )
        seq = "".join(rng.choice(motifs) for _ in range(rng.randint(40, 120)))
        out.append((seq, motifs, {}))

    # ---------- Block M: single short motif over long sequence (20) ----------
    for _ in range(20):
        motif = rng.choice(ALPHABET) + rng.choice(ALPHABET)  # length-2 motif
        seq = motif * rng.randint(20, 100)
        out.append((seq, [motif], {}))

    # ---------- Block N: many motifs (5-8) (40) ----------
    for _ in range(40):
        motifs = _ordered_unique(
            "".join(rng.choice(ALPHABET) for _ in range(rng.randint(2, 6)))
            for _ in range(rng.randint(5, 8))
        )
        seq = "".join(rng.choice(motifs) for _ in range(rng.randint(5, 15)))
        out.append((seq, motifs, {}))

    return out[:1000]  # cap at exactly 1000


def main(out_path):
    cases = _gen_cases()
    results = []
    for seq, motifs, kwargs in cases:
        try:
            r = decompose_cy(seq, motifs, dict(kwargs))
            results.append({"sequence": seq, "motifs": motifs, "kwargs": kwargs, "result": r})
        except Exception as e:
            # Capture the exception type+message; some inputs may exceed score threshold.
            results.append({"sequence": seq, "motifs": motifs, "kwargs": kwargs,
                            "error": f"{type(e).__name__}: {e}"})
    with open(out_path, "w") as f:
        json.dump(results, f, indent=2, sort_keys=True)
    n_err = sum(1 for r in results if "error" in r)
    print(f"Wrote {len(results)} cases to {out_path} ({n_err} errors)")


if __name__ == "__main__":
    main(sys.argv[1] if len(sys.argv) > 1 else ".opt_harness/baseline.json")