"""Handling of genomic scores statistics.
Currently we support only genomic scores histograms.
"""
from __future__ import annotations
import importlib.util
import json
import pathlib
import sys
from collections import Counter
from dataclasses import dataclass
from typing import IO, Any
import numpy as np
import yaml
from matplotlib import ticker
from gain import logging
from gain.genomic_resources.repository import GenomicResource
from gain.genomic_resources.resource_errors import HistogramError
from gain.genomic_resources.statistics.base_statistic import (
NON_NUMERIC_ERRORS,
PYTHON_NUMBER_TYPES,
Statistic,
as_python_number,
non_numeric_error,
)
from gain.genomic_resources.statistics.chart_style import (
CHART_FIGSIZE,
CHART_LABEL_FONT_SIZE,
)
from gain.genomic_resources.statistics.min_max import MinMaxValue
logger = logging.getLogger(__name__)
#: Which score value types a NUMBER histogram can accumulate, one value at a
#: time. ``bool`` is in: ``numpy`` folds it as 0/1, and a two-bin histogram
#: over a flag is meaningful. ``str`` is not, and a resource pairing the two
#: aborted its entire statistics build in ``np.isnan`` (gain#1285); since
#: gain#1336 that pairing is refused when the score is CONSTRUCTED, by
#: :func:`~gain.genomic_resources.score_resource.refuse_unfoldable_histograms`,
#: which is the only thing that reads this set.
#:
#: It lives here, with the config whose acceptance it describes, rather than
#: in the statistics scan that used to own it -- the refusal moved to the
#: score layer, and the score layer must not import the scan to ask what a
#: histogram can fold.
#:
#: Deliberately WIDER than the scan's ``_BULK_HISTOGRAM_VALUE_TYPES``, and
#: the two are different questions rather than one rule stated twice: this
#: asks what a histogram can fold value by value, that asks what it can fold
#: a whole column of (a bulk read yields a number histogram's column as
#: ``float64``, which a ``bool`` score's column is not). Merging them would
#: widen the vectorized path to a type it cannot read.
NUMBER_HISTOGRAM_VALUE_TYPES = ("float", "int", "bool")
[docs]
@dataclass
class NumberHistogramConfig:
"""Configuration class for number histograms."""
view_range: tuple[float | None, float | None]
number_of_bins: int = 100
x_log_scale: bool = False
y_log_scale: bool = False
x_min_log: float | None = None
plot_function: str | None = None
[docs]
def has_view_range(self) -> bool:
"""Whether both ends of the view range are pinned.
A histogram can only be built once its bin edges are known, so a
config whose range is half-open still needs the score's min/max
before it can be used.
"""
return self.view_range[0] is not None and \
self.view_range[1] is not None
[docs]
def to_dict(self) -> dict[str, Any]:
"""Transform number histogram config to dict."""
result = {
"type": "number",
"view_range": {
"min": self.view_range[0],
"max": self.view_range[1],
},
"number_of_bins": self.number_of_bins,
"x_log_scale": self.x_log_scale,
"y_log_scale": self.y_log_scale,
"x_min_log": self.x_min_log,
}
if self.plot_function is not None:
result["plot_function"] = self.plot_function
return result
[docs]
@staticmethod
def from_dict(parsed: dict[str, Any]) -> NumberHistogramConfig:
"""Build a number histogram config from a parsed yaml file."""
hist_type = parsed.get("type")
if hist_type != "number":
logger.error(
"Invalid configuration type (%s)"
" for number histogram!\n%s",
hist_type, parsed,
)
raise TypeError(
"Invalid configuration for number histogram!\n"
f"{parsed}",
)
yaml_range = parsed.get("view_range", {})
x_min = yaml_range.get("min", None)
x_max = yaml_range.get("max", None)
view_range = (x_min, x_max)
number_of_bins = parsed.get("number_of_bins", 100)
x_log_scale = parsed.get("x_log_scale", False)
y_log_scale = parsed.get("y_log_scale", False)
x_min_log = parsed.get("x_min_log")
plot_function = parsed.get("plot_function")
return NumberHistogramConfig(
view_range, number_of_bins,
x_log_scale, y_log_scale,
x_min_log,
plot_function=plot_function,
)
[docs]
@staticmethod
def default_config(
min_max: MinMaxValue | None,
) -> NumberHistogramConfig:
"""Build a number histogram config from a parsed yaml file."""
if min_max is None:
view_range: tuple[float | None, float | None] = (None, None)
elif min_max.min == min_max.max:
view_range = (min_max.min, min_max.min + 1.0)
else:
view_range = (min_max.min, min_max.max)
number_of_bins = 100
x_log_scale = False
y_log_scale = False
return NumberHistogramConfig(
view_range, number_of_bins, x_log_scale, y_log_scale)
DEFAULT_DISPLAYED_VALUES_COUNT = 20
[docs]
@dataclass
class CategoricalHistogramConfig:
"""Configuration class for categorical histograms."""
displayed_values_count: int | None = 20
displayed_values_percent: float | None = None
value_order: list[str | int] | None = None
y_log_scale: bool = False
label_rotation: int = 0
plot_function: str | None = None
enforce_type: bool = True
natural_order: bool = False
allow_only_whole_values_y: bool = False
[docs]
def to_dict(self) -> dict[str, Any]:
"""Transform categorical histogram config to dict."""
result: dict[str, Any] = {
"type": "categorical",
"value_order": self.value_order,
"y_log_scale": self.y_log_scale,
"label_rotation": self.label_rotation,
}
if self.displayed_values_count != DEFAULT_DISPLAYED_VALUES_COUNT:
result["displayed_values_count"] = self.displayed_values_count
if self.displayed_values_percent is not None:
result["displayed_values_percent"] = \
self.displayed_values_percent
if self.plot_function is not None:
result["plot_function"] = self.plot_function
return result
[docs]
@staticmethod
def default_config() -> CategoricalHistogramConfig:
"""The config for a score that declares no categorical ``histogram``.
``enforce_type=False`` here does NOT relax value-type checking --
:meth:`CategoricalHistogram.add_value` refuses a non-``str``/``int``
value either way. What the flag gates is
:attr:`CategoricalHistogram.UNIQUE_VALUES_LIMIT`, and it reads
backwards: it is the *default* config that enforces the limit, so a
score nobody declared categorical cannot silently histogram
thousands of distinct values. A curator who writes an explicit
``categorical`` block has asserted the score really is categorical,
and that config carries no limit.
"""
return CategoricalHistogramConfig(enforce_type=False)
[docs]
@staticmethod
def from_dict(parsed: dict[str, Any]) -> CategoricalHistogramConfig:
"""Create categorical histogram config from configuratin dict."""
hist_type = parsed.get("type")
if hist_type != "categorical":
raise TypeError(
"Invalid configuration type for categorical histogram!\n"
f"{parsed}",
)
displayed_values_count = parsed.get(
"displayed_values_count")
displayed_values_percent = parsed.get(
"displayed_values_percent")
if displayed_values_count is not None \
and displayed_values_percent is not None:
raise ValueError(
"Invalid configuration for categorical histogram: "
"displayed_values_count and displayed_values_percent "
"cannot be both set\n"
f"{parsed}",
)
if displayed_values_percent is None and displayed_values_count is None:
displayed_values_count = DEFAULT_DISPLAYED_VALUES_COUNT
value_order = parsed.get("value_order", [])
y_log_scale = parsed.get("y_log_scale", False)
plot_function = parsed.get("plot_function")
label_rotation = parsed.get("label_rotation", 0)
return CategoricalHistogramConfig(
displayed_values_count=displayed_values_count,
displayed_values_percent=displayed_values_percent,
value_order=value_order,
y_log_scale=y_log_scale,
plot_function=plot_function,
label_rotation=label_rotation,
enforce_type=True,
)
[docs]
@dataclass
class NullHistogramConfig:
"""Configuration class for null histograms."""
reason: str
[docs]
def to_dict(self) -> dict[str, Any]:
"""Render this config as the mapping ``from_dict`` reads back."""
return {
"type": "null",
"reason": self.reason,
}
[docs]
@staticmethod
def default_config() -> NullHistogramConfig:
"""A null config for a caller with no reason of its own to give.
``reason`` is required rather than optional, so that every null
histogram on a summary page can say why it is null; this is the
placeholder for the few call sites that genuinely have nothing to
add.
"""
return NullHistogramConfig("Unspecified reason")
[docs]
@staticmethod
def from_dict(parsed: dict[str, Any]) -> NullHistogramConfig:
"""Create Null histogram from configuration dict."""
hist_type = parsed.get("type")
if hist_type != "null":
raise TypeError(
"Invalid configuration type for null histogram!\n"
f"{parsed}",
)
reason = parsed.get("reason", "Unspecified reason")
return NullHistogramConfig(
reason=reason,
)
[docs]
class NumberHistogram(Statistic):
"""Class to represent a histogram."""
type = "number_histogram"
def __init__(
self, config: NumberHistogramConfig,
bins: np.ndarray | None = None,
bars: np.ndarray | None = None):
super().__init__("histogram", "Collects values for histogram.")
logger.debug("number histogram config: %s", config)
assert isinstance(config, NumberHistogramConfig)
self.config = config
self.out_of_range_values: list[float] = []
self.out_of_range_bins: list[int] = [0, 0]
self.min_value: float = np.nan
self.max_value: float = np.nan
self.choose_bin_index = self.choose_bin_lin
if self.config.x_log_scale:
self.choose_bin_index = self.choose_bin_log
if self.config.x_log_scale and self.config.x_min_log is None:
raise ValueError(
"Invalid histogram configuration, missing x_min_log",
)
if self.config.view_range[0] is None or \
self.config.view_range[1] is None or \
np.isnan(self.config.view_range[0]) or \
np.isnan(self.config.view_range[1]):
logger.error(
"unexpected min/max value: [%s, %s]",
self.config.view_range[0], self.config.view_range[1])
raise ValueError(
"unexpected min/max value:"
f"[{self.config.view_range[0]}, "
f"{self.config.view_range[1]}]")
self.view_range: tuple[float, float] = (
self.config.view_range[0], self.config.view_range[1])
if bins is not None and bars is not None:
self.bins = bins
self.bars = bars
elif bins is None and bars is None:
if self.config.x_log_scale:
assert self.config.x_min_log is not None
self.bins = np.array([
self.config.view_range[0],
* np.logspace(
np.log10(self.config.x_min_log),
np.log10(self.config.view_range[1]),
self.config.number_of_bins,
)])
self._rstep = (self.config.number_of_bins - 1) / \
(np.log10(self.view_max())
- np.log10(self.config.x_min_log))
else:
self.bins = np.linspace(
self.config.view_range[0],
self.config.view_range[1],
self.config.number_of_bins + 1,
)
if (self.view_max() - self.view_min()) <= 0:
self._rstep = 0
else:
self._rstep = self.config.number_of_bins / \
(self.view_max() - self.view_min())
self.bars = np.zeros(self.config.number_of_bins, dtype=np.int64)
assert not np.any(np.isnan(self.bins)), ("nan bins", self.config)
elif self.bins is None or self.bars is None:
raise ValueError(
"Cannot instantiate histogram with only bins or only bars!",
)
[docs]
def view_min(self) -> float:
"""The low edge of the first bin.
This is the histogram's *view* range, not the score's observed
minimum: values below it are counted in ``out_of_range_bins``
rather than binned.
"""
return self.view_range[0]
[docs]
def view_max(self) -> float:
"""The high edge of the last bin.
The counterpart of :meth:`view_min`; values above it are counted
as out of range.
"""
return self.view_range[1]
[docs]
def merge(self, other: Statistic) -> None:
"""Merge two histograms."""
assert isinstance(other, NumberHistogram)
assert self.bins is not None
assert self.bars is not None
assert other.bins is not None
assert other.bars is not None
assert np.allclose(self.bins, other.bins, rtol=1e-5), \
(self.bins, other.bins)
self.bars += other.bars
self.out_of_range_bins[0] += other.out_of_range_bins[0]
self.out_of_range_bins[1] += other.out_of_range_bins[1]
if np.isnan(self.min_value):
self.min_value = min(other.min_value, self.min_value)
else:
self.min_value = min(self.min_value, other.min_value)
if np.isnan(self.max_value):
self.max_value = max(other.max_value, self.max_value)
else:
self.max_value = max(self.max_value, other.max_value)
[docs]
def values_domain(self) -> str:
"""The observed value range, rendered for the summary page.
Unlike :meth:`view_min` / :meth:`view_max` this reports the values
actually seen, so a score whose configured view range is wider than
its data still shows the narrower true extent.
"""
return f"[{self.min_value:0.3g}, {self.max_value:0.3g}]"
[docs]
def add_value(
self, value: float | np.generic | None, count: int = 1,
) -> None:
"""Add value to the histogram.
``np.generic`` is in the signature because a numpy scalar is a real
caller's value, not a curiosity (gain#1338).
"""
# ``np.isnan`` is what refuses a value it cannot read as a number,
# and it raises BEFORE the allow-list below can be reached -- which
# is why that refusal, written for exactly the ``str`` case, could
# never fire for it, and a nullified score's reason read ``ufunc
# 'isnan' not supported`` (gain#1312). Re-wording numpy's complaint
# is the whole fix; the skip itself is untouched, so ``None`` is
# still an NA cell rather than a contract breach.
#
# Catching beats pre-checking the type here, and measurably: this
# runs per value of every record, and an ``isinstance`` ahead of the
# skip costs it ~24% (``str | bytes`` builds a union object on every
# call), where a try/except that does not fire costs ~3%.
try:
if value is None or np.isnan(value):
return
except NON_NUMERIC_ERRORS as err:
raise non_numeric_error(value, "number histogram") from err
# Reached only by values ``np.isnan`` accepted -- text and
# ``Decimal`` never get this far, the skip above refuses those. A
# Python number folds as-is; everything else goes through the rule
# the min/max twin shares, which says why a numpy scalar folds as
# ``item()`` and why a complex or an array does not (gain#1338,
# gain#1358). The isinstance stays inline because this runs per
# value of every record and a Python float is what the scan folds.
if not isinstance(value, PYTHON_NUMBER_TYPES):
value = as_python_number(value, "number histogram")
self.min_value = min(value, self.min_value)
self.max_value = max(value, self.max_value)
index = self.choose_bin_index(value)
if index < 0:
logger.warning(
"out of range %s value %s", self.view_range, value)
tindex = index + 2
self.out_of_range_bins[tindex] += count
return
self.bars[index] += count
[docs]
def add_batch(
self, values: np.ndarray, weights: np.ndarray,
) -> None:
"""Add a batch of ``(value, weight)`` pairs, vectorized.
Bit-for-bit equivalent to calling :meth:`add_value` over each pair in
order: the same bin selection for both x-scales (the truncation of
``choose_bin_lin`` / ``choose_bin_log`` and the clamp to
``number_of_bins - 1``), the same below/above ``out_of_range_bins``
split, the same ``min_value``/``max_value`` tracking, and the same
nan-skip. This is the hot path for the statistics scan, where a
per-value Python call dominates the cost.
The equivalence covers dtype as well as arithmetic: the ``float64``
coercion below is what ``add_value`` reproduces by normalizing a
numpy scalar, so a narrow column -- ``float32``, ``float16``,
``bool`` -- folds the same through either arm (gain#1338).
Where it does not hold: this arm also folds values ``add_value``
refuses outright, because ``asarray`` converts them silently -- a
``complex`` loses its imaginary part, and an object array of
``Decimal`` converts. Nothing in a scan produces either; they are
reachable only by calling this directly.
Both x-scales are vectorized. The log one is bit-exact for the same
reason the linear one is and one more: ``np.log10`` returns the same
float for a scalar as for that scalar inside an array, at every array
width numpy dispatches differently on -- checked over many decades by
test_add_batch_matches_add_value_loop_log_fuzz, which varies batch size
precisely to cover numpy's scalar-loop and SIMD kernels.
"""
values = np.asarray(values, dtype=np.float64)
weights = np.asarray(weights, dtype=np.int64)
finite = ~np.isnan(values)
if not finite.any():
return
values = values[finite]
weights = weights[finite]
# ``add_value`` seeds min/max at nan and folds each value in with
# ``min(value, self.min_value)``; over a set that is just the extremum.
batch_min = float(values.min())
batch_max = float(values.max())
self.min_value = batch_min if np.isnan(self.min_value) \
else min(self.min_value, batch_min)
self.max_value = batch_max if np.isnan(self.max_value) \
else max(self.max_value, batch_max)
below = values < self.view_min()
above = values > self.view_max()
self.out_of_range_bins[0] += int(weights[below].sum())
self.out_of_range_bins[1] += int(weights[above].sum())
in_range = ~(below | above)
idx = self._bin_indices(values[in_range])
# Clamp to the last bin, exactly as choose_bin_lin/choose_bin_log do.
np.minimum(idx, self.config.number_of_bins - 1, out=idx)
np.add.at(self.bars, idx, weights[in_range])
def _bin_indices(self, values: np.ndarray) -> np.ndarray:
"""Vectorized, unclamped bin index for values already known in range.
The array counterpart of :meth:`choose_bin_lin` /
:meth:`choose_bin_log`, minus the below/above tests its caller has
already applied. Both scalar forms truncate with ``int(...)``, which
rounds toward zero; the offsets here are non-negative, so that is
floor, which ``astype(np.int64)`` reproduces.
"""
if not self.config.x_log_scale:
linear: np.ndarray = (
(values - self.view_min()) * self._rstep).astype(np.int64)
return linear
assert self.config.x_min_log is not None
x_min_log = self.config.x_min_log
# choose_bin_log puts everything under x_min_log in bin 0 -- the bin
# spanning view_min..x_min_log, which is why the log formula's index
# starts at 1 -- and log-bins the rest. Zeros give the first case for
# free, so only the second is computed.
indices = np.zeros(values.shape, dtype=np.int64)
logged = values >= x_min_log
indices[logged] = ((
np.log10(values[logged]) - np.log10(x_min_log))
* self._rstep).astype(np.int64) + 1
return indices
[docs]
def choose_bin_lin(self, value: float) -> int:
"""Compute bin index for a passed value for linear x-scale."""
if value < self.view_min():
return -2
if value > self.view_max():
return -1
index = int((value - self.view_min()) * self._rstep)
return min(index, self.config.number_of_bins - 1)
[docs]
def choose_bin_log(self, value: float) -> int:
"""Compute bin index for a passed value for log x-scale."""
assert self.config.x_log_scale
assert self.config.x_min_log is not None
if value < self.view_min():
return -2
if value > self.view_max():
return -1
if value < self.config.x_min_log:
return 0
index = int(
(np.log10(value) - np.log10(self.config.x_min_log))
* self._rstep) + 1
return min(index, self.config.number_of_bins - 1)
[docs]
def to_dict(self) -> dict[str, Any]:
"""Render this histogram as the mapping :meth:`from_dict` reads back.
Bin edges and bar counts are carried as plain lists rather than
arrays, so the result is JSON-serialisable as it stands.
"""
return {
"config": self.config.to_dict(),
"bins": self.bins.tolist(),
"bars": self.bars.tolist(),
"out_of_range_bins": self.out_of_range_bins,
"min_value": float(self.min_value),
"max_value": float(self.max_value),
}
[docs]
def serialize(self) -> str:
"""Render this histogram as the JSON stored in the resource."""
return json.dumps(self.to_dict(), indent=2)
[docs]
def plot(
self,
outfile: IO,
score_id: str,
y_axis_label: str | None = None,
small_values_description: str | None = None,
large_values_description: str | None = None,
) -> None:
"""Plot histogram and save it into outfile."""
# pylint: disable=import-outside-toplevel
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
width = self.bins[1:] - self.bins[:-1]
fig, ax = plt.subplots(figsize=CHART_FIGSIZE, tight_layout=True)
ax.bar(
x=self.bins[:-1], height=self.bars,
log=self.config.y_log_scale,
width=width,
align="edge")
if self.config.x_log_scale:
ax.set_xscale("log")
if (
small_values_description is not None
and large_values_description is not None
):
sec = ax.secondary_xaxis(location=0)
if self.config.x_log_scale:
left_location = self.bins[1]
else:
left_location = self.bins[0]
right_location = self.bins[-1]
sec.set_ticks(
[
left_location,
right_location,
],
labels=[
f"\n{small_values_description}",
f"\n{large_values_description}",
],
wrap=True,
color="gray",
style="italic",
fontsize=CHART_LABEL_FONT_SIZE,
)
ax.set_xlabel(f"\n{score_id}", fontsize=CHART_LABEL_FONT_SIZE)
ax.set_ylabel(
"count" if y_axis_label is None else y_axis_label,
fontsize=CHART_LABEL_FONT_SIZE)
ax.grid(axis="y")
ax.grid(axis="x")
fig.savefig(outfile)
plt.close(fig)
[docs]
@staticmethod
def from_dict(data: dict[str, Any]) -> NumberHistogram:
"""Build a number histogram from a dict."""
config = NumberHistogramConfig.from_dict(data["config"])
hist = NumberHistogram(
config,
bins=np.array(data.get("bins")),
bars=np.array(data.get("bars")),
)
hist.min_value = data.get("min_value", np.nan)
hist.max_value = data.get("max_value", np.nan)
hist.out_of_range_bins = data.get("out_of_range_bins", [0, 0])
return hist
[docs]
@staticmethod
def deserialize(content: str) -> NumberHistogram:
"""Rebuild a number histogram from :meth:`serialize` output."""
data = json.loads(content)
return NumberHistogram.from_dict(data)
[docs]
class HistogramStatisticMixin:
"""Mixin for creating statistics classes with histograms."""
[docs]
@staticmethod
def get_histogram_file(score_id: str) -> str:
return f"histogram_{score_id}.yaml"
[docs]
@staticmethod
def get_histogram_image_file(score_id: str) -> str:
return f"histogram_{score_id}.png"
[docs]
class NullHistogram(Statistic):
"""Class for annulled histograms."""
type = "null_histogram"
def __init__(self, config: NullHistogramConfig | None) -> None:
super().__init__(
"null_histogram", "Used for invalid/annulled histograms",
)
if config is None:
config = NullHistogramConfig.default_config()
self.reason = config.reason
[docs]
def add_value(
self, value: Any, count: int = 1, # ruff: ignore[unused-method-argument]
) -> None:
"""Discard the value.
A null histogram counts nothing by design, so that a statistics
build can feed every score the same way without first asking
whether this one has a histogram.
"""
# pylint: disable=unused-argument
return
[docs]
def merge(self, other: Any) -> None: # ruff: ignore[unused-method-argument]
"""Do nothing: there are no counts to fold together.
Merging is a no-op rather than an error so that a parallel build
can reduce its partial results uniformly, null histograms included.
"""
return
[docs]
def to_dict(self) -> dict[str, Any]:
"""Render this histogram as the mapping :meth:`from_dict` reads back.
Only the config survives a round trip, because the reason is the
whole of a null histogram's content.
"""
return {
"config": {
"type": "null",
"reason": self.reason,
},
}
[docs]
def values_domain(self) -> str:
"""Report that there is no domain, in the other kinds' place."""
return "NO DOMAIN"
# pylint: disable=unused-argument
[docs]
def plot(self, _outfile: IO, _score_id: str) -> None:
"""Draw nothing, leaving ``outfile`` untouched.
The caller is expected to render ``reason`` in place of the image,
rather than to link an image file this never wrote.
"""
return
[docs]
def serialize(self) -> str:
"""Render this histogram as the JSON stored in the resource."""
return json.dumps(self.to_dict(), indent=2)
[docs]
@staticmethod
def from_dict(data: dict[str, Any]) -> NullHistogram:
"""Build a null histogram from a dict."""
config = data["config"]
hist_type = config.get("type")
if hist_type != "null":
raise TypeError(
f"Invalid configuration type for null histogram!\n{data}",
)
reason = config.get("reason", "")
return NullHistogram(NullHistogramConfig(reason=reason))
[docs]
@staticmethod
def deserialize(content: str) -> NullHistogram:
"""Rebuild a null histogram from :meth:`serialize` output."""
data = json.loads(content)
return NullHistogram.from_dict(data)
[docs]
class CategoricalHistogram(Statistic):
"""Class for categorical data histograms."""
type = "categorical_histogram"
UNIQUE_VALUES_LIMIT = 100
# pylint: disable=too-few-public-methods
def __init__(
self,
config: CategoricalHistogramConfig,
counter: dict[str | int, int] | None = None,
*,
truncated: bool = False,
unique_values: int | None = None,
total_count: int | None = None,
):
super().__init__(
"categorical_histogram",
"Collects values for categorical histogram.",
)
self.config = config
self.enforce_type = config.enforce_type
if counter is not None:
self._counter = Counter(counter)
else:
self._counter = Counter()
self.truncated = truncated
self._unique_values = unique_values
self._total_count = total_count
self.y_log_scale = config.y_log_scale
[docs]
def add_value(
self, value: str | int | None,
count: int = 1,
) -> None:
"""Add a value to the categorical histogram.
Returns true if successfully added and false if failed.
Will fail if too many values are accumulated.
"""
if value is None:
return
if not isinstance(value, str | int):
raise TypeError(
"Only string or int values can be added categorical "
f"histogram; bad <{value}>",
)
self._counter[value] += count
if not self.enforce_type and \
len(self._counter) > CategoricalHistogram.UNIQUE_VALUES_LIMIT:
raise HistogramError(
f"Too many unique values {len(self._counter)} "
f"for categorical histogram.",
)
[docs]
def add_batch(
self, values: np.ndarray, weights: np.ndarray,
) -> None:
"""Add a batch of ``(value, weight)`` pairs, vectorized.
Equivalent to calling :meth:`add_value` over each pair in order: the
same ``None`` skip, the same per-value counts, the same ``TypeError``
naming the first value that is neither ``str`` nor ``int``, and the
same ``HistogramError`` when the batch takes the histogram past
``UNIQUE_VALUES_LIMIT``. This is the hot path for the statistics
scan, where a per-value Python call dominates the cost.
The limit's message reports ``UNIQUE_VALUES_LIMIT + 1`` because that
is the count :meth:`add_value` always raises at: it tests after every
single add, so the first add that exceeds the limit is the one that
raises, and the histogram holds exactly one value too many. Which
values the counter holds when it raises is not otherwise observable --
a histogram that raises is replaced by a ``NullHistogram`` carrying
the message.
A batch containing BOTH a value of an unusable type and enough new
values to trip the limit reports the type failure, wherever the two
sit relative to each other; the per-record path reports whichever
comes first. The distinction is unreachable through the scan, which
batches a ``str`` score's column -- every cell of which is a ``str``
or the ``None`` this skips.
"""
cells = np.asarray(values, dtype=object).tolist()
if weights.size and bool(np.all(weights == 1)):
# The C-level count, for the kinds that weigh every record 1.
batch = Counter(cells)
batch.pop(None, None)
else:
batch = Counter()
for cell, weight in zip(cells, weights.tolist(), strict=True):
if cell is None:
continue
batch[cell] += weight
for cell in batch:
if not isinstance(cell, str | int):
raise TypeError(
"Only string or int values can be added categorical "
f"histogram; bad <{cell}>",
)
self._counter.update(batch)
if not self.enforce_type and \
len(self._counter) > CategoricalHistogram.UNIQUE_VALUES_LIMIT:
raise HistogramError(
f"Too many unique values "
f"{CategoricalHistogram.UNIQUE_VALUES_LIMIT + 1} "
f"for categorical histogram.",
)
[docs]
def merge(self, other: Statistic) -> None:
"""Merge with other histogram."""
assert isinstance(other, CategoricalHistogram)
assert self.config == other.config
if self.truncated or other.truncated:
raise HistogramError(
"Can not merge a truncated categorical histogram sidecar; "
"merge needs the full histogram values.",
)
# pylint: disable=protected-access
self._counter += other._counter # ruff: ignore[private-member-access]
if not self.enforce_type and \
len(self._counter) > CategoricalHistogram.UNIQUE_VALUES_LIMIT:
raise HistogramError(
f"Can not merge categorical histograms; "
f"too many unique values {len(self._counter)}")
@property
def raw_values(self) -> dict[str | int, int]:
"""Every counted value with its count, unordered.
This is the histogram's own content, as distinct from
:attr:`display_values`, which is the ordered subset the summary
page draws. A custom ``plot_function`` receives the histogram
itself and so can read either.
Note that on a histogram loaded back from a truncated sidecar
(:meth:`serialize_truncated`) the counter holds only the values
that sidecar carried, so this returns the truncated set --
:attr:`unique_values` and :attr:`total_count` are the ones that
still describe the full data.
"""
return dict(self._counter)
@property
def unique_values(self) -> int:
"""Number of distinct values, across truncation.
On a truncated instance this is the distinct-value count of the full
histogram the sidecar was derived from, not the length of the top-N
counter it carries.
"""
if self._unique_values is not None:
return self._unique_values
return len(self._counter)
@property
def total_count(self) -> int:
"""Sum of all value counts, across truncation.
On a truncated instance this is the total of the full histogram the
sidecar was derived from, not the sum of the top-N counter it carries.
"""
if self._total_count is not None:
return self._total_count
return sum(self._counter.values())
@property
def display_values(self) -> dict[str | int, int]:
"""Return categorical histogram display values in order.
A truncated instance carries exactly the values its config selected
for display at serialization time, so they are returned verbatim --
re-running the selection against the truncated counter would compute
percentages and orderings against the wrong totals.
"""
if self.truncated:
return dict(self._counter)
values = {}
if self.config.value_order:
for key in self.config.value_order:
values[key] = self._counter[key]
if len(values) < len(self._counter):
raise ValueError(
"misconfigured categorical histogram value_order",
f"{self.config.value_order} < {self._counter.keys()}")
return values
if self.config.displayed_values_percent is not None:
total = sum(self._counter.values())
displayed = 0
other = 0
displayed_percent = self.config.displayed_values_percent
for key, count in self._counter.most_common():
if 100.0 * displayed / total < displayed_percent:
values[key] = count
displayed += count
else:
other += count
if other > 0:
values["Other"] = other
return values
ordering = None
if self.config.natural_order:
ordering = sorted(
self._counter, key=lambda x: (isinstance(x, str), x))
else:
ordering = [key for key, _ in self._counter.most_common()]
for key in ordering[:self.config.displayed_values_count]:
values[key] = self._counter[key]
if self.config.displayed_values_count is not None and \
len(self._counter) > self.config.displayed_values_count:
other = 0
for key in ordering[self.config.displayed_values_count:]:
other += self._counter[key]
if other > 0:
values["Other"] = other
return values
[docs]
def values_domain(self) -> str:
"""Render the displayed values, noting truncation when present."""
domain = ", ".join(str(k) for k in self.display_values)
if self.truncated and len(self._counter) < self.unique_values:
return (
f"{domain} "
f"(top {len(self._counter)} of {self.unique_values} values)"
)
return domain
[docs]
def to_dict(self) -> dict[str, Any]:
"""Render this histogram as the mapping :meth:`from_dict` reads back.
Carries *every* counted value, not the displayed subset; the
truncated companion form is :meth:`serialize_truncated`.
"""
return {
"config": self.config.to_dict(),
"values": dict(self._counter),
}
[docs]
def serialize(self) -> str:
"""Render the full histogram as the JSON stored in the resource."""
return json.dumps(self.to_dict(), indent=2)
[docs]
def serialize_truncated(self) -> str:
"""Serialize the truncated sidecar form of this histogram.
The sidecar carries the config, the values this histogram's config
selects for display, and the ``unique_values``/``total_count``
totals of the full histogram, marked with ``"truncated": true``. It
is the small, always-readable companion of a full histogram whose
values file may be absent from a checkout (DVC-tracked, not pulled).
The carried values follow the display selection so the sidecar
renders what the full histogram would: every ``value_order`` key
with its real count for an ordered config, the values covering
``displayed_values_percent`` for a percent config, and the
``displayed_values_count`` most common values otherwise.
"""
values: dict[str | int, int]
if self.config.value_order:
# A JSON round-trip stringifies the counter's keys while the
# config may order the values by their original (e.g. int)
# keys, so a histogram loaded back from disk falls back to
# the stringified key -- it must serialize the same counts a
# live one does.
values = {
key: self._counter.get(key, self._counter.get(str(key), 0))
for key in self.config.value_order}
elif self.config.displayed_values_percent is not None:
values = dict(self._counter.most_common(len(self.display_values)))
else:
top = self.config.displayed_values_count \
or DEFAULT_DISPLAYED_VALUES_COUNT
values = dict(self._counter.most_common(top))
return json.dumps({
"config": self.config.to_dict(),
"values": values,
"truncated": True,
"unique_values": self.unique_values,
"total_count": self.total_count,
}, indent=2)
[docs]
@staticmethod
def from_dict(data: dict[str, Any]) -> CategoricalHistogram:
"""Build a categorical histogram from a dict.
Reads both forms :meth:`to_dict` and :meth:`serialize_truncated`
produce: the truncation flag and the two totals are absent from a
full histogram and default accordingly.
"""
config = CategoricalHistogramConfig.from_dict(data["config"])
return CategoricalHistogram(
config,
data.get("values"),
truncated=data.get("truncated", False),
unique_values=data.get("unique_values"),
total_count=data.get("total_count"),
)
[docs]
@staticmethod
def deserialize(content: str) -> CategoricalHistogram:
"""Rebuild a categorical histogram from :meth:`serialize` output."""
data = json.loads(content)
return CategoricalHistogram.from_dict(data)
[docs]
def plot(
self,
outfile: IO,
score_id: str,
*,
y_axis_label: str | None = None,
small_values_description: str | None = None,
large_values_description: str | None = None,
) -> None:
"""Plot histogram and save it into outfile."""
# pylint: disable=import-outside-toplevel
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
display_values = self.display_values
values = [str(k) for k in display_values]
counts = list(display_values.values())
fig, ax = plt.subplots(figsize=CHART_FIGSIZE, tight_layout=True)
ax.bar(
x=values,
height=counts,
tick_label=[str(v) for v in values],
log=self.config.y_log_scale,
align="center",
)
if len(values) == 1:
ax.set_xlim(-2.5, 2.5)
if self.config.allow_only_whole_values_y:
ax.yaxis.set_major_locator(ticker.MaxNLocator(integer=True))
if small_values_description is not None and \
large_values_description is not None:
sec = ax.secondary_xaxis(location=0)
sec.set_ticks(
[
0,
len(values) - 1,
],
labels=[
f"\n{small_values_description}",
f"\n{large_values_description}",
],
wrap=True,
color="gray",
style="italic",
fontsize=CHART_LABEL_FONT_SIZE,
)
ax.set_xlabel(f"\n{score_id}", fontsize=CHART_LABEL_FONT_SIZE)
ax.set_ylabel(
"count" if y_axis_label is None else y_axis_label,
fontsize=CHART_LABEL_FONT_SIZE)
label_angle = self.config.label_rotation % 360
if self.config.label_rotation < 0:
label_angle = 360 + label_angle
if label_angle != 0:
ax.set_xticklabels(
[str(v) for v in values],
rotation=self.config.label_rotation,
ha="center",
va="top",
rotation_mode="default",
fontsize=CHART_LABEL_FONT_SIZE,
)
fig.savefig(outfile)
plt.close(fig)
[docs]
def build_histogram_config(
config: dict[str, Any] | None) -> HistogramConfig | None:
"""Create histogram config form configuration dict."""
if config is None:
return None
if "histogram" in config:
hist_config = config["histogram"]
if "type" not in hist_config:
return NullHistogramConfig(
f"Missing histogram type in configuration {config}")
hist_type = hist_config["type"]
else:
return None
if hist_type == "number":
return NumberHistogramConfig.from_dict(hist_config)
if hist_type == "categorical":
return CategoricalHistogramConfig.from_dict(hist_config)
if hist_type == "null":
return NullHistogramConfig.from_dict(hist_config)
return NullHistogramConfig(f"Invalid histogram configuration {config}")
[docs]
def build_default_histogram_conf(
value_type: str, **kwargs: Any,
) -> NumberHistogramConfig | CategoricalHistogramConfig | NullHistogramConfig:
"""Build default histogram config for given value type."""
if value_type in ["int", "float"]:
min_max = kwargs.get("min_max")
return NumberHistogramConfig.default_config(min_max)
if value_type == "str":
return CategoricalHistogramConfig.default_config()
return NullHistogramConfig(
"No histogram configured and no default config available for type"
f"{value_type}",
)
[docs]
def build_empty_histogram(
config: HistogramConfig,
) -> NumberHistogram | CategoricalHistogram | NullHistogram:
"""Create an empty histogram from a deserialize histogram dictionary."""
# pylint: disable=broad-except
try:
if isinstance(config, NumberHistogramConfig):
return NumberHistogram(config)
if isinstance(config, CategoricalHistogramConfig):
return CategoricalHistogram(config)
if isinstance(config, NullHistogramConfig):
return NullHistogram(config)
return NullHistogram(NullHistogramConfig(
"Could not match histogram config type",
))
except BaseException as err:
logger.warning(
"Failed to create empty histogram from config", exc_info=True)
return NullHistogram(NullHistogramConfig(
f"Failed to create empty histogram from config: {err}"))
[docs]
def load_histogram(
resource: GenomicResource, filename: str,
) -> Histogram:
"""Load and return a histogram in a resource.
On an error or missing histogram, an appropriate NullHistogram is returned.
"""
try:
with resource.open_raw_file(filename) as infile:
content = infile.read()
except FileNotFoundError:
# Handled, not fatal: a null histogram is returned. So no
# traceback -- during a repair run this is an expected consequence
# of statistics that have not been built yet, and its traceback was
# the only one the user saw, pointing away from the real fault
# (gain#364). WARNING rather than DEBUG keeps a genuinely missing
# histogram on an otherwise healthy resource visible.
logger.warning(
"unable to load histogram file: %s; file not found", filename)
return NullHistogram(NullHistogramConfig(
"Histogram file not found.",
))
if filename.endswith(".yaml"):
hist_data = yaml.safe_load(content)
elif filename.endswith(".json"):
hist_data = json.loads(content)
else:
logger.error(
"Invalid histogram file format: %s", filename)
return NullHistogram(NullHistogramConfig(
"Invalid histogram file format.",
))
config = hist_data["config"]
hist_type = config["type"]
try:
if hist_type == "number":
return NumberHistogram.from_dict(hist_data)
if hist_type == "categorical":
return CategoricalHistogram.from_dict(hist_data)
if hist_type == "null":
return NullHistogram.from_dict(hist_data)
return NullHistogram(NullHistogramConfig("Invalid histogram type"))
except BaseException: # pylint: disable=broad-except
logger.exception(
"Failed to deserialize histogram from %s",
filename,
)
return NullHistogram(NullHistogramConfig(
"Failed to deserialize histogram.",
))
[docs]
def truncated_histogram_filename(histogram_filename: str) -> str:
"""Return the truncated-sidecar filename for a histogram filename.
``statistics/histogram_cell.json`` maps to
``statistics/truncated/histogram_cell.json``: a directory rather
than a name suffix, because the score-id part of the histogram
filename is arbitrary (a score may itself be named
``foo_truncated``) while no score id can contain ``/``.
"""
directory, sep, basename = histogram_filename.rpartition("/")
return f"{directory}{sep}truncated/{basename}"
[docs]
def save_histogram(
resource: GenomicResource,
filename: str,
histogram: Histogram,
) -> None:
"""Save histogram into a resource."""
if not filename.endswith(".json"):
logger.error(
"Invalid histogram file format: %s; "
"histograms are stored in JSON format only", filename)
raise ValueError(
f"Invalid histogram file format: <{filename}>; "
f"histograms are stored in JSON format only",
)
with resource.open_raw_file(filename, mode="wt") as outfile:
outfile.write(histogram.serialize())
HistogramConfig = \
NullHistogramConfig | CategoricalHistogramConfig | NumberHistogramConfig
Histogram = NullHistogram | CategoricalHistogram | NumberHistogram
def _import_from_string(module_name: str, source_code: str) -> Any:
spec = importlib.util.spec_from_loader(module_name, loader=None)
assert spec is not None
module = importlib.util.module_from_spec(spec)
exec(source_code, module.__dict__) # ruff: ignore[exec-builtin] pylint: disable=exec-used
sys.modules[spec.name] = module
return module
[docs]
def plot_histogram(
res: GenomicResource,
image_filename: str,
hist: Histogram,
score_id: str,
small_values_desc: str | None = None,
large_values_desc: str | None = None,
) -> None:
"""Plot histogram and save it into the resource."""
if isinstance(hist, NullHistogram):
return
if hist.config.plot_function is None:
with res.open_raw_file(image_filename, mode="wb") as outfile:
hist.plot(
outfile,
score_id,
small_values_description=small_values_desc,
large_values_description=large_values_desc,
)
return
plot_file, plot_function_name = hist.config.plot_function.split(":")
with res.open_raw_file(
plot_file,
mode="rt",
) as srcfile:
source_code = srcfile.read()
plot_module_name = str(
pathlib.Path(res.resource_id) /
pathlib.Path(plot_file).with_suffix("")).replace("/", ".")
plot_module = _import_from_string(
plot_module_name, source_code)
func = getattr(plot_module, plot_function_name)
with res.open_raw_file(
image_filename,
mode="wb",
) as outfile:
func(
outfile,
hist,
score_id,
small_values_desc,
large_values_desc,
)
return