"""Reducing a region to one value per requested score.
The machinery the per-kind aggregating reads share --
:meth:`~.position.PositionScore.get_scores_in_region_agg`,
:meth:`~.allele.AlleleScore.get_allele_scores_in_region_agg` and
:meth:`~.fragment.FragmentScore.get_fragment_scores_overlapping_region_agg`
-- apart from any score class: resolving a caller's queries to
``(score_id, aggregator)`` pairs, building a fresh aggregator for each, and
folding a stream of fetched segments into one value per request.
:func:`fold_into_bins` is the same reduction per grid bin rather than per
region: a dense block, one row a bin and one column an aggregator, with the
empty-bin rule (:data:`EMPTY_BIN_VALUES`) stated there and nowhere else.
Resolving and building are separate steps here because they are separate
questions. :func:`score_def_for` and :func:`resolve_aggregator_name`
answer "is this request answerable, and by what", which a caller may ask
without meaning to read -- since gain#1131
:meth:`~.position.PositionScore.resolve_aggregation_queries` asks exactly
that when an annotation pipeline loads, so a misconfigured attribute is
refused before any annotatable arrives. :func:`build_region_aggregator`
answers "give me somewhere to accumulate", which only a read needs, and
which a read needs freshly every time.
Nothing here knows what a score KIND is. The one thing that differs
between kinds -- how many times a record's value counts -- reaches
:func:`fold_region_segments` as ``weigh``, the kind's own
:meth:`~.base.GenomicScore.record_weight`. The fold applies what it is
handed and does not re-derive or cross-check it: the rule is stated once
per kind, on the kind (see
``test_the_weight_rule_is_stated_once_per_kind``). Whether the segments
were clipped to the query window before they got here is settled the same
way, by the kind's own read, so the fold carries no flag saying which kind
it serves.
Every BACKEND feeds one path here -- this generic weighted stream over
fetched records -- and a backend-specific fast path was considered and
rejected to keep it that way; see
``.out-of-scope/bigwig-stats-pushdown.md``. That is a claim about
backends, not about the package: :class:`~.position.PositionScore` folds
its own segments for the aggregated plane, and converging the two is
gain#1027's remaining work, not this module's promise.
Every refusal here is of a REQUEST, never of a definition: a score the
resource does not define, an aggregator a caller names that does not
build, a ``bool`` score asked for with no aggregator named. The
definitions a request is resolved against were held to their own rules
when the score was built (a ``scores:`` ``aggregator:`` that cannot be
built is refused there, through ``score_configuration_error``), so by the
time a request arrives every definition is one the resource can reduce
by. That is why these refusals carry the request family's address --
which score, of which resource -- and not the configuration prefix: the
file to edit is the caller's, and ``Invalid configuration: <resource>``
would send the reader to a yaml that is fine. The undefined-score
sentence is not worded here: every surface that refuses that mistake
shares it, from
:func:`~gain.genomic_resources.resource_errors.undefined_scores_message`.
"""
from __future__ import annotations
import math
from collections.abc import Callable, Iterable, Sequence
import numpy as np
from gain.genomic_resources.resource_errors import undefined_scores_message
from gain.genomic_resources.score_def import GenomicScoreDef, ScoreValue
from gain.utils.regions import calc_bin_index
from ..aggregators import Aggregator, ScoreAggregationQuery, to_float64
# How a query surface tells a caller to name an aggregator the score has
# no default for. It lives HERE, with the rule it is appended to, so the
# whole sentence a caller sees is written in this module. A fragment:
# lowercase, no leading punctuation, spliced after the rule's "; ".
QUERY_AGGREGATOR_REMEDY = "name one on the query"
#: The value that one aggregator reduces a region to. The aggregator sets
#: the type: ``list`` and ``most_common`` give a list, ``value_count`` gives
#: a dict that maps each value to its count, and every other aggregator
#: gives a score value.
AggregateValue = ScoreValue | list[ScoreValue] | dict[ScoreValue, int]
[docs]
def resolve_aggregation_queries(
queries: Sequence[ScoreAggregationQuery] | None,
*,
score_definitions: dict[str, GenomicScoreDef],
all_scores: list[str],
resource_id: str,
) -> list[tuple[str, str]]:
"""Resolve kind-neutral queries to ``(score_id, aggregator NAME)`` pairs.
Asks the two questions every aggregation request asks
(:func:`score_def_for`, :func:`resolve_aggregator_name`), and expands
``None`` to every score with its own default. Stops at the NAME, as
:meth:`~.position.PositionScore.resolve_aggregation_queries` does, so
a caller can ask whether a query list is answerable without reading --
an annotator asks exactly that as its pipeline loads. Building the
accumulators is the READ's business, per call
(:func:`build_region_aggregator`).
Shared by the fragment and allele folding reads, which is what
promoted it here from the fragment kind: this package keeps a
derivation on the one kind that needs it and moves it here when a
second does. The position kind keeps its own resolver because its
query carries a third field, ``none_value_replacement``, that has to
be judged between the two questions -- see that method for why the
two must not be merged. A pair has nowhere to put that field, and
the signature keeps its query out: a
:class:`~..aggregators.PositionScoreAggregationQuery` is a sibling of
the neutral query, not a subclass, so it does not type-check here --
see :class:`~..aggregators.ScoreAggregationQuery` for why (gain#1302).
"""
if queries is None:
queries = [
ScoreAggregationQuery(score_id) for score_id in all_scores
]
requests = []
for query in queries:
score_def = score_def_for(
query.score,
score_definitions=score_definitions,
resource_id=resource_id)
requests.append((
query.score,
resolve_aggregator_name(
query.aggregator, score_def,
resource_id=resource_id)))
return requests
[docs]
def score_def_for(
score_id: str,
*,
score_definitions: dict[str, GenomicScoreDef],
resource_id: str,
) -> GenomicScoreDef:
"""The definition an aggregation request names, refusing an unknown one.
The first of the two questions every aggregation request asks, and the
one place aggregation refuses when the answer is no. Whether the
request arrived as a bare score id, as a ``(score_id, aggregator)``
pair, as a
:class:`~gain.genomic_resources.aggregators.ScoreAggregationQuery` or
as a
:class:`~gain.genomic_resources.aggregators.PositionScoreAggregationQuery`
changes nothing about it: the resource either defines that score or it
does not, and the caller is told which ones it has either way.
"""
score_def = score_definitions.get(score_id)
if score_def is None:
raise ValueError(undefined_scores_message(
resource_id, [score_id], score_definitions))
return score_def
[docs]
def resolve_aggregator_name(
aggregator: str | None,
score_def: GenomicScoreDef,
*,
resource_id: str,
) -> str:
"""The aggregator to reduce a score with: the caller's, else its own.
The second question, and the one statement of the rule that a score
with neither is refused. Every surface takes the aggregator as a
field on its query, so every refusal ends with
:data:`QUERY_AGGREGATOR_REMEDY`. The neutral resolver and the
position score's resolver share the whole sentence, and
``test_both_surfaces_state_the_missing_default_rule_identically``
pins that.
The score is named by ``score_def`` rather than beside it:
``score_definitions`` is keyed by ``score_id`` at every construction
path, so a separate argument would be one the caller could contradict.
"""
resolved = aggregator or score_def.aggregator
if resolved is None:
# Every score has a value type, so the only way to get here is a
# type whose class default is deliberately None -- ``bool``, which
# has no meaningful reduction to pick for the caller. Name one
# and it works.
raise ValueError(
f"score {score_def.score_id!r} of resource {resource_id!r} "
f"has no default aggregator for value type "
f"{score_def.value_type!r}; {QUERY_AGGREGATOR_REMEDY}")
return resolved
[docs]
def distinct_score_ids(score_ids: Iterable[str]) -> list[str]:
"""The DISTINCT ids among ``score_ids``, in the order asked for.
One fetch serves every aggregation request, so the same list must both
name what is fetched and index the values that come back -- which is
why the ORDER is part of the answer, and why every aggregating read in
this package derives it here. A second spelling that ordered the
scores differently would not fail; it would have every aggregator
quietly reading its neighbour's column.
Takes the ids rather than the requests they came off, because a
request is shaped differently on each surface -- a pair here, a
:class:`~gain.genomic_resources.aggregators.PositionScoreAggregationQuery`
resolved to a triple on the position score's plane -- and none of that
is what the derivation is about. Each surface projects its own shape
at its call site; the request list gets a named projection,
:func:`request_score_ids`, because it is the one shape TWO readers
(the allele and fragment kinds' folding reads) project, and each hands
the result on to the fold.
Note this is the aggregating reads' derivation, not a package-wide
one: nothing else in the package fetches one list and indexes it
twice. The position annotator's point read asks one score per
attribute, a score named twice included, and pairs the answers back
by position, so it has nothing to dedupe.
"""
return list(dict.fromkeys(score_ids))
[docs]
def request_score_ids(requests: list[tuple[str, str]]) -> list[str]:
""":func:`distinct_score_ids` of a request list's scores.
Derived ONCE per aggregating read and handed to both ends of it: the
reader names the scores to fetch with it and passes the same list to
:func:`fold_region_segments` as ``score_ids``, which says why one
list serves both. The position score's plane projects its own shape
inline instead -- one reader, nothing to hand on.
"""
return distinct_score_ids(score_id for score_id, _ in requests)
[docs]
def build_region_aggregator(
score_id: str, aggregator: str, *, resource_id: str,
) -> Aggregator:
"""Build a FRESH aggregator, naming the resource if it cannot.
Fresh per call, not reused: an aggregator is a mutable accumulator
and explicitly not thread-safe, and every caller in gain builds one
per fold (see
:class:`~gain.genomic_resources.aggregators.Aggregator`). A score
holds none between reads.
``Aggregator.build`` raises a bare ``KeyError('mediann')`` for an
unknown name, saying nothing about which score asked for it. The
name refused is the CALLER's: a definition's own ``aggregator:`` was
built once when the score was constructed
(``score_def.refuse_unbuildable_aggregators``), so a resolved
request never carries a definition-side spelling that fails here.
"""
try:
return Aggregator.build(aggregator)
except (KeyError, ValueError, TypeError) as err:
raise ValueError(
f"score {score_id!r} of resource {resource_id!r} asks "
f"for aggregator {aggregator!r}, which is not valid: "
f"{err}") from err
[docs]
def build_region_aggregators(
requests: Sequence[tuple[str, str]], *, resource_id: str,
) -> list[Aggregator]:
"""One FRESH aggregator per request, parallel to the request list.
:func:`build_region_aggregator` over a resolved request list -- the
shape every folding read hands :func:`fold_region_segments`, built
once per READ and never held on the score, for the reason that
function gives. Fresh accumulators remove one hazard, not the class
of them: two concurrent region reads of one open score still share
the table's line iterator, which is each read's "one live region read
at a time" rule and not this function's promise.
"""
return [
build_region_aggregator(
score_id, aggregator, resource_id=resource_id)
for score_id, aggregator in requests
]
[docs]
def fold_region_segments(
segments: Iterable[tuple[int, int, Sequence[ScoreValue]]],
aggregators: list[Aggregator],
requests: Sequence[tuple[str, str]],
*,
score_ids: Sequence[str],
weigh: Callable[[int, int], int],
) -> list[AggregateValue]:
"""Fold one region read into one value per request.
``segments`` is a stream of ``(left, right, values)`` as
:meth:`~.base.GenomicScore.fetch_region_segments_scores` yields it, with
``values`` positional and parallel to ``score_ids`` -- the list the
caller derived (:func:`request_score_ids`) to name what it fetched,
handed on here so the fold indexes exactly the columns the fetch
carried; one derivation cannot disagree with itself, where two (the
fold once derived its own) had to be kept agreeing. That is how a
request finds its column: two requests for one score share the fetch
and keep separate accumulators. Any
``Sequence`` of values will do -- the fold only ever indexes
``values[column]`` -- so a kind that hands over tuples folds exactly
as one that hands over lists.
``aggregators`` is parallel to ``requests`` and built by the CALLER,
which is what lets an invalid aggregator name be refused before the
region is read at all: a read that refuses an unknown contig when it is
called would otherwise mask a misspelled aggregator behind that
complaint.
``weigh`` is the caller's per-kind weight rule, applied as handed: it
turns a segment's span into the number of times that record's value
counts. It is not second-guessed here, and neither is the stream --
whether a record was first cut down to the query window is the
caller's business, settled before the segments arrive. A kind that
counts a record once counts it wherever the point it collapses to
falls, window or not.
"""
column_of = {score_id: i for i, score_id in enumerate(score_ids)}
targets = [
(aggregator, column_of[score_id])
for aggregator, (score_id, _) in zip(
aggregators, requests, strict=True)
]
for left, right, values in segments:
weight = weigh(left, right)
for aggregator, column in targets:
aggregator.add(values[column], weight)
return [aggregator.get_final() for aggregator in aggregators]
#: What a cell of :func:`fold_into_bins` holds when its aggregator received
#: no non-null value, per aggregator: ``0`` where nothing summed or counted
#: IS zero, ``NaN`` where there is nothing to answer. The keys are also the
#: whole set of aggregators the fold accepts -- the ones whose answer is a
#: float. This is the fragment surfaces' rule (a bin with no fragment has
#: nothing to fold), and it differs on purpose from the position-score
#: binner's, which stores NaN for an uncovered bin unless the track's
#: ``none_value_replacement`` says otherwise; that key does not exist on
#: the fragment kind. The two rules agree on ``count``: its aggregator
#: answers ``0`` for no values, so the position-score binner stores 0.0
#: for an uncovered bin of a ``count`` track too.
EMPTY_BIN_VALUES: dict[str, float] = {
"count": 0.0,
"sum": 0.0,
"mean": math.nan,
"max": math.nan,
"min": math.nan,
"median": math.nan,
"product": math.nan,
}
[docs]
def fold_into_bins(
records: Iterable[tuple[int, int, ScoreValue]],
*,
start: int,
end: int,
bin_size: int,
aggregators: Sequence[str],
) -> np.ndarray:
"""Fold position-sorted records into a dense block of global-grid bins.
``records`` are ``(position, column, value)``: each value is added to
a fresh aggregator of its ``column`` in the bin holding its
``position``. ``aggregators`` names one aggregator per column, which
is how a caller gives each column its meaning -- one per query for a
score's binned read, one per group for a binner.
The answer is a ``float64`` array of shape ``(n_bins, n_columns)``.
Row ``i`` is bin ``calc_bin_index(bin_size, start) + i`` of the grid
anchored at position 1, so ``n_bins`` is
``calc_bin_index(bin_size, end) - calc_bin_index(bin_size, start) + 1``
and blocks of adjacent regions split on a bin edge stack into the
block of their union. Bins are not clipped and no bounds are
returned: an edge row is a whole grid bin, holding only the records
whose position lies in ``[start, end]``.
**The empty-bin rule.** Null values are skipped, as every aggregator
skips them, and a cell whose aggregator received no non-null value
holds :data:`EMPTY_BIN_VALUES` for that aggregator: ``0`` for
``count`` and ``sum``, ``NaN`` for ``mean``, ``max``, ``min``,
``median`` and ``product``. Those seven are the whole set this fold
accepts; any other aggregator cannot answer a float and is refused
before the first record is read, as is a ``bin_size`` below 1.
**Preconditions.** ``records`` are sorted by position, and every
position lies in ``[start, end]``. The fold holds only the current
bin's aggregators, so it consumes the stream once, in order; a record
that breaks either precondition is refused when it is reached. Values
are whatever the aggregators accept -- numbers, or anything at all for
``count``.
"""
refused = [
name for name in aggregators if name not in EMPTY_BIN_VALUES]
if refused:
raise ValueError(
f"a binned fold cannot aggregate with {', '.join(refused)}: "
f"only {', '.join(EMPTY_BIN_VALUES)} answer a float per bin")
if bin_size < 1:
raise ValueError(
f"bins of size {bin_size} asked for; a bin holds at least one "
f"position")
first_bin = calc_bin_index(bin_size, start)
n_bins = calc_bin_index(bin_size, end) - first_bin + 1
block = np.empty((n_bins, len(aggregators)), dtype=np.float64)
block[:] = [EMPTY_BIN_VALUES[name] for name in aggregators]
def flush(row: int, accumulators: list[Aggregator]) -> None:
for column, accumulator in enumerate(accumulators):
final = accumulator.get_final()
if final is None:
continue
# An exact int (``sum`` of an int score) past the float range
# saturates, as a float result already does.
block[row, column] = to_float64(final)
row = -1
previous = start
accumulators: list[Aggregator] = []
for position, column, value in records:
if not start <= position <= end:
raise ValueError(
f"a record at {position} is outside the binned region "
f"[{start}, {end}]")
if position < previous:
raise ValueError(
f"a record at {position} arrived after one at {previous}; "
f"a binned fold needs its records sorted by position")
previous = position
record_row = calc_bin_index(bin_size, position) - first_bin
if record_row != row:
if row >= 0:
flush(row, accumulators)
row = record_row
accumulators = [Aggregator.build(name) for name in aggregators]
accumulators[column].add(value)
if row >= 0:
flush(row, accumulators)
return block