from __future__ import annotations
import logging
import warnings
from typing import cast
import numpy as np
import pandas as pd
import scipy as sp
from tqdm import tqdm
try:
import plotly.graph_objects as go
except ImportError: # pragma: no cover - exercised via env-without-plotly
from process_improve._extras import _MissingExtra
go = _MissingExtra("plotly", "plotting") # type: ignore[assignment]
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.utils.validation import check_is_fitted
from ..multivariate.methods import PCA, MCUVScaler
from .alignment_helpers import backtrack_optimal_path, distance_matrix
from .data_input import check_valid_batch_dict, dict_to_wide, melted_to_dict
logger = logging.getLogger(__name__)
epsqrt = np.sqrt(np.finfo(float).eps)
def _varies_within_any_batch(batches: dict[str, pd.DataFrame], column: object) -> bool:
"""
Report whether ``column`` takes more than one value inside at least one batch.
A column that is flat through every batch carries no trajectory for the alignment to
work with. Every batch is checked rather than only the first, so a tag that happens to
be flat in the first batch but moves in a later one is kept.
"""
return any(_has_spread(batch[column].to_numpy()) for batch in batches.values() if column in batch)
def _has_spread(values: np.ndarray) -> bool:
"""Report whether a numeric array holds more than one distinct usable value."""
usable = values[np.isfinite(values)]
# min against max rather than counting distinct values: one pass, no hashing, and the
# question is only whether the column moves at all.
return bool(usable.size and usable.min() != usable.max())
def _resolve_columns_to_align(
batches: dict[str, pd.DataFrame],
columns_to_align: list | pd.Index | None,
caller: str,
) -> list | pd.Index:
"""
Validate the batch container and resolve which columns to operate on.
The scaling functions take batches as a dict keyed by batch identifier. A single
wide DataFrame holding every batch used to get past the column-resolution branch
and then fail several lines later inside the loop, because ``DataFrame.items()``
yields ``(column, Series)`` pairs rather than ``(batch_id, frame)``. The error
that surfaced named a ``Series`` the caller never created. Reject the
unsupported container here instead, and say how to convert it. (#560)
Parameters
----------
batches : dict[str, pd.DataFrame]
Batch data, in the standard format (keyed by batch identifier).
columns_to_align : list, pd.Index, or None
Checked when given; resolved from the numeric columns of the first batch when
``None``.
caller : str
Name of the calling function, used in the error messages.
Returns
-------
list or pd.Index
The columns to operate on.
Raises
------
TypeError
If ``batches`` is a DataFrame, or is not a dict.
ValueError
If ``batches`` is an empty dict, so there is no batch to take columns from; if a
batch holds no numeric column at all; or if an explicit ``columns_to_align``
names a column that is not numeric.
Notes
-----
Only numeric columns are aligned. Resolving from every column swept in whatever else
a batch frame carried, most often the batch identifier that
:func:`~process_improve.batch.data_input.melted_to_dict` leaves in place: constant
within a batch, so it contributed a zero range, and meaningless to scale. A column
that is not numeric is not a trajectory, and what to do with it (carry it, drop it,
encode it) is the caller's decision, not something to guess at here.
A numeric column that holds one value throughout every batch is dropped for the same
reason. That is what an identifier column looks like, and the identifier
:func:`~process_improve.batch.data_input.melted_to_dict` leaves in place is usually a
number, so a dtype test alone does not catch it. Including it is not harmless: on the
dryer data it takes almost no weight itself (0.000051) but moves the others
substantially, ``JacketTemperatureSP`` from 0.132 to 0.472, because it joins the
distance the alignment minimises and the normalisation that follows.
An explicit ``columns_to_align`` naming a non-numeric column raises rather than
quietly dropping it, since the caller asked for it by name. A constant column named
explicitly is left alone: constant over this particular set of batches does not mean
constant in general, and the caller may know better.
"""
if isinstance(batches, pd.DataFrame):
raise TypeError(
f"{caller} expects `batches` as a dict of per-batch DataFrames, keyed by batch "
f"identifier; got a single {type(batches).__name__} holding every batch, which is not "
"supported yet (tracked on #199). Split it per batch first, for example: "
"`dict(tuple(df.groupby(batch_col)))`."
)
if not isinstance(batches, dict):
raise TypeError(
f"{caller} expects `batches` as a dict of per-batch DataFrames, keyed by batch "
f"identifier; got {type(batches).__name__}."
)
if not batches:
raise ValueError(
f"{caller} cannot resolve `columns_to_align` from an empty `batches` dict; "
"pass `columns_to_align` explicitly, or supply at least one batch."
)
first_batch = batches[next(iter(batches))]
if columns_to_align is not None:
non_numeric = {
str(column): str(first_batch[column].dtype)
for column in columns_to_align
if column in first_batch and not pd.api.types.is_numeric_dtype(first_batch[column])
}
if non_numeric:
listed = ", ".join(f"{name!r} ({dtype})" for name, dtype in sorted(non_numeric.items()))
raise ValueError(
f"{caller} can only align numeric columns, but `columns_to_align` names: {listed}. "
"Drop them from `columns_to_align`; a non-numeric column is not a trajectory, and "
"carrying, encoding or discarding it is yours to decide."
)
return columns_to_align
numeric = first_batch.select_dtypes(include="number").columns
varying = [column for column in numeric if _varies_within_any_batch(batches, column)]
if not varying:
raise ValueError(
f"{caller} found no column to align in the first batch, whose columns are "
f"{list(first_batch.columns)}. A column must be numeric and must vary within at least "
"one batch. Pass `columns_to_align` explicitly, or supply batches holding at least one "
"numeric trajectory."
)
return varying
#: Quantile pair behind each ``settings["robust_range"]`` choice in :func:`determine_scaling`.
_ROBUST_RANGE_QUANTILES = {"q98-q02": (0.02, 0.98), "iqr": (0.25, 0.75)}
def _warn_about_collapsed_ranges(collapsed: dict[str, int], n_batches: int) -> None:
"""
Report tags whose range was zero and was replaced by 1.0, once for the whole call.
A zero range means the tag held a single value across the batch, so the substitution
leaves it unscaled. That is recoverable, but the caller should know: a constant tag
carries nothing for the alignment to work with, and an identifier column swept in by
the default ``columns_to_align`` shows up here too.
"""
if not collapsed:
return
listed = ", ".join(
f"{tag!r} ({count} of {n_batches} batches)"
for tag, count in sorted(collapsed.items(), key=lambda item: (-item[1], item[0]))
)
warnings.warn(
f"determine_scaling found tags with a zero range and substituted 1.0, leaving them "
f"unscaled: {listed}. A tag that holds one value across a batch carries no information "
f"for alignment; exclude it with `columns_to_align`.",
category=UserWarning,
stacklevel=3,
)
[docs]
def determine_scaling(
batches: dict[str, pd.DataFrame],
columns_to_align: list | pd.Index | None = None,
settings: dict | None = None,
) -> pd.DataFrame:
"""
Scales the batch data according to the variable ranges.
Parameters
----------
batches : dict[str, pd.DataFrame]
Batch data, in the standard format (keyed by batch identifier).
columns_to_align : list, optional
The column names (tags) to be scaled. If ``None``, the columns of the
first batch are used.
settings : dict, optional
Optional overrides:
``"robust"`` (bool, default ``True``)
Switches between a robust per-batch range and the raw ``max - min``.
``"robust_range"`` (str, default ``"q98-q02"``)
Which robust range to use when ``robust`` is True, either
``"q98-q02"`` or ``"iqr"`` (q75 - q25, the interquartile range that
:func:`~process_improve.batch.features.f_iqr` computes). Ignored when
``robust`` is False.
The two are not interchangeable, which is the answer to the question the
old TODO here posed. The IQR spans the middle half of a batch; q98 - q02
spans nearly all of it. On a Gaussian tag the second is about 3.05 times
the first, but a batch trajectory is not Gaussian, so the ratio varies by
tag: measured per tag on the bundled data it runs from 1.02 to 4.23
(dryer) and 1.21 to 2.95 (nylon). Switching therefore re-weights the tags
against each other rather than rescaling them together. The IQR also
collapses to zero more often, because a tag that holds one value for more
than half of a batch has no interquartile spread at all:
``DifferentialPressure`` collapses in 21 of the 71 dryer batches under the
IQR against 16 under q98 - q02. Prefer the IQR when the tags carry
excursions you want the scaling to ignore; the default otherwise.
Returns
-------
range_scalers: DataFrame
J rows, 2 columns: column 1 = the per-tag range (approximately
``q98 - q02`` when ``settings["robust"]`` is True, else raw
``max - min``); column 2 = the per-tag minimum. Both columns are
aggregated across batches with the median when ``robust=True`` and
the mean otherwise, but the per-batch minimum itself is always the
raw ``batch.min(axis=0)``, not a quantile.
"""
# A fit / transform wrapper over this trio now exists as :class:`BatchScaler`, which
# also accepts the melted-DataFrame input these functions reject (#199).
default_settings: dict = {"robust": True, "robust_range": "q98-q02"}
if settings:
default_settings.update(settings)
settings = default_settings
if settings["robust_range"] not in _ROBUST_RANGE_QUANTILES:
raise ValueError(
f"settings['robust_range']={settings['robust_range']!r} is not recognized; expected "
f"one of {sorted(_ROBUST_RANGE_QUANTILES)}."
)
columns_to_align = _resolve_columns_to_align(batches, columns_to_align, "determine_scaling")
collector_rnge = []
collector_mins = []
collapsed: dict[str, int] = {}
for batch in batches.values():
if settings["robust"]:
lower, upper = _ROBUST_RANGE_QUANTILES[settings["robust_range"]]
rnge = batch[columns_to_align].quantile(upper) - batch[columns_to_align].quantile(lower)
else:
rnge = batch[columns_to_align].max() - batch[columns_to_align].min()
rnge = cast("pd.Series", rnge)
is_zero = rnge.to_numpy() == 0
for tag in rnge.index[is_zero]:
collapsed[str(tag)] = collapsed.get(str(tag), 0) + 1
rnge[is_zero] = 1.0
collector_rnge.append(rnge)
# Restricted to `columns_to_align`, as the range above is. Taking the minimum over
# every column raised `TypeError: Cannot convert [...] to numeric` the moment a
# batch carried a non-numeric column, which a batch identifier of string type
# always is (#197). With integer identifiers it did not raise, but it still put
# rows in the result for columns that were never scaled, leaving `Range` NaN
# against a real `Minimum` for each of them.
collector_mins.append(batch[columns_to_align].min(axis=0))
_warn_about_collapsed_ranges(collapsed, len(batches))
if settings["robust"]:
scalings = pd.concat(
[
pd.DataFrame(collector_rnge).median(),
pd.DataFrame(collector_mins).median(),
],
axis=1,
)
else:
scalings = pd.concat(
[pd.DataFrame(collector_rnge).mean(), pd.DataFrame(collector_mins).mean()],
axis=1,
)
scalings.columns = ["Range", "Minimum"]
return scalings
[docs]
def apply_scaling(
batches: dict[str, pd.DataFrame],
scale_df: pd.DataFrame,
columns_to_align: list | pd.Index | None = None,
) -> dict:
"""Scales the batches according to the information in the scaling dataframe.
Parameters
----------
batches : dict[str, pd.DataFrame]
The batches, in standard format.
scale_df : pd.DataFrame
The scaling dataframe, from `determine_scaling`.
columns_to_align : list, pd.Index, or None, optional
Which columns of each batch to scale. Columns outside this list are
dropped from the output. When ``None`` (the default) the columns of
the first batch are used, matching :func:`determine_scaling`.
Returns
-------
dict
The scaled batch data. Each value carries only the ``columns_to_align``
columns, in that order.
"""
columns_to_align = _resolve_columns_to_align(batches, columns_to_align, "apply_scaling")
out = {}
for batch_id, batch in batches.items():
out[batch_id] = batch[columns_to_align].copy()
for tag, column in out[batch_id].items():
tag = cast("str", tag)
minimum = cast("float", scale_df.loc[tag, "Minimum"])
scale_range = cast("float", scale_df.loc[tag, "Range"])
out[batch_id][tag] = (column - minimum) / scale_range
return out
[docs]
def reverse_scaling(
batches: dict[str, pd.DataFrame],
scale_df: pd.DataFrame,
columns_to_align: list | pd.Index | None = None,
) -> dict:
"""Reverse the scaling applied by `apply_scaling`.
Parameters
----------
batches : dict[str, pd.DataFrame]
The scaled batches, in standard format.
scale_df : pd.DataFrame
The scaling dataframe, from :func:`determine_scaling`.
columns_to_align : list, pd.Index, or None, optional
Which columns of each batch to un-scale. Columns outside this list
are dropped from the output. When ``None`` (the default) the columns
of the first batch are used, matching :func:`apply_scaling`.
Returns
-------
dict
The un-scaled batch data.
"""
columns_to_align = _resolve_columns_to_align(batches, columns_to_align, "reverse_scaling")
out = {}
for batch_id, batch in batches.items():
out[batch_id] = batch[columns_to_align].copy()
for tag, column in out[batch_id].items():
tag = cast("str", tag)
minimum = cast("float", scale_df.loc[tag, "Minimum"])
scale_range = cast("float", scale_df.loc[tag, "Range"])
out[batch_id][tag] = column * scale_range + minimum
return out
[docs]
class DTWresult:
"""Result class."""
def __init__( # noqa: PLR0913
self,
synced: np.ndarray | pd.DataFrame,
penalty_matrix: np.ndarray,
md_path: np.ndarray, # multi-dimensional path through distance mesh D = penalty_matrix
warping_path: np.ndarray,
distance: float,
normalized_distance: float,
):
self.synced = synced
self.penalty_matrix = penalty_matrix
self.md_path = md_path
self.warping_path = warping_path
self.distance = distance
self.normalized_distance = normalized_distance
[docs]
class BatchScaler(TransformerMixin, BaseEstimator):
"""
Range-scale batch trajectories, as a fit / transform estimator.
Wraps the three functions :func:`determine_scaling`, :func:`apply_scaling` and
:func:`reverse_scaling` in the estimator shape the rest of this package uses, so
batch preprocessing composes with :class:`~sklearn.pipeline.Pipeline` and with
``clone`` / ``get_params`` / ``set_params`` the way
:class:`~process_improve.multivariate.MCUVScaler` already does (#199). The three
functions remain public and unchanged; this adds a way to carry the fitted scaling
around as one object instead of threading a ``scale_df`` through every call.
Each tag is mapped to roughly ``[0, 1]`` by subtracting a per-tag minimum and
dividing by a per-tag range, both aggregated across the batches seen in
:meth:`fit`. That is a different normalisation from mean-centring to unit variance:
it preserves the shape of a trajectory within its own operating range, which is what
the alignment distance needs.
Parameters
----------
columns_to_align : list or None, optional
The tags to scale. ``None`` (the default) takes the columns of the first batch,
matching :func:`determine_scaling`.
batch_col : str or None, optional
When set, :meth:`fit` and :meth:`transform` also accept a single melted
DataFrame holding every batch, and split it on this column. This is the
DataFrame input case #199 asked for: the functions reject a DataFrame outright
and tell the caller to split it themselves.
robust : bool, optional
Use a robust per-batch range (the default) rather than ``max - min``.
robust_range : str, optional
Which robust range: ``"q98-q02"`` (the default) or ``"iqr"``. See
:func:`determine_scaling` for what separates them.
Attributes
----------
scale_df_ : pd.DataFrame
The fitted scaling, exactly as :func:`determine_scaling` returns it: a ``Range``
and a ``Minimum`` per tag.
columns_to_align_ : list
The tags actually scaled, resolved during :meth:`fit`.
n_features_in_ : int
Number of tags scaled, for the sklearn contract.
Examples
--------
>>> scaler = BatchScaler(columns_to_align=["Temperature"]) # doctest: +SKIP
>>> scaled = scaler.fit_transform(batches) # doctest: +SKIP
>>> original = scaler.inverse_transform(scaled) # doctest: +SKIP
A melted frame works when the batch column is named:
>>> scaler = BatchScaler(batch_col="batch_id") # doctest: +SKIP
>>> scaled = scaler.fit_transform(melted_frame) # doctest: +SKIP
"""
def __init__(
self,
columns_to_align: list | None = None,
batch_col: str | None = None,
robust: bool = True,
robust_range: str = "q98-q02",
) -> None:
self.columns_to_align = columns_to_align
self.batch_col = batch_col
self.robust = robust
self.robust_range = robust_range
def _as_batches(self, X: dict | pd.DataFrame) -> dict:
"""Accept either the dict format or, when ``batch_col`` is set, a melted frame."""
if isinstance(X, pd.DataFrame):
if self.batch_col is None:
raise TypeError(
f"{type(self).__name__} was given a single DataFrame but no `batch_col`, so it "
"cannot tell which column identifies the batch. Pass "
f"{type(self).__name__}(batch_col=...), or split the frame yourself with "
"`dict(tuple(df.groupby(batch_col)))`."
)
return melted_to_dict(X, batch_id_col=self.batch_col)
return X
[docs]
def fit(self, X: dict | pd.DataFrame, y: object = None) -> BatchScaler: # noqa: ARG002
"""
Determine the per-tag range and minimum from these batches.
``y`` is accepted and ignored, per the sklearn transformer contract.
"""
batches = self._as_batches(X)
settings = {"robust": self.robust, "robust_range": self.robust_range}
self.scale_df_ = determine_scaling(batches, columns_to_align=self.columns_to_align, settings=settings)
self.columns_to_align_ = list(
_resolve_columns_to_align(batches, self.columns_to_align, f"{type(self).__name__}.fit")
)
self.n_features_in_ = len(self.columns_to_align_)
return self
[docs]
def align_with_path(md_path: np.ndarray, batch: pd.DataFrame) -> pd.DataFrame:
"""Align a batch to the reference using the DTW path.
Where several samples of ``batch`` map to the same reference index (a
compression in the warping path), the synced value for that index is the
average of those batch samples. The running ``temp`` accumulator is therefore
seeded with the first batch sample for the current index - the same value
assigned to ``synced`` row 0 just below - not with a reference row. A former
``initial_row`` argument seeded it from the reference row (in one caller) or
from an out-of-space batch index (in the other), which mixed an unrelated row
into the row-0 average (#197).
Non-numeric columns are carried through rather than averaged. A batch frame from
:func:`~process_improve.batch.data_input.melted_to_dict` still holds its identifier
column, and averaging a label is meaningless even when it happens to be a number:
with string identifiers it raised ``TypeError: unsupported operand type(s) for /``,
and with integer ones it silently wrote the mean of the identifier into the aligned
frame. Such a column is constant within a batch, so the first value is taken (#197).
"""
numeric = batch.select_dtypes(include="number")
passthrough = batch.columns.difference(numeric.columns, sort=False)
row = 0
nr = md_path[:, 0].max() + 1 # to account for the zero-based indexing
synced = pd.DataFrame(np.zeros((nr, numeric.shape[1])), columns=numeric.columns)
synced.iloc[row, :] = numeric.iloc[md_path[0, 1], :]
temp: pd.Series | np.ndarray = numeric.iloc[md_path[0, 1], :]
for idx in np.arange(1, md_path.shape[0]):
if md_path[idx, 0] != md_path[idx - 1, 0]:
row += 1
synced.iloc[row, :] = temp = numeric.iloc[md_path[idx, 1], :]
else:
# More than one batch sample maps to this reference index (a compression in
# the warping path), so the synced value is the average of those samples.
# Pinned by tests/batch/test_dtw_align_with_path.py.
temp = np.vstack((temp, numeric.iloc[md_path[idx, 1], :]))
synced.iloc[row, :] = np.nanmean(temp, axis=0)
for column in passthrough:
# Constant within a batch, so every row gets the same value and the column keeps
# its own dtype instead of being coerced into the float frame.
synced[column] = batch[column].iloc[0]
return pd.DataFrame(synced[batch.columns])
[docs]
def dtw_core(
test: pd.DataFrame,
ref: pd.DataFrame,
weight_matrix: np.ndarray,
band: object = None,
) -> DTWresult:
"""
Compute DTW alignment of test batch against reference batch.
``band`` is an optional constraint on the warping path, resolved by
:func:`~process_improve.batch.alignment_helpers.resolve_band`. The default
``None`` places no constraint.
"""
nt = test.shape[0] # 'test' data; will be align to the 'reference' data
nr = ref.shape[0]
if test.shape[1] != ref.shape[1]:
raise ValueError(
f"test and ref must have the same number of columns; "
f"got test.shape[1]={test.shape[1]}, ref.shape[1]={ref.shape[1]}."
)
D = distance_matrix(test.values, ref.values, weight_matrix, band=band)
md_path, distance = backtrack_optimal_path(D)
warping_path = np.zeros(nr)
for idx in range(nr):
warping_path[idx] = md_path[np.where(md_path[:, 0] == idx)[0][-1], 1]
# Now align the `test` batch:
synced = align_with_path(md_path=md_path, batch=test)
return DTWresult(
synced,
D,
md_path,
warping_path,
distance,
normalized_distance=distance / (nr + nt),
)
#: Huber's tuning constant, the value giving 95% efficiency at the Gaussian.
_HUBER_CUTOFF = 1.345
#: Smallest batch weight :func:`_batch_weights` will return, so no batch is ever silenced.
_MIN_BATCH_WEIGHT = 1e-3
def _batch_weights(aligned_batches: dict, batch_weighting: str) -> np.ndarray:
"""
Weight each batch by how well it aligned, for the variable-weight update.
Every batch contributes equally to the variable weights under ``"equal"``, so one
badly aligned batch inflates the summed deviation of whichever variables it misfits
and depresses their weights for every other batch. ``"huber"`` downweights such a
batch in proportion to how far its alignment distance sits from the rest.
The distance used is each batch's ``normalized_distance``, which divides by the
summed path length and so compares across batches of unequal duration. It is turned
into a robust z-score against the median and the MAD of the batch set, then passed
through Huber's function: weight 1 inside the cutoff, falling off as ``1 / |z|``
outside it.
Huber rather than a redescending function (Tukey's bisquare, say) because it never
reaches zero. A batch that is downweighted pulls the average trajectory away from
itself, which makes it look worse on the next iteration; a weight that can reach
zero turns that feedback into a one-way door, where a batch excluded once can never
return. The weights are recomputed from scratch every iteration and floored at
``_MIN_BATCH_WEIGHT`` for the same reason.
Parameters
----------
aligned_batches : dict
The :class:`DTWresult` of each batch from the current iteration.
batch_weighting : str
``"equal"`` (every batch weight 1.0) or ``"huber"``.
Returns
-------
np.ndarray
One weight per batch, in the iteration order of ``aligned_batches``, scaled so
they average 1.0. That scaling keeps the accumulated deviations on the magnitude
they had under equal weighting, so the existing relative floor on them, and the
convergence tolerance, keep their meaning.
"""
n_batches = len(aligned_batches)
if batch_weighting == "equal" or n_batches == 0:
return np.ones(n_batches)
distances = np.array([result.normalized_distance for result in aligned_batches.values()], dtype=float)
finite = np.isfinite(distances)
if not finite.any():
return np.ones(n_batches)
median = float(np.median(distances[finite]))
# 1.4826 scales the MAD to estimate the standard deviation of a Gaussian.
mad = 1.4826 * float(np.median(np.abs(distances[finite] - median)))
if mad <= epsqrt:
# The batches are indistinguishable on this measure (or over half of them share
# one distance exactly), so there is nothing to tell apart. MAD can be zero on
# data that does vary, which is why this returns equal weights rather than
# dividing by a floored MAD and manufacturing a spread.
return np.ones(n_batches)
z_scores = np.abs(distances - median) / mad
weights = np.where(z_scores <= _HUBER_CUTOFF, 1.0, _HUBER_CUTOFF / np.maximum(z_scores, epsqrt))
weights = np.where(finite, weights, _MIN_BATCH_WEIGHT)
weights = np.maximum(weights, _MIN_BATCH_WEIGHT)
return weights / float(np.mean(weights))
def _accumulate_deviations(
aligned_batches: dict,
average_batch: pd.DataFrame,
weighting: str,
n_columns: int,
batch_weighting: str = "equal",
) -> np.ndarray:
"""
Sum each variable's deviation from the average trajectory, across all batches.
The reciprocal of this becomes the variable's alignment weight, so a variable that
tracks the average trajectory consistently earns a large weight.
Parameters
----------
aligned_batches : dict
The :class:`DTWresult` of each batch from the current iteration. Every
``synced`` frame sits on the reference grid, so all batches contribute the same
number of rows and no length correction is needed.
average_batch : pd.DataFrame
The average trajectory of the current iteration.
weighting : str
``"quadratic"`` for the sum of squared deviations, which makes the reciprocal an
inverse-variance (precision) weight and matches the Mahalanobis form of the
weighted DTW distance; ``"absolute"`` for the sum of absolute deviations, which
is less sensitive to one badly aligned batch but is not a precision, and changes
both the fixed point and the number of iterations to reach it.
n_columns : int
Number of columns being aligned.
batch_weighting : str, optional
How much each batch contributes, resolved by :func:`_batch_weights`. The default
``"equal"`` gives every batch a weight of exactly 1.0, which is the behaviour
this function had before batch weighting existed.
Returns
-------
np.ndarray
Row vector, one accumulated deviation per column.
"""
accumulated = np.zeros((1, n_columns))
square = weighting == "quadratic"
batch_weights = _batch_weights(aligned_batches, batch_weighting)
for weight, result in zip(batch_weights, aligned_batches.values(), strict=True):
deviation = result.synced - average_batch
term = np.power(deviation, 2) if square else np.abs(deviation)
accumulated = accumulated + weight * np.nansum(term, axis=0)
return accumulated
[docs]
def one_iteration_dtw(
batches_scaled: dict,
refbatch_sc: pd.DataFrame,
weight_matrix: np.ndarray,
settings: dict | None = None,
) -> tuple[dict, pd.DataFrame]:
"""Perform one iteration of the DTW alignment algorithm."""
default_settings: dict = {"show_progress": True, "subsample": 1, "band": None}
if settings:
default_settings.update(settings)
settings = default_settings
aligned_batches = {}
average_batch = refbatch_sc.copy().reset_index(drop=True) * 0.0
successful_alignments = 0
for batch_id, batch in tqdm(batches_scaled.items(), disable=not (settings["show_progress"])):
try:
# see Kassidas, page 180
batch_subset = batch.iloc[:: int(settings["subsample"]), :]
result = dtw_core(batch_subset, refbatch_sc, weight_matrix=weight_matrix, band=settings["band"])
average_batch = average_batch + result.synced
aligned_batches[batch_id] = result
successful_alignments += 1
except ValueError as exc: # noqa: PERF203
# Chain the cause: a band constraint that is too narrow explains itself, and
# that explanation is the actionable part. `from None` discarded it.
raise ValueError(f"Failed on batch {batch_id}: {exc}") from exc
average_batch = average_batch / successful_alignments
return aligned_batches, average_batch
def _validate_time_axis(settings: dict) -> None:
"""
Check the resampled time axis is well formed, before any alignment work is done.
The axis has to hold more than one point for the interpolation to have anything to
interpolate along, and both values have to be positive for ``np.arange`` to produce
an increasing axis at all.
"""
maximum = float(settings["interpolate_time_axis_maximum"])
delta = float(settings["interpolate_time_axis_delta"])
if delta <= 0 or maximum <= 0:
raise ValueError(
f"settings['interpolate_time_axis_maximum'] and ['interpolate_time_axis_delta'] must both "
f"be positive; got maximum={maximum} and delta={delta}."
)
if delta >= maximum:
raise ValueError(
f"settings['interpolate_time_axis_delta']={delta} must be smaller than "
f"['interpolate_time_axis_maximum']={maximum}, so that the resampled axis has more than "
"one point."
)
def _validate_dtw_settings(settings: dict) -> None:
"""
Reject an unusable ``batch_dtw`` setting once, before any alignment work is done.
Each of these would otherwise surface far from its cause: an unrecognized string
silently taking the other branch, or a band constraint failing once per batch.
"""
if settings["weighting"] not in {"quadratic", "absolute"}:
raise ValueError(
f"settings['weighting']={settings['weighting']!r} is not recognized; expected "
"'quadratic' (the default, and the published method) or 'absolute'."
)
if settings["batch_weighting"] not in {"equal", "huber"}:
raise ValueError(
f"settings['batch_weighting']={settings['batch_weighting']!r} is not recognized; "
"expected 'equal' (the default) or 'huber'."
)
_validate_time_axis(settings)
band = settings["band"]
if not (band is None or callable(band) or isinstance(band, np.ndarray)):
raise TypeError(
f"settings['band']={band!r} is not a band constraint; expected None, an "
"(n_test, 2) array of row bounds, or a callable of (n_test, n_ref). See "
"`alignment_helpers.sakoe_chiba` and `alignment_helpers.itakura`."
)
[docs]
def batch_dtw( # noqa: C901, PLR0915
batches: dict[str, pd.DataFrame],
columns_to_align: list,
reference_batch: str,
settings: dict | None = None,
) -> dict:
"""
Synchronize, via iterative DTW, with weighting.
Algorithm: Kassidas et al. (2004): https://doi.org/10.1002/aic.690440412
Parameters
----------
batches : dict[str, pd.DataFrame]
Batch data, in the standard format.
columns_to_align : list
Which columns to use during the alignment process. The others are aligned, but
get no weight, and therefore do not influence the objective function.
reference_batch : str
Which key in the `batches` is the reference batch to use.
settings : dict
Default settings are::
{
"maximum_iterations": 25, # stops here, even if not converged
"tolerance": 0.1, # convergence tolerance
"robust": True, # use robust scaling
"show_progress": True, # show progress
"subsample": 1, # use every sample
"weighting": "quadratic", # "quadratic" or "absolute"; see below
"batch_weighting": "equal", # "equal" or "huber"; see below
"band": None, # warping-path constraint; see below
"interpolate_time_axis_maximum": 100, # resample time axis to this scale
"interpolate_time_axis_delta": 1, # resolution of the resampled axis
"interpolate_method": "cubic", # any scipy.interpolate.interp1d method
}
The default settings resample the time axis to 100 points, starting at 0 and ending
at 99, so each point is one percent of the batch's duration however long the batch
actually ran. Lower the delta for a finer axis (``0.5`` gives 200 points, ``0.25``
gives 400) or change the maximum for a different scale. The delta no longer has to
divide the maximum exactly: values such as ``0.3`` or ``7`` used to fail an
assertion.
``weighting`` selects how a variable's deviation from the average trajectory is
accumulated before the weight is taken as its reciprocal:
``"quadratic"`` (the default)
The sum of squared deviations, as in Kassidas et al. The reciprocal is then
an inverse-variance (precision) weight, which is what the weighted distance
in :func:`~process_improve.batch.alignment_helpers.distance_matrix` expects:
that distance is a Mahalanobis form, quadratic in the deviations.
``"absolute"``
The sum of absolute deviations. Less sensitive to a single badly aligned
batch, but the reciprocal is no longer a precision, so the weighted distance
loses its Mahalanobis reading. It does not simply flatten the weighting:
on the bundled dryer data the ratio of largest to smallest weight rose from
2.8 to 6.0 and the iteration count from 2 to 3, so both the fixed point and
the path to it differ. Offered for comparison; it is not the published
method, and the effect on your own data should be measured rather than
assumed.
``batch_weighting`` decides how much each batch contributes to the variable
weights. Under ``"equal"`` (the default) every batch counts the same, so one badly
aligned batch inflates the summed deviation of whichever variables it misfits and
depresses their weights for every other batch. ``"huber"`` weights each batch by
Huber's function applied to the robust z-score of its ``normalized_distance``,
against the median and MAD of the batch set: weight 1 inside a cutoff of 1.345,
falling off as ``1 / |z|`` beyond it, then rescaled to average 1.0.
Huber rather than a redescending function because it never reaches zero. A
downweighted batch pulls the average trajectory away from itself, so it looks
worse on the next iteration; a weight that could reach zero would make that a
one-way door. The weights are recomputed from scratch each iteration and floored,
so a batch that recovers is counted again.
``band`` constrains the warping path: the reference rows each test sample may
map to. ``None`` (the default) places no constraint. Pass an ``(n_test, 2)``
array of half-open row bounds, or a callable of ``(n_test, n_ref)`` returning
one, since the two lengths differ from batch to batch and are not known until
each pair is aligned::
from process_improve.batch.alignment_helpers import sakoe_chiba, itakura
settings = {"band": sakoe_chiba(window=0.1)} # 10% of the batch duration
settings = {"band": itakura(max_slope=2.0)}
A constraint speeds up the dynamic programme, from ``O(n_ref * n_test)`` to
``O(window * n_test)``, which matters because every batch is re-aligned on
every iteration. It also changes the result: a corridor that excludes the true
warp changes the aligned trajectories, so the iterated average converges to a
different fixed point. Widen it until the alignment stops changing.
Returns
-------
dict
Various outputs relevant to the alignment, keyed by ``scale_df``,
``aligned_batch_objects``, ``aligned_batch_dfdict``, ``last_average_batch``,
``weight_history`` and ``distances``.
``distances`` is a DataFrame indexed by batch identifier, with the ``Distance``
and ``Normalized distance`` of each batch to the reference on the final
iteration. ``Normalized distance`` divides by the summed path length, so it is
comparable across batches of unequal duration. Use it to see which batches
aligned poorly, for instance ``outputs["distances"]["Normalized distance"]
.nlargest(5)``.
Notation
--------
I = number of batches: index = i
i = index for the batches
J = number of tags (columns in each batch)
j = index for the tags
k = index into the rows of each batch, the samples: 0 ... k ... K_i
"""
default_settings: dict = dict(
maximum_iterations=25, # maximum iterations (stops here, even if not converged)
tolerance=0.1, # convergence tolerance
robust=True, # use robust scaling
show_progress=True, # show progress
subsample=1, # use every sample
weighting="quadratic", # how to accumulate deviations: "quadratic" or "absolute"
batch_weighting="equal", # how much each batch contributes: "equal" or "huber"
band=None, # warping-path constraint; None places none. See `sakoe_chiba`, `itakura`.
interpolate_time_axis_maximum=100, # interpolates everything to be on this scale
interpolate_time_axis_delta=1,
interpolate_method="cubic", # any method from scipy.interpolate.interp1d allowed
)
if settings:
default_settings.update(settings)
settings = default_settings
if settings["maximum_iterations"] < 3:
raise ValueError(
f"At least 3 iterations are required; got maximum_iterations={settings['maximum_iterations']}."
)
_validate_dtw_settings(settings)
if reference_batch not in batches:
raise KeyError(f"`reference_batch` was not found in the dict of batches; got {reference_batch!r}.")
if not check_valid_batch_dict(
{k: v[columns_to_align] for k, v in batches.items()},
no_nan=True,
):
raise ValueError("One or more batches in the input dict failed validation.")
scale_df = determine_scaling(batches=batches, columns_to_align=columns_to_align, settings=settings)
batches_scaled = apply_scaling(batches, scale_df, columns_to_align)
refbatch_sc = batches_scaled[reference_batch].iloc[:: int(settings["subsample"]), :]
weight_vector = np.ones(refbatch_sc.shape[1])
weight_matrix = np.diag(weight_vector)
weight_history = np.zeros_like(weight_vector) * np.nan
average_batch = None
delta_weight: np.floating | np.ndarray = np.linalg.norm(weight_vector)
iter_step = 0
while (np.linalg.norm(delta_weight) > settings["tolerance"]) and (iter_step <= settings["maximum_iterations"]):
if settings["show_progress"]:
print(f"Iter = {iter_step} and norm = {np.linalg.norm(delta_weight)}") # noqa: T201
iter_step += 1
logger.debug(
"batch_dtw: iteration %d, weight-delta norm=%g (tolerance=%g)",
iter_step,
float(np.linalg.norm(delta_weight)),
settings["tolerance"],
)
weight_matrix = np.diag(weight_vector)
weight_history = np.vstack((weight_history, weight_vector.copy()))
if iter_step > 3:
refbatch_sc = average_batch
aligned_batches, average_batch = one_iteration_dtw(
batches_scaled=batches_scaled,
refbatch_sc=refbatch_sc,
weight_matrix=weight_matrix,
settings=settings,
)
next_weights = _accumulate_deviations(
aligned_batches,
average_batch,
settings["weighting"],
refbatch_sc.shape[1],
settings["batch_weighting"],
)
# Kassidas: each variable's weight is inversely proportional to its
# summed squared deviation from the average trajectory, so a variable
# that aligns consistently (SSQ ~ 0) must get a LARGE weight. The
# previous guard substituted the magic value 10000 for a near-zero
# SSQ, handing the best-aligned variables a weight of ~1e-4 - the
# exact opposite - and the constant was scale-dependent. Floor the
# SSQ (relative to the largest observed SSQ) instead, so the weight
# stays large but finite. The floor is relative to the largest observed value,
# so it holds for either `weighting` choice despite the SSQ-era name.
ssq_floor = max(epsqrt, 1e-6 * float(np.max(next_weights)))
next_weights = 1.0 / np.maximum(next_weights, ssq_floor)
weight_vector = (next_weights / np.sum(next_weights) * len(columns_to_align)).ravel()
# If change in delta_weight is small, we terminate early; no need to fine-tune excessively.
delta_weight = np.diag(weight_matrix) - weight_vector # old - new
# OK, the weights are found: now use the last iteration's result to get back to original
# scaling for the trajectories
weight_history = weight_history[1:, :]
aligned_df_collection: list[pd.DataFrame] = []
new_time_axis = np.arange(
0,
settings["interpolate_time_axis_maximum"],
settings["interpolate_time_axis_delta"],
)
for batch_id, result in tqdm(
aligned_batches.items(),
desc="Interpolating",
disable=not (settings["show_progress"]),
):
synced = align_with_path(
result.md_path,
batches[batch_id].iloc[:: int(settings["subsample"]), :],
)
# Resample the trajectories of the aligned data now along this sequence, whose
# endpoints are taken from the target axis rather than recomputed. Deriving them
# separately assumed the delta divided the maximum exactly: `np.arange` stops at
# the last multiple below the maximum, while `maximum - delta` does not, so a
# delta of 0.3 or 7 put the two axes' endpoints in different places. Two asserts
# caught that as a bare AssertionError, and `python -O` strips asserts, so under
# optimisation it silently extrapolated instead (#197). Sharing the endpoints
# makes any delta work and leaves nothing to assert.
sequence = np.linspace(new_time_axis[0], new_time_axis[-1], synced.shape[0])
synced_interpolated = pd.DataFrame()
for column in synced:
if column in ["batch_id", "_sequence_"]:
continue
interp_column = sp.interpolate.interp1d(
sequence,
synced[column],
kind=settings["interpolate_method"],
assume_sorted=True,
)
synced_interpolated[column] = interp_column(new_time_axis)
# Pop in an extra column at the start of the df
synced_interpolated.insert(0, "batch_id", batch_id)
synced_interpolated.set_index(new_time_axis)
# Overwrite existing dataframe with this, unscaled, and interpolated dataframe.
result.synced = synced_interpolated
aligned_df_collection.append(synced_interpolated)
# Make the batch_id label consistent
aligned_df = pd.concat(aligned_df_collection)
aligned_df["batch_id"] = aligned_df["batch_id"].astype(type(batch_id))
last_average_batch = reverse_scaling(dict(avg=cast("pd.DataFrame", average_batch)), scale_df)["avg"]
aligned_batch_dfdict = melted_to_dict(aligned_df, batch_id_col="batch_id")
# Each DTWresult already carries its distance to the reference, so the per-batch
# distances need no extra computation: they are the final iteration's values.
distances = pd.DataFrame(
[
{
"batch_id": batch_id,
"Distance": result.distance,
"Normalized distance": result.normalized_distance,
}
for batch_id, result in aligned_batches.items()
]
).set_index("batch_id")
return dict(
scale_df=scale_df,
aligned_batch_objects=aligned_batches,
aligned_batch_dfdict=aligned_batch_dfdict,
last_average_batch=last_average_batch,
weight_history=pd.DataFrame(weight_history, columns=columns_to_align),
distances=distances,
)
[docs]
def resample_to_reference(
batches: dict[str, pd.DataFrame],
columns_to_align: list,
reference_batch: str,
settings: dict | None = None,
) -> dict:
"""Resamples all `batches` (only the `columns_to_align`) to the duration of batch with
identifier `reference`.
Parameters
----------
batches : dict[str, pd.DataFrame]
Batch data, in the standard format.
columns_to_align : list
Which columns to use. Others are ignored.
reference_batch : str
Which key in the `batches` is the reference batch.
settings : dict, optional
[description], by default None
Returns
-------
dict
Batch data, in the standard format.
"""
default_settings = {
"interpolate_kind": "cubic", # must be a valid "scipy.interpolate.interp1d" `kind`
}
if settings:
default_settings.update(settings)
settings = default_settings
out = {}
target_time = np.arange(0, batches[reference_batch].shape[0])
target_time = target_time / target_time[-1]
for batch_id, batch in batches.items():
to_resample = np.arange(0, batch.shape[0])
to_resample = to_resample / to_resample[-1]
out_df = {}
for column, series in batch.items():
if column in columns_to_align:
out_df[column] = sp.interpolate.interp1d(
to_resample,
series.values,
copy=False,
kind=settings["interpolate_kind"],
)(target_time)
out[batch_id] = pd.DataFrame(out_df)
return out
[docs]
def find_average_length(batches: dict[str, pd.DataFrame], settings: dict | None = None) -> str:
"""
Find the batch in `batches` with the average length.
Parameters
----------
batches : dict[str, pd.DataFrame]
Batch data, in the standard format.
settings : dict
Default settings are::
{"robust": True} # use robust (median) average batch length
Returns
-------
One of the dictionary keys from `batches`.
"""
default_settings = {
"robust": True, # use robust metric to calculate average batch length
}
if settings:
default_settings.update(settings)
settings = default_settings
batch_lengths = pd.Series({batch_id: df.shape[0] for batch_id, df in batches.items()})
if settings["robust"]:
# If multiple batches of the median length, return the last one.
try:
median_match = np.where((batch_lengths == batch_lengths.median()).to_numpy())[0][-1]
return cast("str", batch_lengths.index[int(median_match)])
except IndexError:
# Very exceptional: if batch_lengths.median() is computed as the average of 2 numbers
# and therefore the median length batch doesn't actually exist
return find_average_length(batches, settings=dict(robust=False))
else:
closest = int((batch_lengths - batch_lengths.mean()).abs().argmin())
return cast("str", batch_lengths.index[closest])
[docs]
def find_reference_batch(
batches: dict[str, pd.DataFrame],
columns_to_align: list,
settings: dict | None = None,
) -> str | list[str]:
"""
Find a reference batch. Assumes NO missing data.
Starts with the average duration batch; resamples (simple interpolation) of all batches to
that duration. Unfolds that resampled data. Does PCA on the wide, unfolded data. Fits,
by default, 4 components. Excludes all batches with Hotelling's T2 > 90% limit. Refits PCA
with 4 components. Finds the batch which has the multivariate combination of scores which are
the smallest (i.e. closest to the model center) and ensures this batch has SPE < 50% of the
model limit.
Parameters
----------
batches : dict[str, pd.DataFrame]
Batch data, in the standard format.
columns_to_align : list
Which columns to use. Others are ignored.
settings : dict, optional
Default settings are::
{
"robust": True, # use robust scaling
"subsample": 1, # use every sample
"method": "pca_most_average", # most average batch from a crude PCA
"n_components": 4,
"number_of_reference_batches": 1, # only a single batch returned
}
Returns
-------
str or list[str]
When ``settings["number_of_reference_batches"] == 1`` (the default), a
single dictionary key from ``batches`` is returned. When more than one
reference batch is requested, a list of that many keys is returned,
ordered from most to least central in the PCA model.
"""
default_settings: dict[str, int | float | str | bool] = {
"robust": True, # use robust scaling
"subsample": 1, # use every sample
"method": "pca_most_average",
"n_components": 4,
"number_of_reference_batches": 1,
}
if isinstance(settings, dict):
default_settings.update(settings)
settings = default_settings
if not isinstance(columns_to_align, list):
raise TypeError(f"`columns_to_align` must be a list of column names; got {type(columns_to_align).__name__}.")
if not check_valid_batch_dict({k: v[columns_to_align] for k, v in batches.items()}):
raise ValueError("One or more batches in the input dict failed validation.")
# Starts with the average duration batch.
initial_reference_id = find_average_length(batches, settings)
# Resamples (simple interpolation) of all batches to that duration.
resampled = resample_to_reference(
batches,
columns_to_align,
reference_batch=initial_reference_id,
settings=settings,
)
# Unfolds that resampled data.
basewide = dict_to_wide(resampled)
# Does PCA on the wide, unfolded data. A=4
scaler = MCUVScaler().fit(basewide)
mcuv = scaler.fit_transform(basewide)
# You can't fit more components than rows in the matrix.
n_components = min(int(np.floor(settings["n_components"])), basewide.shape[0] - 1)
pca_first = PCA(n_components=n_components).fit(mcuv)
# Excludes all batches with Hotelling's T2 > 90% limit.
hotellings_t2_limit_90 = pca_first.hotellings_t2_limit(0.90)
to_keep = pca_first.hotellings_t2_.iloc[:, -1] < hotellings_t2_limit_90
# Refits PCA with A=4 on a subset of the batches, to avoid biasing the PCA model too much.
basewide = basewide.loc[to_keep, :]
scaler = MCUVScaler().fit(basewide)
mcuv = scaler.fit_transform(basewide)
pca_second = PCA(n_components=n_components).fit(mcuv)
# Finds batch with scores; and ensures this batch has SPE < 50% of the model limit.
metrics = pd.DataFrame(
{
"HT2": pca_second.hotellings_t2_.iloc[:, -1],
"SPE": pca_second.spe_.iloc[:, -1],
}
)
metrics = metrics.sort_values(by=["HT2", "SPE"])
requested = int(settings["number_of_reference_batches"])
if requested < 1:
raise ValueError(f"number_of_reference_batches must be >= 1, got {requested}.")
if requested > len(metrics):
# SEC-13 (#261): without this guard the cutoff-relaxation loop below
# walks past ``conf_level=1.0``, which trips an ``assert`` inside
# ``spe_calculation`` and was -O-strippable. Fail fast with a clear
# message instead.
raise ValueError(
f"number_of_reference_batches={requested} exceeds the number of "
f"candidate batches ({len(metrics)}); cannot select that many."
)
# Boolean-mask indexing instead of a DataFrame.query() expression string
# built by f-string (no expression is assembled or evaluated).
spe_metrics = metrics[metrics["SPE"] < pca_second.spe_limit(conf_level=0.5)]
# SEC-13 (#261): bound the cutoff strictly below 1.0. If the loop runs out
# of headroom before enough batches pass, give up with a clear error
# rather than tripping the inner ``assert conf_level < 1.0`` (which is
# ``python -O``-strippable).
start_cutoff = 0.5
max_cutoff = 0.95
while spe_metrics.shape[0] < requested and start_cutoff <= max_cutoff:
spe_metrics = metrics[metrics["SPE"] < pca_second.spe_limit(conf_level=start_cutoff)]
start_cutoff += 0.05
if spe_metrics.shape[0] < requested:
raise ValueError(
f"Could not find {requested} reference batches even at "
f"conf_level={max_cutoff:.2f}; only {spe_metrics.shape[0]} "
"batches passed the SPE cutoff."
)
if requested == 1:
return spe_metrics.index[0] # returns a single entry from the index
return spe_metrics.index[0:requested].to_list()
[docs]
def unfold_blocks(
blocks: dict[str, dict],
*,
initial_conditions: pd.DataFrame | None = None,
group_by_batch: bool = False,
) -> dict[str, pd.DataFrame]:
"""Unfold several aligned batch blocks batchwise, ready for a multi-block model (#193).
:func:`~process_improve.batch.data_input.dict_to_wide` already unfolds one
block of aligned batches into a one-row-per-batch matrix. What it cannot do
is keep several blocks side by side and *separate*, which is what
:meth:`~process_improve.multivariate.methods.MBPCA.fit` and
:meth:`~process_improve.multivariate.methods.MBPLS.fit` want: they take a
``dict[str, pd.DataFrame]`` and preprocess each block on its own.
``BatchPCA`` and ``BatchPLS`` unfold too, but they concatenate the
initial-conditions block onto the trajectories to make a single wide frame,
because the model underneath them is single-block. Here the blocks stay
apart, so a block's own variance decides its weight in the fit rather than
its column count deciding it by accident.
Parameters
----------
blocks : dict[str, dict]
One entry per block. Keys are block names, carried through to the
result. Values are standard batch-data dictionaries: keys are batch
identifiers, values are per-batch dataframes with identical numeric
columns and the same number of rows within a block. Blocks may have
different numbers of columns and different trajectory lengths from one
another; they only have to describe the same batches.
initial_conditions : pd.DataFrame, optional
One row per batch, indexed by batch identifier: measurements taken
before the batch ran, which have no time axis. Added as its own block
under the name ``"initial_conditions"``, not glued onto a trajectory
block.
group_by_batch : bool, optional
Passed to :func:`dict_to_wide` for every block. ``False`` (default)
orders each block's columns ``(tag, sequence)``; ``True`` swaps them to
``(sequence, tag)``.
Returns
-------
dict[str, pd.DataFrame]
One wide dataframe per block, every one sharing the same row index in
the same order, so row *i* is the same batch in every block.
Raises
------
ValueError
If ``blocks`` is empty, if a block name collides with the
``initial_conditions`` block, or if the blocks do not all cover exactly
the same batches.
Examples
--------
>>> wide = unfold_blocks({"spectra": spectra, "process": process}) # doctest: +SKIP
>>> MBPCA(n_components=2).fit(wide) # doctest: +SKIP
"""
if not blocks:
raise ValueError("At least one block is required.")
ic_name = "initial_conditions"
if initial_conditions is not None and ic_name in blocks:
raise ValueError(
f"Block name {ic_name!r} is reserved for the `initial_conditions` argument; "
"rename that block or pass its data as `initial_conditions`."
)
unfolded: dict[str, pd.DataFrame] = {}
for name, batches in blocks.items():
check_valid_batch_dict(batches, no_nan=True)
unfolded[name] = dict_to_wide(batches, group_by_batch=group_by_batch)
# Every block must describe the same batches. Comparing as sets first gives a
# message naming what is missing where; the reindex below then puts them all
# in one order, since `dict_to_wide` sorts by batch id and a caller's dicts
# need not have been built in the same order.
reference_name, reference = next(iter(unfolded.items()))
expected = set(reference.index)
for name, wide in unfolded.items():
if set(wide.index) != expected:
missing = expected - set(wide.index)
extra = set(wide.index) - expected
raise ValueError(
f"Every block must cover the same batches. Block {name!r} differs from "
f"{reference_name!r}: missing {sorted(missing, key=str)}; "
f"unexpected {sorted(extra, key=str)}."
)
if initial_conditions is not None:
unfolded[ic_name] = _validated_initial_conditions(initial_conditions, reference.index)
# The single place row alignment happens. The trajectory blocks already agree,
# because `dict_to_wide` pivots on batch_id and so comes back sorted; the
# initial-conditions block is the one that arrives in the caller's order and has
# to be moved. Reindexing all of them keeps the guarantee in one line and stops
# it depending on `dict_to_wide` continuing to sort.
order = reference.index
return {name: wide.reindex(order) for name, wide in unfolded.items()}
def _validated_initial_conditions(initial_conditions: pd.DataFrame, batch_ids: pd.Index) -> pd.DataFrame:
"""Check the initial-conditions block covers exactly ``batch_ids``, and is numeric and complete.
Row *order* is not this function's business: ``unfold_blocks`` puts every
block into the reference order on the way out, this one included.
"""
if not isinstance(initial_conditions, pd.DataFrame):
raise TypeError(
"initial_conditions must be a pandas DataFrame indexed by batch identifier; "
f"got {type(initial_conditions).__name__}."
)
if set(initial_conditions.index) != set(batch_ids):
missing = set(batch_ids) - set(initial_conditions.index)
extra = set(initial_conditions.index) - set(batch_ids)
raise ValueError(
"initial_conditions must have exactly one row per batch. "
f"Missing batch ids: {sorted(missing, key=str)}; unmatched extra ids: {sorted(extra, key=str)}."
)
z_wide = initial_conditions
if z_wide.select_dtypes(include="number").shape[1] != z_wide.shape[1]:
raise ValueError("All initial_conditions columns must be numeric.")
if z_wide.isna().to_numpy().sum() > 0:
raise ValueError("No missing values allowed in initial_conditions.")
return z_wide