The Scoring Engine

The compute chapter shows compute$scores() as two lines that return a matrix. That is the right level for most analyses. This chapter is for when the two lines are not enough on their own — a model library large enough that the call needs a resource estimate first, a job that fails with an out-of-memory error, or a cluster allocation you have to size before submitting anything. It documents the mechanism: what compute$scores()'s memory argument actually controls, what compute$scores.plan() tells you before you commit to a real run, and what determines a batch.

What a score request has to resolve

Every model contributes its own list of variants. compute$scores() first collates the union of variants across every model being scored, then matches each one to a genotype record by chromosome and position, accepting either allele orientation. PLINK2's --score then multiplies dosages by effect sizes and sums them into one PRS per model, using the score file's effect-allele column to handle orientation — no dosage is ever flipped in R.

Two things follow from doing this the way PolyGenius does. First, this whole path is driven through study$genotypes$scores(), which is why the object you call compute$scores() on needs a model library and at least one genotype attached, not handed to it as loose arguments. Second, the union is scored either fused — every model in one PLINK pass over one score file — or through a batched via-extract path, and which one you get depends on the workload and the knobs described below.

Why sharing a backbone bounds scoring cost

The previous section described a union over "every model being scored" as if that were a simple pass. It is worth being precise about what makes it simple, because the alternative — the naive way to compute M scores against a shared genotype file — is a per-model loop, and a per-model loop over a large library is the single biggest cost trap in this engine's history.

compute$scores() never scans the model library one model at a time. Every shape you can pass as models — a PGS, a list of them, a raw data.frame, a PGSLibrary — is normalized into one backbone-backed PGSLibrary before any scoring logic runs at all. By the time the union is being built, every model's variants are already rows in the shared dictionary; there is no separate per-model variant table left to scan. See Data Architecture for what that dictionary and its weight store actually look like — this section is about what that shape buys you at scoring time.

Mechanically: compute.scores() routes to scores.model.cursor() (R/genotype-source-scores-resolve.R), which populates cursor$backbone and cursor$model.slots with no per-model materialization — a model is a slot index into the backbone, not a table. GenotypeSourceScores$compute() reads that field and dispatches to scores.index.backbone() (R/genotype-source-scores-store.R), which reads each selected model's (dict.slot, beta) pair straight off backbone$weights$get(model.slots[[i]]). No per-model (chr, position, nea, ea) table is ever built, and no codec call happens on this path at all — the packed variant key described later in this chapter (scores.pack.key(), R/genotype-source-scores-key.R) belongs to the non-backbone table-cursor path only. The backbone-backed path dispenses with keying for dedup purposes entirely: scores.index.backbone() gathers the selected dict.slots into a dictionary-wide hit-flag array — one raw() byte per dictionary row (not sub-byte bit-packing, so 4x smaller than the same array as logical()'s 4-byte elements, not the ~32x a true bit-per-flag structure would give). A variant already flagged by an earlier model just flips a flag that is already set. It is never re-inserted, re-keyed, or given a second row, no matter how many of the M models reference it. dict.rank (see Data Architecture) gives this array's row order, but the array itself is scoped to one call: the per-model (variant.idx, beta) records below use a compacted variant.idx numbering just the variants this call actually selected, not dict.rank itself.

That flag array is the whole mechanism behind the cost claim: with sharing, building the union and the per-model records costs O(V_dict + sum of model sizes), never O(M x V). What's held in memory is one dictionary-wide flag array plus one compact (variant.idx, beta) record per model — never a matrix sized models x variants. This is also why the file-backed regime exists at all: above file.backed.threshold (a cheap backbone$n.variants() x M upper bound, checked before any real indexing by scores.file.backed.decision()), the same union/record derivation goes out-of-core via scores.index.backbone.blocks(), walking a FileBlocks/MemoryBlocks shard window at a time instead of holding every model's column resident — a streaming version of the same O(V + M) shape, not a different algorithm.

This is the same class of blowup this design was built to eliminate: a per-model loop that redundantly re-fetches or re-decodes shared state turns an O(M) cost into O(M^2)/O(M x V). An earlier per-model-loop design — the shape backbone-sharing replaced — hit exactly this, measured: a benchmark cell went from 721.7s to 8068.2s (more than 10x slower) between M=300 and M=1000 models, a mere 3.3x increase in model count. That number is the old design's failure, not a re-measurement of the current flag-array-and-record design on the same workload — a benchmark comparing the two at the same M=300/M=1000 sizes is defined (benchmark.scores/, tracking the "M-axis" sweep) but has not completed as of this writing. See .claude/references/scoring-cost-model.md §7 for the full catalog of patterns that turn an innocent-looking per-model loop into that shape, scoring-data-structures.md §1–§2 for the flag-array and record layout in more detail, and scoring-invariants.md's "never build a dense M x V structure" rule for why this is enforced as an invariant rather than left as a current implementation detail.

memory: one argument, three knobs

compute$scores() exposes a single memory argument, memory. It resolves into three concrete numbers the engine actually reads:

• cells — the batch-sizing cap, in cells (V_batch × M_batch), not bytes. This is the knob that decides how many PLINK --score calls a run makes.
• plink — PLINK's own --memory flag, in MB, capping one PLINK child process. A real, concurrent reservation, not a sizing heuristic.
• file.backed — a cell-count threshold (backbone variants × models), not a reservation. Above it, a backbone-backed model library has its weight store read out-of-core instead of held resident; below it, the whole weight matrix stays in memory for the call.

You never set these three directly — max.cells.per.batch, plink.memory and file.backed.threshold do not exist on compute$scores()'s signature at all. memory is what you pass, and it accepts three shapes:

compute$scores(data)                                   # memory = NULL (default)
compute$scores(data, memory = 8)                        # an 8 GiB total budget
compute$scores(data, memory = list(cells = 5e7))        # override one knob directly

memory = NULL falls back to workspace$config$max.memory, the package-wide admission budget in GiB (Inf by default). With no budget in force anywhere, the three knobs take the historical literal values byte-for-byte: cells = 200e6, plink = NULL (PLINK's own default, uncapped), file.backed = 2.5e8.

A bare numeric scalar is a total GiB budget of its own, overriding the package-wide one. A named list can override any subset of total (GiB), plink (MB), cells (integer cell count) and file.backed (integer cell count); anything you supply explicitly wins verbatim over deriving it from total. When a total budget is in force, plink — a genuine concurrent reservation — is subtracted from total before cells is derived from what remains. file.backed is a decision threshold rather than a reservation, so it is always derived from the full, un-subtracted total. This ordering is the one thing worth remembering about memory: a plink reservation shrinks the batch cap you get, but never shrinks the point at which the engine decides to go out-of-core.

workspace$config$max.memory is otherwise not enforced here — scoring never goes through the execution engine's own admission path (see the execution engine for where that path does apply).

Batching: what determines a batch, and what "peak RAM" bounds

Once cells is resolved, the union of scoreable variants is partitioned into contiguous windows of width W = floor(cells / M), where M is the number of models being scored. Each window is one batch, and because variant.idx orders variants genomically, a window is a genomic range — which also happens to cluster a batch's variants by source file. A model whose variants span more than one window has its partial sums added downstream, which is exact because a PRS is a sum.

Two consequences worth knowing. W is sized against the global model count, not how many models actually fall in a given window — so a library of small, region-confined models (a handful of single-locus scores, say) under-fills its windows and produces more, narrower batches than a genome-wide library of the same total variant count would. And lowering cells amplifies that: it caps peak memory, but at the cost of more PLINK passes, so it is a memory/throughput trade, not a free lever.

What "peak RAM" actually bounds is the per-batch score-matrix transient — V_batch × M_batch cells, in both the fused and via-extract regimes, because each strategy still writes one score file per source file (this batch intersected with that file's variants). The fused path is only available when the whole union scores in one batch and no minor-allele-frequency filter is in force; anything else takes the via-extract path. Setting maf.thr > 0 therefore always forfeits fusion, but it does not add a second genotype read — the frequency computation is folded into the same extraction pass that resolves variants, so a large imputed fileset is still read once, not once for --freq and again to extract it.

Gotcha: a model filter forces the in-memory path

model.filter is normally cheap in exactly the way you'd hope: it is pushed down as a mask (scores.model.mask.set(), R/genotype-source-scores-store.R) rather than materializing a filtered copy of the library. models itself stays the original, unfiltered PGSLibrary throughout — the filter travels as a side-channel of masked (dict.slot, beta) columns, not a second copy of the PGS library.

But a supplied filter forces the in-memory regime unconditionally, regardless of how large the unfiltered library is. scores.index.backbone.blocks() — the out-of-core reader — has no filter-pushdown mechanism at all. There is nowhere in that code path to apply a mask before deciding which shard windows to read, so the moment you pass a non-default model.filter, the engine falls back to scores.index.backbone(), the resident path, no matter what file.backed.threshold says.

This is a real practical trap, not a hypothetical one. "Filter down to what you need" is good general advice, and for a small-to-medium library it is also free. But at genome scale it inverts: if you have a library large enough to trip file.backed.threshold on its own (say, tens of thousands of models against a multi-million-row dictionary), and you filter it down to a handful of models you actually want scored, that filtered call does not stream. It holds the filtered result fully resident in memory — the in-memory path runs on whatever the filter lets through, not on the full unfiltered library, but it is still the in-memory path, with no windowed reads at all. If you are scoring a small subset of a very large library and are relying on file.backed to keep memory bounded, a supplied model.filter silently removes that guarantee.

Be careful trusting compute$scores.plan() here, too: its file-backed decision (est.for() in R/compute-scores-plan.R) is computed from dims$n.variants.dict, which is always backbone$n.variants() — the full, unfiltered dictionary size — with no filter.is.default check anywhere in that path. A plan run with a model.filter supplied can still report strategy = "file.backed", because the plan's threshold check does not know that a supplied filter forces the real call into the in-memory regime. The plan and the real dispatch can disagree on this specific point; treat a filtered plan's "file.backed" strategy as unconfirmed until you've checked memory on an actual run.

compute$scores.plan(): an estimate, never a score matrix

compute$scores.plan() is the dedicated way to find out what a real call would cost before running it. It takes the same study, maf.thr, model.filter, models and memory arguments as compute$scores(), and never runs PLINK and never returns a score matrix — only ever a resource estimate, or a table of them.

plan <- compute$scores.plan(data)
plan$strategy       # "fused", "via.extract" or "file.backed"
plan$n.batches
plan$peak.ram.mb

The single-estimate return carries peak.ram.mb, n.batches, n.groups, strategy, n.models, n.variants.union, n.variants.total, n.samples, the est.* cost-model breakdown terms, and a memory element — list(requested, resolved, source) — recording exactly what your memory argument resolved to and whether each of the three knobs was "literal", "derived" or "supplied". strategy is the one place "file.backed" appears as a value in its own right: a real compute$scores() call's returned strategy only ever reports "fused"/"via.extract", because which indexer built the variant union is invisible to the caller once scoring is done.

exact controls how the variant union is counted: exact = TRUE (the default) computes the true scored union with one streaming pass over the models and no genotype access; exact = FALSE uses the cheap dense upper bound n.variants.union = n.variants.total — useful for a first-glance check, pessimistic whenever models share variants.

sweep lets you compare several memory budgets against the same measured workload in one call, without re-measuring it per candidate:

compute$scores.plan(data, sweep = c(4, 8, 16, 32))
compute$scores.plan(data, sweep = list(2, list(cells = 5e7), list(total = 8, plink = 2000)))

This returns a data.table, one row per candidate, with label, strategy, n.batches, peak.ram.mb, cells, plink.memory and file.backed.threshold — the table to read before choosing an allocation, rather than guessing and re-submitting.

No wall-clock estimate is provided anywhere in this output. The underlying cost model was fitted against measured peak memory, not against measured wall time, so n.batches is the closest available time-correlated proxy; treat it that way rather than looking for a duration that was never fitted.

The variant-fate record

Every scored matrix carries a fate record: what became of each variant in the requested union — scored, missing, maf.filtered, no.freq or collision — together with whether its alleles were flipped and whether it is strand-ambiguous. It sits at diagnostics(scores, "fate"), beside the coverage diagnostics.

scores <- compute$scores(data, maf.thr = 0.01)
fate   <- diagnostics(scores, "fate")
 
fate                          # counts, plus the first few variants
as.data.table(fate)           # one row per requested variant
as.data.table(fate)[status != "scored"]   # only the ones that did not make it

Printing shows the counts and a short head. as.data.table() and as.data.frame() expand the record into one row per requested variant per genotype, carrying chr, position, nea, ea, status, matched, flipped, strand.ambiguous, freq and n.records.

There is no retention setting to choose. Every requested variant keeps its detail, because the record stores a single byte per variant per genotype beside one copy of the variant identities, rather than the ten-column table it expands into. That matters because the record scales with the size of the requested union — not with how many models you scored — and it travels inside the object every time you save it. Measured at 1.1 million requested variants it serializes to 19 MB against an 8.3 MB score matrix, and to 174 MB at 10 million; the table it expands to would be 62 MB and 568 MB.

Two things are reconstructed rather than stored, so they cost nothing: matched, which is simply status != "missing", and n.records, which is 0 or 1 except at a collision. A scored variant's freq is not kept either — at the default maf.thr = 0 no frequency pass runs, so there is none to keep — while every excluded variant keeps its own.

With several genotype filesets the counts are over unique variants, and a variant scored in any one of them counts as scored: the multi-fileset score is a sum, so it contributes. fate$counts.by.genotype holds the per-genotype breakdown.

What this does not do

This page is about PolyGenius's own scoring mechanism, not PLINK2's --score flag in general — for PLINK2 mechanics beyond what PolyGenius drives, that is the PLINK2 documentation and plink2-expert's territory, not this one.

The scoring path does not go through the execution engine described in the execution engine chapter — no rule matching, no scheduler, no store-backed caching of a score matrix itself. workspace$config$max.memory is reused only as memory's own fallback default (see above); the engine's own memory-admission check, which refuses to schedule a task whose declared requirement exceeds the budget, belongs to generation and never runs for a compute$scores() call.

compute$scores.plan() never validates that the estimate it returns is achievable on your machine — it predicts what a real call would need, it does not check what your machine has. Compare peak.ram.mb against your own allocation yourself.

It hoists nothing across calls. Scoring several cohorts is one compute$scores() call per study, and each study owns its own model library, so there is no shared state a second call could reuse even in principle. Every call re-resolves its PGS library — .compute.scores.materialize() and scores.model.mask.set() run from scratch, and variants.for.PGSLibraryView() reads straight off the weight store with no result cache anywhere in PGSLibrary/PGSLibraryView/PolyGeniusBackbone to short-circuit a repeat. When a call trips the file-backed regime it also re-reads the weight store window by window, because scores.index.backbone.blocks() lives inside a single call's body.

That cost is per call and scales with the number of cohorts, but it only bites once the library itself is large enough to be file-backed — the whole point of file.backed.threshold. A handful of cohorts against a library in the hundreds of models, the common case, stays in the resident-weight-store regime, where per-call resolution is cheap.

Where to go next

Computing scores and covariates is the everyday entry point this chapter supports. The execution engine covers the separate scheduling and caching layer that generation (not scoring) runs through.