tests/testthat/test-perf-regression.R

In misha: Toolkit for Analysis of Genomic Data

# Performance regression tests.
#
# These guard against the kind of subtle slowdowns that a C++ optimization
# audit accidentally introduced — they don't need to be precise benchmarks,
# only catch order-of-magnitude regressions in setup/per-track overhead.
#
# Helper functions and the budget rule live in helper-perf.R. Each test
# embeds a single `baseline_ms` measured on the lab box. Budget for a test is
# `baseline_ms * 3`; the 3x factor keeps noisy NFS hardware from false-firing
# while still catching genuine >3x slowdowns.
#
# Pick `baseline_ms` as the FASTEST wall-clock you can reproduce for the
# operation across versions you care about. If a future fix makes an
# operation faster, lower the number (do not raise for convenience — the
# whole point is to keep the bar tight). If a planned change makes it
# legitimately slower, justify it before relaxing the baseline.
#
# OPT-IN: skipped by default. They allocate tracks, do real I/O, and depend
# on wall-clock timing — running them inside a parallel `devtools::test()`
# under load would make timings noisy and slow the whole suite.
#
#     MISHA_PERF_TESTS=true R -e "devtools::test(filter='perf-regression')"

create_isolated_test_db()

# Shared fixtures (created once, reused across tests via gtrack.* lookups).
.perf_setup <- function(n_dense_tracks = 30, dense_bin_size = 50L) {
    chrom_len <- gintervals.all()[gintervals.all()$chrom == "chr1", "end"]
    full_iv <- gintervals(1, 0, chrom_len)

    dense_tracks <- vapply(seq_len(n_dense_tracks), function(i) {
        nm <- random_track_name(prefix = "test")
        gtrack.create_dense(nm, "perf fixture", full_iv,
            values = runif(1), binsize = dense_bin_size, defval = 0
        )
        nm
    }, character(1))

    sparse_track <- random_track_name(prefix = "test")
    sparse_intervs <- gintervals(1, seq(0, 1e6, by = 100), seq(50, 1e6 + 50, by = 100))
    sparse_intervs <- sparse_intervs[sparse_intervs$end <= chrom_len, ]
    gtrack.create_sparse(sparse_track, "perf fixture", sparse_intervs,
        values = runif(nrow(sparse_intervs))
    )

    list(
        dense_tracks = dense_tracks,
        sparse_track = sparse_track,
        chrom1_len = chrom_len,
        full_chr1 = full_iv,
        small_chr1 = gintervals(1, 1e6, 1e6 + 1000)
    )
}

.perf_cleanup <- function(fix) {
    for (t in fix$dense_tracks) try(gtrack.rm(t, force = TRUE), silent = TRUE)
    try(gtrack.rm(fix$sparse_track, force = TRUE), silent = TRUE)
}

# --- Tests ---

test_that("gextract setup over many dense tracks stays fast (MAP_POPULATE regression)", {
    # Bug history: v5.6.11 added MmapFile with MAP_POPULATE which forced
    # eager paging-in of every per-chrom file at every mmap. v5.6.17 removed
    # MAP_POPULATE and switched the per-chrom-per-track validation loops to a
    # metadata-only path. This test would have caught both.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 30)
    on.exit(.perf_cleanup(fix), add = TRUE)

    measured <- time_op(gextract(fix$dense_tracks, fix$small_chr1, iterator = 50))
    expect_perf_baseline(measured, "many-dense setup", baseline_ms = 112)
})

test_that("gextract single dense track full chr1 stays fast (bin-scan path)", {
    # Guards the inner read loop in GenomeTrackFixedBin::read_interval — the
    # hot path that the perf audit's mmap zero-copy was supposed to speed up.
    # On chr1 (~250Mbp) at binsize=50 this scans ~5M bins.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)

    measured <- time_op(gextract(fix$dense_tracks[1], fix$full_chr1, iterator = 50))
    expect_perf_baseline(measured, "dense full-chr scan", baseline_ms = 2005)
})

test_that("gextract sparse track full chr1 stays fast", {
    # Guards GenomeTrackSparse::read_interval — uses a separate code path
    # (BufferedFile + load-into-memory) that didn't benefit from MmapFile but
    # could regress independently if the perf audit touched its hot loop.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 0)
    on.exit(.perf_cleanup(fix), add = TRUE)

    measured <- time_op(gextract(fix$sparse_track, fix$full_chr1))
    expect_perf_baseline(measured, "sparse full-chr", baseline_ms = 20)
})

test_that("gextract vtrack with avg + window stays fast", {
    # Realistic vtrack pattern: avg over a +/-N window. Touches the
    # iterator-modifier code path that the perf audit's hash-key changes
    # (BackendKey struct) sit on top of.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)
    on.exit(try(gvtrack.rm("v_avg"), silent = TRUE), add = TRUE)

    gvtrack.create("v_avg", fix$dense_tracks[1], func = "avg")
    gvtrack.iterator("v_avg", sshift = -250, eshift = 250)

    measured <- time_op(gextract("v_avg", fix$small_chr1, iterator = 50))
    expect_perf_baseline(measured, "vtrack avg+window", baseline_ms = 18)
})

test_that("gextract vtrack with LSE sliding window stays fast", {
    # LSE is the function Tamar uses for motif scoring. Worth its own test
    # because the LSE state machine in GenomeTrackFixedBin maintains a sliding
    # log-sum-exp across bins and is touched by several perf-audit changes.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)
    on.exit(try(gvtrack.rm("v_lse"), silent = TRUE), add = TRUE)

    gvtrack.create("v_lse", fix$dense_tracks[1], func = "lse")
    gvtrack.iterator("v_lse", sshift = -100, eshift = 100)

    measured <- time_op(gextract("v_lse", fix$small_chr1, iterator = 50))
    expect_perf_baseline(measured, "vtrack lse+window", baseline_ms = 18)
})

test_that("gquantiles vtrack LSE windowed full-chr scan stays fast (mode-3 sliding-bypass regression)", {
    # Bug history: v5.6.7 (commit 1cbfa801, "C++ optimization audit") added a
    # single-function fast path (mode 3) in GenomeTrackFixedBin that intercepted
    # single-function LSE *before* the mode-1 reducer path. Mode 3 recomputes the
    # log-sum-exp from scratch over the whole window on every output bin, bypassing
    # the incremental sliding-window LSE that mode 1 maintains. For a windowed lse
    # vtrack scanned bin-by-bin genome-wide - Tamar's motif-energy quantile
    # workload (compute_genomewide_motif_quantiles) - this was ~7.7x slower at a
    # 40-bin window (5111ms vs 667ms here) and grows with the window width.
    #
    # Two reasons the small_chr1 LSE test above could not catch it: (1) a ~20-bin
    # scope is dominated by fixed overhead, hiding the per-bin sliding advantage;
    # (2) it uses gextract, which materializes every output row back to R and so
    # dwarfs the C++ inner-loop cost. This test uses gquantiles (a streaming
    # reduction, exactly her workload) over full chr1 (~5M bins) with a wide
    # 40-bin window, so the per-bin LSE cost dominates and the regressed path
    # lands well outside the 3x budget.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)
    on.exit(try(gvtrack.rm("v_lse_full"), silent = TRUE), add = TRUE)

    gvtrack.create("v_lse_full", fix$dense_tracks[1], func = "lse")
    gvtrack.iterator("v_lse_full", sshift = -1000, eshift = 1000)

    measured <- time_op(gquantiles("v_lse_full", percentiles = 0.99, intervals = fix$full_chr1, iterator = 50))
    expect_perf_baseline(measured, "vtrack lse+window full-chr", baseline_ms = 667)
})

test_that("gscreen on many dense tracks stays fast", {
    # gscreen is the second hottest entry point after gextract for Tamar's
    # workloads. Same setup-cost characteristic as gextract — the validation
    # loops go through the same per-track path.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 30)
    on.exit(.perf_cleanup(fix), add = TRUE)

    expr <- paste(fix$dense_tracks, collapse = " > 0 & ")
    expr <- paste(expr, "> 0")

    measured <- time_op(gscreen(expr, intervals = fix$small_chr1, iterator = 50))
    expect_perf_baseline(measured, "gscreen many-dense", baseline_ms = 109)
})

test_that("gsummary on dense full chr1 stays fast", {
    # Reduction path — single-function fast-path territory in
    # GenomeTrackFixedBin. The audit added a "mode 3" single-function fast
    # path; this test would catch a regression there.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)

    measured <- time_op(gsummary(fix$dense_tracks[1], intervals = fix$full_chr1))
    expect_perf_baseline(measured, "gsummary dense full-chr", baseline_ms = 268)
})

test_that("gquantiles on dense full chr1 stays fast", {
    # Streaming-percentile path (StreamPercentiler) — the audit templated its
    # comparator for inlining. Different code path than the reducer/agg paths.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)

    measured <- time_op(gquantiles(fix$dense_tracks[1], percentiles = c(0.1, 0.5, 0.9), intervals = fix$full_chr1))
    expect_perf_baseline(measured, "gquantiles dense full-chr", baseline_ms = 429)
})

test_that("2D gextract on rect track full chr1xchr1 stays fast", {
    # 2D path is independent of the 1D hot paths above. The perf audit
    # touched StatQuadTree (uint8_t instead of bool) — that would show up here.
    skip_unless_perf()
    rect_tracks <- gtrack.ls("test\\.rects$")
    if (length(rect_tracks) == 0) skip("no rects track in test_db")

    iv2d <- gintervals.2d("chr1", 0, 1e7, "chr1", 0, 1e7)

    measured <- time_op(gextract(rect_tracks[1], intervals = iv2d, iterator = c(1e5, 1e5)))
    expect_perf_baseline(measured, "2D rect extract", baseline_ms = 21)
})

# ---- gmultitasking.strategy regression ---------------------------------
#
# The auto-strategy heuristic must NOT slow down small / few-track queries.
# These tests pin baselines for queries where 'auto' should pick 'tiles'
# and where it should pick 'tracks' — both must stay fast.

test_that("single-track gextract with auto strategy stays fast (must pick tiles)", {
    # 1-track gextract has nothing to track-split; auto must keep using tiles.
    # If a future change accidentally routes this through mclapply the wall
    # would jump because of fork+merge overhead.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 1)
    on.exit(.perf_cleanup(fix), add = TRUE)

    old <- options(gmultitasking.strategy = "auto")
    on.exit(options(old), add = TRUE)
    measured <- time_op(gextract(fix$dense_tracks[1], fix$full_chr1, iterator = 50))
    expect_perf_baseline(measured, "auto-strategy single-track full-chr", baseline_ms = 2005)
})

test_that("few-track small-iterator gextract with auto stays on tiles (no track-parallel overhead)", {
    # 5 tracks × 1000 intervals = 5000 < 1e6 threshold ⇒ auto picks tiles.
    # Re-uses the small_chr1 setup-overhead baseline.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 5)
    on.exit(.perf_cleanup(fix), add = TRUE)

    old <- options(gmultitasking.strategy = "auto")
    on.exit(options(old), add = TRUE)
    measured <- time_op(gextract(fix$dense_tracks, fix$small_chr1, iterator = 50))
    # Generous baseline: 5 tracks should be a fraction of the 30-track 112ms cost.
    expect_perf_baseline(measured, "auto-strategy few-track small-iter setup", baseline_ms = 60)
})

test_that("many-track large-iterator gextract: track-parallel >= tile-parallel", {
    # The flagship case: many tracks × many intervals on a warm-cache test_db.
    # Even on warm cache the track-parallel path should not be SLOWER than
    # tile-parallel — if it is, the heuristic or merge step has regressed.
    skip_unless_perf()
    fix <- .perf_setup(n_dense_tracks = 30)
    on.exit(.perf_cleanup(fix), add = TRUE)

    chrom_len <- fix$chrom1_len
    # Build an iterator scope that crosses the 1e6 threshold so 'auto' picks
    # tracks: 30 tracks × 50,000 intervals = 1.5e6 ≥ 1e6.
    starts <- seq.int(0, chrom_len - 100L, length.out = 50000L)
    starts <- as.integer(starts) - (as.integer(starts) %% 50L)
    big_scope <- gintervals(1, starts, starts + 50L)
    big_scope <- big_scope[!duplicated(big_scope$start), ]

    old <- options(gmultitasking.strategy = "tiles")
    on.exit(options(old), add = TRUE)
    t_tiles <- time_op(gextract(fix$dense_tracks, intervals = big_scope, iterator = big_scope))

    options(gmultitasking.strategy = "tracks")
    t_tracks <- time_op(gextract(fix$dense_tracks, intervals = big_scope, iterator = big_scope))

    # Track-parallel must not be more than 2× slower on warm cache. (It's
    # expected to win by ~5× cold; this guards against catastrophic merge or
    # mclapply-overhead regressions on the warm path.)
    msg <- sprintf(
        "tiles: %.1fms, tracks: %.1fms (ratio %.2fx)",
        t_tiles * 1000, t_tracks * 1000, t_tracks / t_tiles
    )
    if (t_tracks > t_tiles * 2) {
        fail(msg)
    } else {
        succeed(msg)
    }
})

test_that("an overlapping 1D iterator does not cost more than the same intervals disjoint", {
    skip_unless_perf()
    withr::local_options(gmultitasking = FALSE, gmax.data.size = 1e9)

    # Passing overlapping intervals as a 1D iterator runs a merge test over the whole
    # scope. The test used to walk back over every interval of a block for every scope
    # interval that touched it - fine while the intervals of a block sit near each other
    # (sliding windows), quadratic once one long interval spans the block, which is what a
    # domain with peaks nested inside it is. Measured on the lab box before the fix: 0.76s
    # against 0.09s for the same peaks without the container, and 8.7s cold.
    #
    # Stated as a ratio against the disjoint control rather than an absolute budget: both
    # sides extract the same peaks, so everything except the merge test cancels.
    k <- 40000
    st <- seq(0, by = 200, length.out = k)
    peaks <- gintervals(1, st, st + 100)
    nested <- gintervals(1, c(0, st), c(k * 200, st + 100))

    t_disjoint <- time_op(gextract("test.fixedbin", peaks, iterator = peaks))
    t_nested <- time_op(suppressWarnings(gextract("test.fixedbin", nested, iterator = nested)))

    msg <- sprintf(
        "disjoint: %.0fms, one container over %d nested peaks: %.0fms (ratio %.2fx)",
        t_disjoint * 1000, k, t_nested * 1000, t_nested / t_disjoint
    )
    if (t_nested > t_disjoint * 3) {
        fail(msg)
    } else {
        succeed(msg)
    }
})

test_that("an overlapping iterator that reaches back past the walk cap costs no step", {
    skip_unless_perf()
    withr::local_options(gmultitasking = FALSE, gmax.data.size = 1e9)

    # The merge test walks back over the intervals that reach into each scope interval, and
    # hands the query to a tree once the walk would run past MAX_WALKED_INTERVALS (128). The
    # first version of that hand-off reported one index at a time, which cost about two node
    # visits per reaching interval where the walk cost one sequential read - so crossing the
    # cap put a step in the wall clock: measured on the lab box, 150k windows at reach 200
    # took 2.36x the same windows at reach 129, against 1.30x for the walk alone. The tree
    # now reports runs of consecutive indices, and the same measurement gives 1.24x.
    #
    # Both sides are the same shape and the same row count; only the window width differs
    # (200bp against 129bp), so the ratio is close to the extra data and nothing else. The
    # test goes quiet rather than false-firing if MAX_WALKED_INTERVALS is ever raised above
    # 200, which is the safe direction.
    n <- 150000L
    st <- seq(0, by = 1, length.out = n)
    dense <- gintervals(1, st, st + 200) # reaches back 199 intervals: past the cap
    subcap <- gintervals(1, st, st + 129) # reaches back 128: the deepest walk-only case

    t_dense <- time_op(gextract("test.fixedbin", dense, iterator = dense))
    t_subcap <- time_op(gextract("test.fixedbin", subcap, iterator = subcap))

    msg <- sprintf(
        "reach 129 (walk): %.0fms, reach 200 (tree): %.0fms (ratio %.2fx)",
        t_subcap * 1000, t_dense * 1000, t_dense / t_subcap
    )
    if (t_dense > t_subcap * 1.8) {
        fail(msg)
    } else {
        succeed(msg)
    }
})