Add glm_pred fused GLM predictor virtual track

NAMESPACE

@@ -99,6 +99,12 @@ export(gintervals.union)
 export(gintervals.update)
 export(giterator.cartesian_grid)
 export(giterator.intervals)
+export(glm_batch_quantiles)
+export(glm_extract_features)
+export(glm_pred.create)
+export(glm_pred.info)
+export(glm_pred.ls)
+export(glm_pred.rm)
 export(glookup)
 export(gpartition)
 export(gquantiles)
@@ -192,11 +198,15 @@ importFrom(utils,write.table)
 useDynLib(misha,C_gcis_decay)
 useDynLib(misha,C_gcompute_strands_autocorr)
 useDynLib(misha,C_gextract)
+useDynLib(misha,C_glm_extract_features)
 useDynLib(misha,C_gpartition)
 useDynLib(misha,C_gquantiles)
+useDynLib(misha,C_gquantiles_multi)
 useDynLib(misha,C_gsample)
 useDynLib(misha,C_gscreen)
+useDynLib(misha,C_gscreen_multi)
 useDynLib(misha,C_gsegment)
+useDynLib(misha,C_gsummary_multi)
 useDynLib(misha,C_gwilcox)
 useDynLib(misha,garrays_import)
 useDynLib(misha,gbins_quantiles)

NEWS.md

@@ -1,3 +1,12 @@
+# misha (development version)
+
+### New features
+
+* New `glm_pred` virtual track type: a fused per-position GLM predictor over an arbitrary number of motif tracks with multi-axis strata defined by `selector_tracks` + `selector_breaks`. `glm_pred.create()` writes the model, `glm_pred.ls/info/rm` manage it, and the track is consumed via the usual `gextract` / `gscreen` paths. Weights/bias/interaction weights are labeled multi-axis arrays (`dim = c(N, K_1, ..., K_M)`); single-selector still accepts a plain `matrix(N, K_1)`.
+* Added `glm_extract_features()`: single C++ call that replaces the R training-side pipeline (chunked `gextract` -> scale -> pivot -> logistic transforms -> GC interactions) with one pass over motif + GC tracks, returning the full `n_peaks x n_features` matrix. Multithreaded with per-chromosome handle reuse.
+* Added `glm_batch_quantiles()`: parallel exact genome-wide quantiles for a batch of motif tracks (e.g. p=0.9999 per-motif caps for `glm_extract_features`). Bypasses `gquantiles` / virtual tracks; ~1.7x faster per track.
+* `gquantiles`, `gsummary`, `gscreen` accept a vector of expressions and return a `data.frame` for multi-track calls (single-expression matrix returns unchanged). An optimized C++ fast path runs automatically for bare-track / simple-vtrack shapes in `gsummary` / `gscreen`; `gquantiles` opts in via `fast = TRUE` (>=5x faster for extreme quantiles on a 10-motif benchmark).
+
 # misha 5.10.2
 
 * **Performance fix:** a single-function `lse` (or `sum`/`exists`/`size`) virtual track scanned over a sliding window (`gvtrack.iterator(sshift=, eshift=)`) genome-wide is fast again. Since 5.6.7 a single-function "fast path" recomputed the windowed reduction from scratch on every step, bypassing the incremental sliding-window path; the common motif-energy quantile workload (windowed `lse` + `gquantiles`/`gscreen`) was ~2.5x slower (worse with wider windows). Such single-function reducers now keep the sliding-window path. Output is unchanged.

R/batch-dispatch.R

@@ -0,0 +1,327 @@
+# Internal helpers for the batched fast-path dispatcher used by gsummary
+# (Phase 3), gscreen (Phase 4), and gquantiles (Phase 6).
+#
+# Contract: detect_fast_path(exprs, iterator, intervals, band) returns
+# either NULL (not fast-path eligible, caller should use slow path) or a
+# list with fields:
+#   $tracks   character vector of underlying track names
+#   $func     "lse" | "avg" | "sum" | "max" | "min"
+#   $sshift   integer
+#   $eshift   integer
+#
+# Window semantics: the C++ scan computes
+#   window = [c + sshift, c + iterator_step + eshift]
+# at each iterator step c, mirroring the slow-path
+# Iterator_modifier1D::transform convention
+#   [interv.start + sshift, interv.end + eshift]
+# where the iterator interval has width iterator_step. Bare tracks pass
+# sshift = eshift = 0 so the window is exactly the iterator step.
+#
+# Preconditions for eligibility:
+#   1. Each expr is a bare track name OR a vtrack wrapping a single
+#      source track with func in {avg, sum, max, min, lse}.
+#   2. All exprs share the same (func, sshift, eshift) tuple.
+#   3. iterator is a fixed integer step (not a track-based iterator).
+#   4. band is NULL.
+#   5. intervals is either NULL, ALLGENOME, or a 1D intervals data.frame.
+
+.describe_single_expr <- function(e) {
+    # Bare track?
+    if (gtrack.exists(e)) {
+        info <- tryCatch(gtrack.info(e), error = function(err) NULL)
+        if (is.null(info)) {
+            return(NULL)
+        }
+        if (!identical(info$type, "dense") && !identical(info$type, "sparse")) {
+            return(NULL)
+        }
+        # Bare track → window equals one iterator step (sshift=eshift=0).
+        # `default_iterator` is the natural step when the caller passes
+        # iterator=NULL: bin_size for a dense track (one bin per step). A
+        # sparse track has no bin grid and its implicit (iterator=NULL)
+        # iterator is the irregular sparse iterator (one position per stored
+        # interval), which the grid-based fast path does NOT reproduce. So a
+        # bare sparse track is only eligible with an EXPLICIT iterator;
+        # default_iterator = NA forces that (mirrors the vtrack rule).
+        if (identical(info$type, "sparse")) {
+            default_it <- NA_integer_
+        } else {
+            bsz <- info$bin.size
+            if (is.null(bsz) || !is.numeric(bsz)) bsz <- 1L
+            default_it <- as.integer(bsz)
+        }
+        return(list(
+            track = e, func = "avg", sshift = 0L, eshift = 0L,
+            default_iterator = default_it
+        ))
+    }
+    # Virtual track?
+    if (exists("GVTRACKS", envir = misha:::.misha)) {
+        gwd <- get("GWD", envir = misha:::.misha)
+        vts <- get("GVTRACKS", envir = misha:::.misha)[[gwd]]
+        if (!is.null(vts) && e %in% names(vts)) {
+            v <- vts[[e]]
+            if (!is.character(v$src) || length(v$src) != 1) {
+                return(NULL)
+            }
+            if (!gtrack.exists(v$src)) {
+                return(NULL)
+            }
+            if (!is.character(v$func) || length(v$func) != 1) {
+                return(NULL)
+            }
+            if (!(v$func %in% c("avg", "sum", "max", "min", "lse"))) {
+                return(NULL)
+            }
+            if (is.null(v$itr) || !identical(v$itr$type, "1d")) {
+                return(NULL)
+            }
+            sshift <- as.integer(v$itr$sshift)
+            eshift <- as.integer(v$itr$eshift)
+            if (is.na(sshift) || is.na(eshift)) {
+                return(NULL)
+            }
+            return(list(
+                track = v$src, func = v$func,
+                sshift = sshift, eshift = eshift,
+                default_iterator = NA_integer_
+            ))
+        }
+    }
+    NULL
+}
+
+# TRUE iff every scope interval sits on the iterator grid (origin 0), so the
+# fast path's global-grid + point-membership scan matches the slow path's
+# clipped fixed-bin iterator. A non-aligned start or a non-aligned end that is
+# NOT the chromosome end would shift or truncate a partial bin differently
+# between the two paths, so such scopes must fall back to the slow path. An end
+# equal to the chromosome size is allowed: both paths clip the final window to
+# the chromosome boundary identically.
+.scope_grid_aligned <- function(iv, it_int) {
+    allg <- get("ALLGENOME", envir = misha:::.misha)[[1]]
+    chrom_size <- stats::setNames(as.numeric(allg$end), as.character(allg$chrom))
+    starts <- as.numeric(iv$start)
+    ends <- as.numeric(iv$end)
+    csz <- chrom_size[as.character(iv$chrom)]
+    all(starts %% it_int == 0) &&
+        all((ends %% it_int == 0) | (!is.na(csz) & ends == csz))
+}
+
+.detect_fast_path <- function(exprs, iterator, intervals, band) {
+    if (!is.character(exprs) || length(exprs) == 0) {
+        return(NULL)
+    }
+    if (!is.null(band)) {
+        return(NULL)
+    }
+    # Cheap lexical gate: skip the gtrack.info / .gvtrack.get calls when
+    # any expression contains operators or whitespace.
+    if (!all(.looks_like_bare_name(exprs))) {
+        return(NULL)
+    }
+
+    infos <- lapply(exprs, .describe_single_expr)
+    if (any(vapply(infos, is.null, logical(1)))) {
+        return(NULL)
+    }
+
+    # iterator is optional when every expression is a bare track with a
+    # known bin size — default to the common bin size (all must match).
+    # Vtracks have default_iterator = NA, so any vtrack forces the user
+    # to specify iterator explicitly.
+    if (is.null(iterator)) {
+        defaults <- vapply(infos, `[[`, integer(1), "default_iterator")
+        if (any(is.na(defaults))) {
+            return(NULL)
+        }
+        if (length(unique(defaults)) != 1) {
+            return(NULL)
+        }
+        it_int <- defaults[1]
+    } else {
+        if (!is.numeric(iterator) || length(iterator) != 1) {
+            return(NULL)
+        }
+        it_int <- as.integer(iterator)
+        if (is.na(it_int) || it_int <= 0) {
+            return(NULL)
+        }
+    }
+
+    funcs <- vapply(infos, `[[`, character(1), "func")
+    sshift <- vapply(infos, `[[`, integer(1), "sshift")
+    eshift <- vapply(infos, `[[`, integer(1), "eshift")
+
+    if (length(unique(funcs)) != 1) {
+        return(NULL)
+    }
+    if (length(unique(sshift)) != 1) {
+        return(NULL)
+    }
+    if (length(unique(eshift)) != 1) {
+        return(NULL)
+    }
+
+    # Intervals: accept 1D data.frame. ALLGENOME is a list of
+    # (1D_df, 2D_df); unwrap to the 1D part. Reject bigset handles
+    # (character strings) and anything 2D.
+    if (!is.null(intervals)) {
+        iv <- intervals
+        if (is.list(iv) && !is.data.frame(iv) && length(iv) == 2 &&
+            is.data.frame(iv[[1]])) {
+            iv <- iv[[1]]
+        }
+        if (!is.data.frame(iv)) {
+            return(NULL)
+        }
+        if (!all(c("chrom", "start", "end") %in% colnames(iv))) {
+            return(NULL)
+        }
+        if ("chrom1" %in% colnames(iv)) {
+            return(NULL)
+        }
+        # Only grid-aligned scopes match the slow path (see .scope_grid_aligned).
+        if (!.scope_grid_aligned(iv, it_int)) {
+            return(NULL)
+        }
+    }
+
+    list(
+        tracks = vapply(infos, `[[`, character(1), "track"),
+        func = funcs[1],
+        sshift = sshift[1],
+        eshift = eshift[1],
+        iterator = it_int
+    )
+}
+
+
+# Screen-specific parser. Each expression must be a single comparison:
+#   "<lhs> <op> <const>"  where op in <, <=, ==, >=, >
+# and <lhs> must satisfy .describe_single_expr (bare track or simple
+# vtrack). Returns either NULL or a list with tracks/func/sshift/eshift/
+# iterator/ops/thresholds fields.
+.detect_screen_fast_path <- function(exprs, iterator, intervals, band) {
+    if (!is.character(exprs) || length(exprs) == 0) {
+        return(NULL)
+    }
+    if (!is.null(band)) {
+        return(NULL)
+    }
+    # Parse "<lhs> <op> <const>" from each expr.
+    rx <- "^\\s*(.+?)\\s*(<=|>=|==|<|>)\\s*([-+0-9.eE]+)\\s*$"
+    m <- regmatches(exprs, regexec(rx, exprs))
+    if (any(vapply(m, function(x) length(x) != 4, logical(1)))) {
+        return(NULL)
+    }
+    lhs <- vapply(m, `[`, character(1), 2)
+    ops <- vapply(m, `[`, character(1), 3)
+    thr <- suppressWarnings(as.numeric(vapply(m, `[`, character(1), 4)))
+    if (any(is.na(thr))) {
+        return(NULL)
+    }
+
+    infos <- lapply(lhs, .describe_single_expr)
+    if (any(vapply(infos, is.null, logical(1)))) {
+        return(NULL)
+    }
+
+    funcs <- vapply(infos, `[[`, character(1), "func")
+    sshifts <- vapply(infos, `[[`, integer(1), "sshift")
+    eshifts <- vapply(infos, `[[`, integer(1), "eshift")
+
+    if (is.null(iterator)) {
+        defaults <- vapply(infos, `[[`, integer(1), "default_iterator")
+        if (any(is.na(defaults))) {
+            return(NULL)
+        }
+        if (length(unique(defaults)) != 1) {
+            return(NULL)
+        }
+        it_int <- defaults[1]
+    } else {
+        if (!is.numeric(iterator) || length(iterator) != 1) {
+            return(NULL)
+        }
+        it_int <- as.integer(iterator)
+        if (is.na(it_int) || it_int <= 0) {
+            return(NULL)
+        }
+    }
+
+    if (length(unique(funcs)) != 1) {
+        return(NULL)
+    }
+    if (length(unique(sshifts)) != 1) {
+        return(NULL)
+    }
+    if (length(unique(eshifts)) != 1) {
+        return(NULL)
+    }
+
+    if (!is.null(intervals)) {
+        iv <- intervals
+        if (is.list(iv) && !is.data.frame(iv) && length(iv) == 2 &&
+            is.data.frame(iv[[1]])) {
+            iv <- iv[[1]]
+        }
+        if (!is.data.frame(iv)) {
+            return(NULL)
+        }
+        if (!all(c("chrom", "start", "end") %in% colnames(iv))) {
+            return(NULL)
+        }
+        if ("chrom1" %in% colnames(iv)) {
+            return(NULL)
+        }
+        # Only grid-aligned scopes match the slow path (see .scope_grid_aligned).
+        if (!.scope_grid_aligned(iv, it_int)) {
+            return(NULL)
+        }
+    }
+
+    list(
+        tracks = vapply(infos, `[[`, character(1), "track"),
+        func = funcs[1],
+        sshift = sshifts[1],
+        eshift = eshifts[1],
+        iterator = it_int,
+        ops = ops,
+        thresholds = thr
+    )
+}
+
+# Map comparison operator strings to the C-side CmpOp enum integers
+# (see ThresholdScreen::CmpOp in BatchScreen.cpp).
+.screen_op_to_int <- function(op_chr) {
+    tbl <- c("<" = 0L, "<=" = 1L, "==" = 2L, ">=" = 3L, ">" = 4L)
+    out <- tbl[op_chr]
+    if (any(is.na(out))) stop("unknown comparison operator", call. = FALSE)
+    as.integer(out)
+}
+
+# Cheap lexical pre-check: an expression that contains operators,
+# whitespace-around-operators, or parentheses cannot be a bare track or
+# vtrack name, so skip the expensive .describe_single_expr path entirely.
+# Matches identifier-like strings (letters, digits, "." and "_"), optionally
+# wrapped in outer whitespace. Anything else falls through immediately.
+.looks_like_bare_name <- function(s) {
+    grepl("^\\s*[A-Za-z0-9._]+\\s*$", s)
+}
+
+# Emit a one-time per-session informational message explaining which path
+# was taken (or why the fast path was declined). Suppressed via
+# options(misha.quiet_dispatch = TRUE).
+.fast_dispatch_msg <- function(fn, reason) {
+    if (isTRUE(getOption("misha.quiet_dispatch"))) {
+        return(invisible())
+    }
+    key <- paste0("misha.dispatch_msg.", fn)
+    if (isTRUE(getOption(key))) {
+        return(invisible())
+    }
+    packageStartupMessage(sprintf("[%s] %s", fn, reason))
+    args <- setNames(list(TRUE), key)
+    do.call(options, args)
+}