# (c) Kevin Dunn, 2010-2026. MIT License. Based on own private work over the years.
"""Shared primitives for the multivariate package (ENG-01).
This leaf module holds the small, dependency-free building blocks used across
the PCA / PLS / TPLS / multiblock implementations: the :data:`DataMatrix` type
alias, the :data:`epsqrt` tolerance, the NIPALS denominator floor helper
:func:`_nz`, and :class:`SpecificationWarning`. Nothing here imports any of the
sibling submodules, so it sits at the base of the dependency graph.
"""
from __future__ import annotations
import functools
import inspect
import warnings
from collections.abc import Callable
from typing import Any, Literal, TypeAlias, cast
import numpy as np
import pandas as pd
from scipy import sparse
from sklearn.model_selection import KFold, RepeatedKFold, RepeatedStratifiedKFold, StratifiedKFold
from sklearn.utils import Bunch
from sklearn.utils.validation import check_array
from process_improve._warnings import user_stacklevel
# Names re-exported to the rest of the package. Declared explicitly so CodeQL
# does not flag the ``DataMatrix`` type alias (only ever referenced in lazy
# annotation strings under ``from __future__ import annotations``) as unused.
__all__ = [
"Q2_MIN_INCREMENT",
"DataMatrix",
"DiagnoseResult",
"NotEnoughVarianceError",
"SelectionRule",
"SpecificationWarning",
"UncentredDataWarning",
"check_response_block",
"components_with_spread",
"constant_columns",
"epsqrt",
"informative_widths",
"text_table",
"unobserved_columns",
"user_stacklevel",
"y_is_exhausted",
"y_reach",
]
#: Component-selection rules supported by ``PLS.select_n_components``.
#:
#: ``"1se"`` (the current default) applies the one-standard-error rule of
#: Breiman, Friedman, Olshen & Stone (1984, *CART* sec.3.4.3), endorsed by
#: Hastie, Tibshirani & Friedman (*The Elements of Statistical Learning*,
#: sec.7.10): pick the smallest component count whose mean cross-validated
#: error is within one standard error of the lowest. Needs per-fold errors,
#: so repeated K-fold is recommended.
#:
#: ``"min"`` returns the component count with the lowest cross-validated error.
#: On data sets where the validated error keeps drifting down by noise-level
#: amounts after the systematic components are exhausted, this routinely runs
#: to (or near) the maximum component count.
#:
#: ``"q2_increment"`` keeps adding components while each one raises the
#: cumulative cross-validated :math:`Q^2` by at least ``min_q2_increase``;
#: it is a Wold's-R-style heuristic with an absolute (not relative) threshold.
#:
#: ``"randomization"`` is Van der Voet's (1994, *Chemom. Intell. Lab. Syst.*
#: 25(2):313-323) randomization test: under the null that two PLS models have
#: the same predictive ability, randomly flip the sign of each observation's
#: paired squared-residual difference between a candidate model and the
#: argmin-RMSECV reference. The smallest model whose *p*-value exceeds
#: ``alpha`` (i.e. fails to reject the null) is chosen - statistically
#: indistinguishable from the best, but more parsimonious.
SelectionRule = Literal["1se", "min", "q2_increment", "randomization"]
#: Default minimum increase in the cross-validated :math:`Q^2` for an extra
#: component to be judged worth keeping under the ``"q2_increment"`` selection
#: rule. ``0.01`` means "must add at least one percentage point of
#: cross-validated explained variance".
Q2_MIN_INCREMENT = 0.01
DataMatrix: TypeAlias = np.ndarray | pd.DataFrame
epsqrt = np.sqrt(np.finfo(float).eps)
def components_with_spread(score_sd: np.ndarray) -> np.ndarray:
"""Return which components hold some of the data's spread: the one rule every Hotelling's T² uses (#1026, #1124).
A component whose training score standard deviation is at most ``epsqrt`` times the largest
holds none of it: the last component of a rank-deficient fit has an sd near 1e-16 times the
first, for instance. Its standardised score is noise over noise for a training row and
enormous for a new row off the model plane, so it adds 0 to T². It is also left out of every
matrix that a T² or a trimmed score regression inverts, which it would make singular.
The floor is relative to the largest sd in any units. It used to be ``epsqrt`` itself while
every sd was below 1, which took every component of data in small units for one with no
spread, and gave them a T² of 0 (#1070).
Parameters
----------
score_sd : np.ndarray of shape (n_components,)
The training score standard deviations (``scaling_factor_for_scores_``, or the square
root of ``explained_variance_``).
Returns
-------
np.ndarray of shape (n_components,), dtype bool
True for each component with spread.
"""
score_sd = np.asarray(score_sd, dtype=float)
# fmax skips a NaN sd, as max(1.0, ...) did: that component is left out, the others are not.
return score_sd > epsqrt * float(np.fmax.reduce(score_sd))
def column_sd(values: np.ndarray | pd.DataFrame) -> np.ndarray:
"""Return the sample standard deviation (``ddof=1``) of each column, missing values skipped, in any units.
A square overflows beyond about 1e154 and underflows to 0 below about 1e-162, so a standard
deviation taken in the data's own units was ``inf`` for data near 1e160 and 0 for data near
1e-160 (#1070). Each column is divided first by the power of two just above its largest
absolute value, which brings it within (-1, 1), and the result multiplied back. Dividing by a
power of two is exact, so a column whose standard deviation did not overflow or underflow
gives the same one, bit for bit.
Parameters
----------
values : np.ndarray or pd.DataFrame of shape (n_samples, n_features)
The data, with NaN where a value is missing.
Returns
-------
np.ndarray of shape (n_features,)
The standard deviations, NaN for a column with fewer than two observed values. One beyond
the largest float is inf, with numpy's overflow warning, as ``numpy.nanstd`` gives it.
"""
values = np.asarray(values, dtype=float)
largest = np.fmax.reduce(np.abs(values), axis=0) # fmax skips NaN, and an all-NaN column, without a warning
_, exponent = np.frexp(np.where(np.isfinite(largest) & (largest > 0), largest, 1.0))
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning) # fewer than two observed values: NaN, which is the answer
spread = np.nanstd(np.ldexp(values, -exponent), axis=0, ddof=1)
# ldexp scales by the power of two without forming it: at 2**1023 and above that power is 2**1024, which
# overflowed to inf and gave the column a standard deviation of inf * 0, NaN (#1149).
return np.ldexp(spread, exponent)
#: Smallest positive float; used to floor NIPALS denominators away from zero.
_DENOM_FLOOR = float(np.finfo(float).tiny)
def _nz(denominator: float) -> float:
"""Floor a non-negative NIPALS denominator away from zero.
Sum-of-squares and vector-norm denominators (``v @ v``, ``norm(v)``) are
non-negative. When a score or loading vector collapses to (near) zero during
NIPALS - a fully-deflated component, or a degenerate / perfectly collinear
block - the denominator is ~0 and the division would yield ``inf``/``nan``
that silently poisons the fitted model. Flooring to the smallest positive
float leaves every well-conditioned value untouched (real denominators are
far larger) while turning the degenerate ``0/0`` into a finite (~0)
projection, since the numerator collapses with the same vector.
"""
return max(_DENOM_FLOOR, denominator)
def _own_axes(frame: pd.DataFrame) -> pd.DataFrame:
"""Return ``frame`` with copies of its row and column labels, and its values unchanged and not copied.
pandas passes an index on by reference: a frame built on the caller's ``X.index`` holds
that very object, so a name given to the caller's axis after the fit renamed the
model's frames too. A model takes its labels from here, or from ``Index.copy``, instead.
"""
owned = frame.copy(deep=False)
owned.index = frame.index.copy()
owned.columns = frame.columns.copy()
return owned
#: The ``missing_data_settings`` keys removed in 2.0, and the constructor parameter that
#: replaced each one. ``md_method`` is not here: it names the imputation algorithm, which
#: is the one axis ``missing_data_settings`` still owns.
_REMOVED_MD_KEYS: dict[str, str] = {"md_tol": "tol", "md_max_iter": "max_iter"}
#: The historical bounds on the NIPALS tolerance, kept verbatim from the per-estimator
#: checks this helper replaced. The ceiling is arbitrary but catches a caller who passed a
#: percentage; the floor is where a relative comparison stops being meaningful.
_TOL_CEILING = 10.0
def check_at_least_one(name: str, value: float | None) -> None:
"""Refuse a count that is not a whole number of at least 1, from ``fit`` (#1071).
``__init__`` stores its arguments unchanged, as sklearn requires, and ``set_params`` and
``clone`` do not pass through any check there, so ``fit`` is the one place every route to
a setting reaches. ``None`` passes: an estimator that accepts it reads it as "as many as
the data allow".
Parameters
----------
name : str
The parameter's name, used in the error message.
value : int or None
The value to check.
Raises
------
ValueError
If ``value`` is not a whole number (NaN included), or is below 1.
"""
check_whole(name, value)
if value is not None and value < 1:
raise ValueError(f"{name} must be >= 1; got {value}.")
def check_whole(name: str, value: float | None) -> None:
"""Refuse a count that is not a whole number; ``None`` passes (#1114).
A count of 2.5 was read as 2, without a word, or failed in a slice with a ``TypeError``.
A float with no fractional part is refused too, as scikit-learn's ``Integral`` constraint
refuses it in an estimator's settings.
Parameters
----------
name : str
The parameter's name, used in the error message.
value : int or None
The value to check.
Raises
------
ValueError
If ``value`` is not an integer, NaN included.
"""
if value is not None and not isinstance(value, (int, np.integer)):
raise ValueError(f"{name} must be a whole number; got {value}.")
class FoldsByClass(int):
"""A number of cross-validation folds, each class shuffled into them on its own: PLS-DA's integer ``cv`` (#1149).
PLS-DA runs PLS's cross-validation on its class-indicator matrix. A plain integer shuffles
the rows into folds whatever their class, so a training fold could hold no row of a class:
every column of its indicator matrix was then constant, and the fold's fit refused it, as
a ``nested_cv`` on twelve rows of two classes did for 7 seeds in 40. This integer is the
number of folds everywhere PLS reads one; where PLS shuffles the rows into folds,
:func:`_shuffled_folds` shuffles each class into them on its own, as
:class:`~sklearn.model_selection.StratifiedKFold` does, with the same seed.
"""
__slots__ = ()
def _shuffled_folds(
n_splits: int, X: pd.DataFrame, Y: pd.DataFrame, *, n_repeats: int = 1, random_state: int | None = None
) -> list[tuple[np.ndarray, np.ndarray]]:
"""Return the folds an integer ``cv`` stands for: the rows shuffled into ``n_splits`` folds, ``n_repeats`` times.
A :class:`FoldsByClass` shuffles each class into the folds on its own, the class of a row
being the column of the indicator matrix ``Y`` that marks it; any other count shuffles the
rows whatever their response.
Parameters
----------
n_splits : int
The number of folds.
X, Y : pd.DataFrame
The rows to split, and their response.
n_repeats : int, default=1
How many times the rows are shuffled into folds afresh.
random_state : int or None
The seed of the shuffle, as scikit-learn's splitters take it.
Returns
-------
list of (np.ndarray, np.ndarray)
The training and held-out positions of each fold, repeat by repeat.
"""
if isinstance(n_splits, FoldsByClass):
classes = np.argmax(np.asarray(Y), axis=1)
if n_repeats == 1:
return list(StratifiedKFold(int(n_splits), shuffle=True, random_state=random_state).split(X, classes))
repeated = RepeatedStratifiedKFold(n_splits=int(n_splits), n_repeats=n_repeats, random_state=random_state)
return list(repeated.split(X, classes))
if n_repeats == 1:
return list(KFold(int(n_splits), shuffle=True, random_state=random_state).split(X))
return list(RepeatedKFold(n_splits=int(n_splits), n_repeats=n_repeats, random_state=random_state).split(X))
def components_to_fit(
n_components: int | None, n_rows: int, n_columns: int, *, columns: str = "the number of columns"
) -> int:
"""Return how many components a fit computes: ``n_components``, cut to what the data can give (#505).
No more components than the smaller of ``n_rows`` and ``n_columns`` can be computed. PCA, PLS
and OPLS cut a larger count to it with this warning; the other members of the family went on
into a raw numpy error, a ``RuntimeWarning``, or components fitted to rounding noise (#1149).
``None`` asks for that many. The constructor parameter stays as the caller set it, as sklearn's
``clone`` needs: the count is the fitted ``n_components_``.
Parameters
----------
n_components : int or None
The count asked for.
n_rows, n_columns : int
The two sizes the count cannot exceed.
columns : str
What ``n_columns`` is, as the warning names it.
Returns
-------
int
The count to fit.
Raises
------
ValueError
If ``n_components`` is not a whole number of at least 1.
Warns
-----
SpecificationWarning
If ``n_components`` is more than the data can give.
"""
check_at_least_one("n_components", n_components)
most = int(min(n_rows, n_columns))
if n_components is None:
return most
if most < n_components:
warnings.warn(
"The requested number of components is more than can be computed from data. The maximum number of "
f"components is the minimum of either the number of rows ({n_rows}) or {columns} ({n_columns}).",
SpecificationWarning,
stacklevel=user_stacklevel(),
)
return most
return int(n_components)
def refuse_empty_block(block: DataMatrix, label: str) -> None:
"""Refuse a block with no columns, naming it (#1077).
scikit-learn's check of a DataFrame with no columns fails inside numpy, with "at least one
array or dtype is required", which names neither the block nor the cause. An array is read
the same way, so a single-block estimator can check its input before scikit-learn does (#1114).
Parameters
----------
block : pd.DataFrame or np.ndarray
The block to check.
label : str
Where the block sits, as ``'y'``, ``'X["name"]'`` or ``'Z["group"]'``.
Raises
------
ValueError
If ``block`` has no columns.
"""
if np.shape(block)[1:] == (0,): # 2-D with no columns; anything else is left to the checks that follow
raise ValueError(f"Block {label} has no columns; a block needs at least one.")
def resolve_loop_settings(
*,
tol: float,
max_iter: int,
missing_data_settings: dict | None,
estimator_name: str = "",
) -> dict:
"""Resolve one estimator's ``tol`` and ``max_iter`` into the mapping its NIPALS loop reads.
Every latent-variable estimator that iterates resolves its loop settings here, so the
rule is defined once instead of five times (#588). The constructor's ``tol`` and
``max_iter`` supply the values, and any other keys of a caller's
``missing_data_settings`` are added alongside them. Both are checked here, at fit
time, so a value set through ``set_params`` or ``clone`` is refused too (#1071).
The keys ``md_tol`` and ``md_max_iter``, which used to override ``tol`` and
``max_iter``, were removed in 2.0 and are refused. They described a convergence
setting as a missing-data setting, which it is not: the loop runs, and needs a
tolerance, whether or not a cell is missing. Other unknown keys pass through, but
these two raise, because ignoring a convergence setting a caller asked for would
change the fit without saying so.
Parameters
----------
tol : float
The estimator's convergence tolerance, relative to the norm of the score vector.
max_iter : int
The estimator's per-component iteration cap.
missing_data_settings : dict or None
The caller's dict, or None. Keys other than the removed two are passed through
untouched, so an estimator with its own keys (``md_method`` on ``PLS``) still sees
them.
estimator_name : str
Used in the error message, so the reader knows which class to change.
Returns
-------
dict
The loop settings under the internal keys ``md_tol`` and ``md_max_iter``, plus any
other keys the caller supplied.
Raises
------
ValueError
If ``missing_data_settings`` holds ``md_tol`` or ``md_max_iter``, if ``max_iter``
is below 1, or if the tolerance is out of bounds.
"""
settings: dict = {"md_tol": tol, "md_max_iter": max_iter}
if isinstance(missing_data_settings, dict):
removed = [key for key in _REMOVED_MD_KEYS if key in missing_data_settings]
if removed:
where = f"{estimator_name}(missing_data_settings=...)" if estimator_name else "missing_data_settings"
keys = " and ".join(repr(key) for key in removed)
replacements = ", ".join(f"{_REMOVED_MD_KEYS[key]}={missing_data_settings[key]!r}" for key in removed)
raise ValueError(
f"{where} no longer takes {keys}, removed in 2.0. Pass {replacements} to the constructor instead."
)
settings.update(missing_data_settings)
check_at_least_one("max_iter", max_iter)
settings["md_max_iter"] = int(settings["md_max_iter"])
if not settings["md_tol"] < _TOL_CEILING:
raise ValueError(f"Tolerance should not be too large; got {settings['md_tol']}.")
if not settings["md_tol"] > epsqrt**1.95:
raise ValueError(f"Tolerance must exceed machine precision; got {settings['md_tol']}.")
return settings
def _reject_sparse(X: object, estimator_name: str) -> None:
"""Raise if ``X`` is a SciPy sparse matrix, naming the remedy that actually helps.
NIPALS centres and scales every column, which destroys sparsity, so there is nothing
to gain from a sparse representation here and every estimator in this package sets
``accept_sparse=False``. sklearn's own message for that says to call ``.toarray()``,
which is the expensive way round when the sparsity came from a one-hot block: the
dense array is materialised anyway, but only after the sparse one has been built.
The usual source is a :class:`~sklearn.compose.ColumnTransformer` whose one-hot block
pushed the result over ``sparse_threshold`` (default 0.3), which flips the *whole*
concatenated output to sparse. Two knobs at that level avoid the round trip
altogether, and the message names both (#399).
Parameters
----------
X : object
The candidate input. Anything SciPy does not consider sparse passes through.
estimator_name : str
Named in the message, so a Pipeline failure says which step rejected the input.
Raises
------
TypeError
If ``X`` is sparse.
"""
if not sparse.issparse(X):
return
msg = (
f"{estimator_name} does not accept sparse input: NIPALS centres and scales every "
"column, which destroys sparsity, so there is nothing to gain from it. This usually "
"arrives from a ColumnTransformer whose one-hot block pushed the concatenated output "
"over `sparse_threshold` (default 0.3). Pass `sparse_threshold=0` to the "
"ColumnTransformer, or `OneHotEncoder(sparse_output=False)`, so it hands over a dense "
"array directly. `X.toarray()` also works, but on a wide one-hot block it builds the "
"sparse matrix first and then the dense one anyway."
)
raise TypeError(msg)
[docs]
class SpecificationWarning(UserWarning):
"""Parent warning class."""
class DiagnoseResult(Bunch):
"""The :class:`~sklearn.utils.Bunch` a ``diagnose`` returns; reading its ``hotellings_t2`` warns (#983).
``cumulative_hotellings_t2`` holds Hotelling's T² in one shape from every estimator: a
DataFrame with a row per diagnosed row and a column per component, column ``a`` holding the
T² over the first ``a`` components, as the fitted T² attribute does. ``hotellings_t2`` was
that DataFrame from PCA and TPLS, and its last column, a Series, from PLS, MBPLS and MBPCA.
It still holds what it did, until 3.0. Reading it by name, as an attribute, with ``[]`` or
with ``get``, gives a :class:`DeprecationWarning`; reading the whole result, as ``dict(result)``
or ``result.items()`` do, does not.
Parameters
----------
estimator : str
The estimator's name under ``process_improve``, such as ``"multivariate.PLS"``, for the
warning.
**fields
The fields of the result.
"""
def __init__(self, estimator: str, **fields: Any) -> None: # noqa: ANN401 - a field of any type
super().__init__(**fields)
# Not a field: Bunch would store an attribute as a key.
self.__dict__["_estimator"] = estimator
def _warn_if_deprecated(self, key: object) -> None:
if key == "hotellings_t2":
# An unpickled result has no _estimator: Bunch.__setstate__ restores no attribute.
estimator = self.__dict__.get("_estimator", "multivariate model")
warnings.warn(
f"process_improve.{estimator}.diagnose(...).hotellings_t2 is deprecated since 2.4.0 and will be "
"removed in 3.0; use .cumulative_hotellings_t2 instead, which holds Hotelling's T² in the same "
"shape from every estimator: a column per component, cumulative, as the fitted T² attribute. "
"Its last column, .cumulative_hotellings_t2.iloc[:, -1], is the T² of the full model.",
category=DeprecationWarning,
stacklevel=user_stacklevel(),
)
def __getitem__(self, key: object) -> Any: # noqa: ANN401 - a field of any type
self._warn_if_deprecated(key)
return super().__getitem__(key)
def get(self, key: object, default: Any = None) -> Any: # noqa: ANN401 - as dict.get
self._warn_if_deprecated(key)
return super().get(key, default)
[docs]
class UncentredDataWarning(SpecificationWarning):
"""Emitted when a model that fits no intercept is handed an un-centred block.
Raised by :meth:`PLS.fit <process_improve.multivariate.PLS.fit>` under
``scale=False``, which centres nothing and fits no intercept, so a block
carrying a non-zero mean displaces every prediction.
It exists as its own class so a caller who fits un-centred data on purpose
can permit *this* diagnostic without going blind to the rest. Under a
``filterwarnings = error`` policy::
# pytest.ini / pyproject.toml, or @pytest.mark.filterwarnings on one test
filterwarnings =
error
ignore::process_improve.multivariate.UncentredDataWarning
Narrower still, and with no global filter at all, is the estimator flag:
``PLS(..., scale=False, warn_on_uncentred=False)`` silences the check for
that one model and leaves every other ``SpecificationWarning`` in force.
Subclasses :class:`SpecificationWarning`, so filters and ``pytest.warns``
assertions written against the parent keep matching it.
"""
class NotEnoughVarianceError(RuntimeError, ValueError):
"""Raised when more components are requested than the data's variance supports.
During NIPALS extraction (PCA / PLS) the deflated data array can run out of
variance before the requested number of components is reached, even when that
count is within ``min(N, K)`` - for example with rank-deficient or perfectly
collinear data. The model cannot compute any further components in that case.
Subclasses :class:`RuntimeError` so existing ``except RuntimeError`` handlers
keep working, while callers that want to react specifically (e.g. trimming the
component count during cross-validation) can catch this narrower type. It is a
:class:`ValueError` as well, because it replaced the ``ValueError`` TPLS raised for
data with no variance and the ``LinAlgError`` (a ``ValueError``) PLS and OPLS
raised for a Y with no covariance with X: an ``except ValueError`` written for
those still catches it.
"""
def check_variance_left(block: str, values: np.ndarray, extracted: int, size: float) -> None:
"""Refuse component ``extracted + 1`` when ``block`` has no variance left, relative to its size.
``size`` is the block's norm before the first component. The floor is relative to it, as OPLS
judges it, so the same data in smaller units are not taken for none (#1070). PLS judges X and Y
this way before every component, and MBPLS judges Y (#1128), so a Y that is fitted exactly is
refused by both, in the same words. Missing cells are left out of the sum.
Raises
------
NotEnoughVarianceError
If ``block`` has no variance left.
"""
if np.sqrt(np.nansum(values**2)) <= epsqrt * size:
raise NotEnoughVarianceError(
f"There is no variance left in the data array for {block}: cannot compute any more components "
f"beyond component {extracted}."
)
def _align_to_fit_features(X: pd.DataFrame, fit_feature_names: pd.Index) -> pd.DataFrame:
"""Validate and align new-data columns against the features seen during ``fit``.
PCA / PLS only checked that data passed to ``transform`` / ``predict`` had the
right *number* of columns. A correctly-shaped frame whose columns are renamed
or reordered would otherwise be projected positionally (PCA, ``X.values @ P``)
or silently label-aligned to all-``NaN`` (PLS, ``X @ direct_weights_``),
producing wrong scores with no error raised (issue #195). This helper makes
that consistency explicit, mirroring scikit-learn's ``feature_names_in_``
handling:
* If both the training data and ``X`` carry string feature names, the *set*
of names must match (otherwise :class:`ValueError`); columns supplied in a
different order are reordered to the training order.
* If the training data had names but ``X`` does not (e.g. a bare ndarray was
passed), the columns are taken to correspond positionally and are labelled
with the training names, so downstream label-aligned arithmetic stays
correct rather than collapsing to ``NaN``.
* If the training data itself had no (string) feature names, there is nothing
to validate and ``X`` is returned unchanged.
The caller is expected to have already validated the column *count*.
Parameters
----------
X : pd.DataFrame
New data passed to ``transform`` / ``predict`` (already coerced to a
DataFrame and count-checked by the caller).
fit_feature_names : pd.Index
The ``X.columns`` captured during ``fit`` (``self._feature_names``).
Returns
-------
pd.DataFrame
``X`` with columns aligned to the training feature order.
Raises
------
ValueError
If both sides carry string names but the sets of names differ.
"""
fit_names = list(fit_feature_names)
if not (fit_names and all(isinstance(name, str) for name in fit_names)):
# Training data had no string feature names (e.g. fitted from an ndarray);
# only the column count is meaningful, which the caller already checked.
return X
new_names = list(X.columns)
if not all(isinstance(name, str) for name in new_names):
# New data carries default positional columns (e.g. came in as an
# ndarray). Assume positional correspondence and label it with the
# training names so label-aligned operations behave correctly.
X = X.copy()
X.columns = pd.Index(fit_names)
return X
if set(new_names) != set(fit_names):
missing = [name for name in fit_names if name not in set(new_names)]
unexpected = [name for name in new_names if name not in set(fit_names)]
raise ValueError(
"Feature names of the data passed to predict/transform do not match "
"those seen during fit. "
f"Missing columns: {missing}; unexpected columns: {unexpected}."
)
if new_names != fit_names:
# Same names, different order: reorder to the training order so the
# positional projection lines up.
X = X[fit_names]
return X
def _float_frame(data: DataMatrix) -> pd.DataFrame:
"""Return ``data`` as a DataFrame of floats, a pandas missing value (pd.NA) as NaN, as ``fit`` reads it (#1144).
An array's rows and columns are numbered from 0, and a Series keeps its name, as
``pd.DataFrame`` gives them. A frame of floats is returned as it is: a copy can change
its memory layout, and with it the order sums run in.
"""
frame = pd.DataFrame(data)
return frame if (frame.dtypes == np.float64).all() else frame.astype(float)
def _component_curve(values: np.ndarray | pd.Series | list[float], name: str) -> np.ndarray:
"""Return a per-component curve as a float array whose entry ``i`` belongs to ``i + 1`` components.
Every selection rule below returns a position plus one, which is a component count
only when the curve starts at one component and skips none. An array or list cannot
say otherwise, so it is taken as that; a count that was not evaluated keeps its place
as NaN. A Series can say otherwise, through its index, and is refused unless the
index is ``1..A``: dropping it would turn a curve indexed from 2 into an answer one
count too small, with no sign of it (#513).
Raises
------
ValueError
If ``values`` is a Series whose index is not ``1..A``.
"""
if isinstance(values, pd.Series) and not values.index.equals(pd.RangeIndex(1, len(values) + 1)):
msg = f"{name} must be indexed by component count 1..{len(values)}; got {values.index.tolist()}."
raise ValueError(msg)
return np.asarray(values, dtype=float)
def _recommend_n_components_q2(
q2_cumulative: np.ndarray | pd.Series | list[float],
*,
min_increment: float = Q2_MIN_INCREMENT,
) -> int:
r"""Recommend a component count from a cumulative cross-validated :math:`Q^2` curve.
Walks outward from one component, keeping component ``a`` only if it lifts
the cumulative :math:`Q^2` by at least ``min_increment`` (with an implied
:math:`Q^2` of ``0`` before any component). The first component that fails
that test - a plateau, a drop, or a non-finite entry - stops the search,
and the recommendation is the last component that passed (never fewer than
one, since a model needs at least one component).
This is a Wold's-R-style heuristic: parsimonious and cheap, but the
threshold is absolute and hand-tuned. Prefer the 1-SE rule
(:func:`_recommend_n_components_one_se`) when per-fold errors are available.
Entry ``i`` of the curve must belong to ``i + 1`` components; see
:func:`_component_curve`.
"""
q2 = _component_curve(q2_cumulative, "q2_cumulative")
recommended = 1
previous = 0.0
for a, value in enumerate(q2, start=1):
if not np.isfinite(value) or (value - previous) < min_increment:
break
recommended = a
previous = float(value)
return recommended
def _recommend_n_components_one_se(
mean_error: np.ndarray | pd.Series | list[float],
se_error: np.ndarray | pd.Series | list[float],
) -> int:
r"""Recommend a component count by the one-standard-error rule.
Given the mean cross-validated error per component count (``mean_error``,
e.g. RMSECV across folds and repeats) and its standard error
(``se_error``), find ``a* = nanargmin(mean_error)`` and return the
*smallest* component count ``a`` whose ``mean_error[a]`` is no larger than
``mean_error[a*] + se_error[a*]``.
The 1-SE rule (Breiman, Friedman, Olshen & Stone 1984, *CART* sec.3.4.3;
Hastie, Tibshirani & Friedman, *ESL* sec.7.10) trades a tiny amount of
fit for a more parsimonious model whose error is statistically
indistinguishable from the best. It is the cheapest robustness upgrade for
integer hyperparameters like the number of latent variables.
Non-finite entries in ``mean_error`` are skipped. The selection always
returns at least ``1``. If ``se_error`` at the minimum is non-finite, zero
or negative, the rule degenerates to the argmin. Entry ``i`` of each curve
must belong to ``i + 1`` components; see :func:`_component_curve`.
"""
mean = _component_curve(mean_error, "mean_error")
se = _component_curve(se_error, "se_error")
if mean.shape != se.shape:
raise ValueError(f"mean_error and se_error must have the same shape; got {mean.shape} and {se.shape}.")
if mean.ndim != 1:
raise ValueError("mean_error must be 1-D.")
finite = np.isfinite(mean)
if not finite.any():
# Caller should already have raised before now; fall through to 1 so
# this helper never returns a sentinel.
return 1
masked = np.where(finite, mean, np.inf)
best = int(np.argmin(masked))
# Floored at zero, so the minimum always meets its own threshold and the search
# below always finds a count.
se_best = float(se[best]) if np.isfinite(se[best]) and se[best] > 0 else 0.0
return int(np.flatnonzero(masked <= masked[best] + se_best)[0]) + 1
def _select_n_components(
rule: SelectionRule,
*,
mean_error: np.ndarray | pd.Series | list[float],
se_error: np.ndarray | pd.Series | list[float] | None = None,
q2_cumulative: np.ndarray | pd.Series | list[float] | None = None,
min_q2_increase: float = Q2_MIN_INCREMENT,
) -> int:
"""Dispatch component selection to one of the supported rules.
See :data:`SelectionRule` for the rule semantics. Raises ``ValueError`` if
a rule is requested without the data it needs (``"1se"`` needs
``se_error``; ``"q2_increment"`` needs ``q2_cumulative``), or if a curve it
reads is a Series not indexed ``1..A`` (see :func:`_component_curve`).
"""
if rule == "min":
mean = _component_curve(mean_error, "mean_error")
if not np.isfinite(mean).any():
return 1
return int(np.nanargmin(mean)) + 1
if rule == "1se":
if se_error is None:
raise ValueError("selection_rule='1se' requires per-component standard errors.")
return _recommend_n_components_one_se(mean_error, se_error)
if rule == "q2_increment":
if q2_cumulative is None:
raise ValueError("selection_rule='q2_increment' requires a cumulative Q^2 curve.")
return _recommend_n_components_q2(q2_cumulative, min_increment=min_q2_increase)
raise ValueError(
f"Unknown selection_rule {rule!r}; expected one of "
f"'1se', 'min', 'q2_increment', 'randomization' "
f"('randomization' is handled by PLS.select_n_components, not this dispatcher)."
)
def _vandervoet_randomization(
per_obs_sse: np.ndarray,
*,
total_rmsecv: np.ndarray,
n_permutations: int = 999,
alpha: float = 0.01,
random_state: int | np.random.Generator | None = None,
) -> tuple[int, np.ndarray]:
"""Van der Voet (1994) randomization test for PLS component selection.
Compares every candidate model against the reference (argmin-RMSECV)
model under the null that the two have the same predictive ability.
For each observation the paired difference of squared residuals
``D_i = sse[a, i] - sse[a*, i]`` is computed; under the null its sign
is random, so the permutation distribution of ``T = sum_i D_i`` is
obtained by flipping each ``D_i``'s sign with probability 1/2 over
``n_permutations`` draws. The *p*-value is the right-tail probability
of seeing a sum as large as the observed one (``T_obs >= T_perm``);
the recommendation is the smallest ``a`` whose ``p > alpha`` -
statistically indistinguishable from the reference, but more
parsimonious.
Parameters
----------
per_obs_sse : np.ndarray of shape (n_components, n_samples)
Out-of-fold per-observation squared total residual at every
component count, summed across Y columns. Rows that are NaN
(observation never held out) are dropped.
total_rmsecv : np.ndarray of shape (n_components,)
Pooled total RMSECV per component count; used to pick the
reference model ``a*`` = ``nanargmin(total_rmsecv) + 1``.
n_permutations : int, default 999
Number of sign-flip permutations.
alpha : float, default 0.01
Significance level. Smaller values are more parsimonious.
random_state : int, numpy.random.Generator or None, optional
Seed or generator for reproducible permutations.
Returns
-------
recommended : int
Smallest 1-based component count with ``p > alpha``.
p_values : np.ndarray of shape (n_components,)
Right-tail *p*-value per candidate; the reference model gets
``1.0`` by construction (paired differences are all zero).
References
----------
Van der Voet, H. (1994). Comparing the predictive accuracy of
models using a simple randomization test. *Chemom. Intell. Lab.
Syst.*, 25(2), 313-323.
"""
a_count = per_obs_sse.shape[0]
a_ref = int(np.nanargmin(total_rmsecv))
rng = np.random.default_rng(random_state)
p_values = np.zeros(a_count)
p_values[a_ref] = 1.0
sse_ref = per_obs_sse[a_ref]
for a in range(a_count):
if a == a_ref:
continue
d = per_obs_sse[a] - sse_ref
# Drop observations with NaN (a custom splitter may have left some
# rows unheld), since the paired difference is undefined there.
d = d[np.isfinite(d)]
if d.size == 0:
p_values[a] = 1.0
continue
t_obs = float(d.sum())
signs = rng.choice([-1.0, 1.0], size=(n_permutations, d.size))
t_perm = (signs * d).sum(axis=1)
# Right-tail probability under the null. Add 1 to both numerator
# and denominator (the "permutation test +1" correction) so the
# p-value is strictly positive even at the extreme.
p_values[a] = float((np.sum(t_perm >= t_obs) + 1) / (n_permutations + 1))
recommended = a_ref + 1 # fall back to the reference if nothing qualifies
for a in range(a_count):
if p_values[a] > alpha:
recommended = a + 1
break
return recommended, p_values
def _equal_weight_r2_total(per_target: np.ndarray) -> np.ndarray:
r"""Pool per-target validated :math:`R^2_Y` with every target weighted equally.
The ``"total"`` column of a validated-:math:`R^2_Y` table is
``1 - sum_m PRESS_m / sum_m TSS_m`` on the *original* Y scale, so a target
whose spread is two orders of magnitude larger than its neighbours' decides
the number almost on its own. Autoscaling Y first fixes that, and costs
nothing to compute here: after mean-centring and unit-variance scaling every
target's TSS is the same, so the pooled ratio collapses to the arithmetic
mean of the per-target values,
.. math::
1 - \frac{\sum_m \mathrm{PRESS}_m / s_m^2}{\sum_m \mathrm{TSS}_m / s_m^2}
= \frac{1}{M} \sum_m \left(1 - \frac{\mathrm{PRESS}_m}{\mathrm{TSS}_m}\right).
Targets with no spread contribute ``NaN`` per-target values and are left out
of the mean rather than dragging it to ``NaN``.
Parameters
----------
per_target : np.ndarray
Validated :math:`R^2` per target, shape ``(A, M)``.
Returns
-------
np.ndarray
The equal-weight pooled value per component count, shape ``(A,)``.
"""
values = np.asarray(per_target, dtype=float)
if values.ndim != 2:
raise ValueError(f"per_target must be 2-D (A x M); got shape {values.shape}.")
with warnings.catch_warnings():
# An all-NaN row (every target constant) means "nothing to pool"; NaN
# is the right answer and the RuntimeWarning that says so is noise.
warnings.simplefilter("ignore", RuntimeWarning)
return np.nanmean(values, axis=1) if values.shape[1] else np.full(values.shape[0], np.nan)
def _model_method(fn: Callable[..., Any]) -> Callable[..., Any]:
"""Wrap a module-level ``fn(model, ...)`` as an introspectable instance method.
ENG-05: estimators (PCA / PLS / ...) expose convenience methods such as
``score_plot``, ``spe_limit`` and ``vip`` that forward to the standalone
functions with ``self`` supplied as the ``model`` argument. Defining them
via this factory at class-body time - rather than binding
``functools.partial`` instances in ``fit`` - means ``help`` and
``inspect.signature`` report the underlying function (minus ``self``), the
fitted model stays picklable, and subclasses can override cleanly.
"""
@functools.wraps(fn)
def method(self: object, *args, **kwargs) -> object:
return fn(self, *args, **kwargs)
return method
def _scale_block_contributions(blocks: dict[str, np.ndarray], scaling: str) -> dict[str, np.ndarray]:
"""Apply a contribution-plot scaling jointly across every block.
The scalings are those of Miller, Swanson and Heckler (1994). Both are taken
over the blocks together, not one block at a time: a super score pools all
the blocks, so scaling each block separately would make bars from different
blocks incomparable, which is the one thing a multi-block contribution plot
is read for.
"""
if scaling == "none":
return blocks
# A missing cell's contribution is NaN, so both scalings skip NaN: NaN in the max or a
# row's sum left the bars unscaled. A largest of 0, or a row whose total is at most epsqrt
# times the largest row's, has nothing to scale and is left as it is. Both floors are relative,
# so bars in small units are scaled as in units of 1; an absolute epsqrt left them unscaled (#1070).
if scaling == "maximum":
largest = max((float(np.nanmax(np.abs(v))) for v in blocks.values() if np.isfinite(v).any()), default=0.0)
scalar = largest if largest > 0 else 1.0
return {name: values / scalar for name, values in blocks.items()}
if scaling == "within":
totals = sum(np.nansum(np.abs(values), axis=1) for values in blocks.values())
per_row = np.where(totals > epsqrt * np.max(totals, initial=0.0), totals, 1.0).reshape(-1, 1)
return {name: values / per_row for name, values in blocks.items()}
msg = f"scaling must be one of 'none', 'maximum' or 'within', got {scaling!r}."
raise ValueError(msg)
def largest_ssq_column(values: np.ndarray) -> np.ndarray:
"""Return a copy of the column of ``values`` with the largest sum of squares, as an ``(N, 1)`` array.
NIPALS seeds its iteration from this column (#195). Any non-degenerate seed converges to
the same component, but the first column can be constant (all zero after preprocessing),
and then the iteration has nothing to follow. The highest-variance column also needs
fewer iterations. Missing cells are skipped in the sum and set to 0 in the returned column.
Parameters
----------
values : np.ndarray
An ``(N, K)`` array, possibly with NaN.
Returns
-------
np.ndarray
The chosen column, ``(N, 1)``, float, NaN replaced by 0.
"""
ssq = np.nansum(values**2, axis=0)
# Autoscaled columns all have the same sum of squares up to round-off, so a plain argmax picks
# whichever is one ulp larger on this platform. Take the first column within round-off of the
# largest instead: the start, and so the iteration count, is then the same everywhere.
start_col = int(np.flatnonzero(ssq >= ssq.max() * (1 - 1e-9))[0])
return np.nan_to_num(values[:, [start_col]].astype(float))
def _cross_size(x_blocks: list[np.ndarray], y: np.ndarray) -> float:
"""Return ``sqrt(sum_b ||X_b' Y||^2)``, missing cells counted as zero: how much of Y the blocks can still reach."""
y0 = np.nan_to_num(y)
return float(np.sqrt(sum(np.sum((np.nan_to_num(x).T @ y0) ** 2) for x in x_blocks)))
def y_reach(x_blocks: list[np.ndarray], y: np.ndarray) -> tuple[float, float]:
"""Return what :func:`y_is_exhausted` compares against: ``||X'Y||`` and Y's sum of squares, before any component."""
return _cross_size(x_blocks, y), float(np.nansum(y**2))
def check_x_explains_y(x_blocks: list[np.ndarray], y: np.ndarray) -> None:
"""Refuse a fit in which Y is uncorrelated with every column of X: it has no first component (#1144).
``X'Y`` is judged against the largest it can be, ``||X|| ||Y||``, so the data's units do not
matter. With ``X'Y`` about zero, no direction of X covaries with Y, and the first weight is 0 / 0:
PLS and OPLS raised a ``LinAlgError`` that blamed collinear columns, and MBPLS fitted every weight
0 and predicted NaN. A component after the first is fitted along X instead
(:func:`y_is_exhausted`), since the ones before it carry the model. A block with no sum of squares
at all is left to the variance checks, which name it. Missing cells count as zero.
Raises
------
NotEnoughVarianceError
If Y is uncorrelated with every column of X.
"""
x_size = np.sqrt(sum(float(np.nansum(x**2)) for x in x_blocks))
if _cross_size(x_blocks, y) < epsqrt * x_size * np.sqrt(float(np.nansum(y**2))):
raise NotEnoughVarianceError("There is nothing in X to explain Y by: Y is uncorrelated with every column of X.")
def y_is_exhausted(
estimator: str, a: int, x_blocks: list[np.ndarray], y: np.ndarray, initial: tuple[float, float]
) -> bool:
"""Return whether component ``a + 1`` has nothing left in Y to explain, warning when it has (#913).
When what is left of Y is uncorrelated with what is left of X (``X'Y`` about zero
after deflation), the next component has no direction to find: NIPALS then changes
sign at every iteration and reaches ``max_iter``, and its warning points at the
iteration settings rather than at the number of components. This happens in designed
experiments, where one component of orthogonal coded factors already gives the
least-squares fit. ``X'Y`` is judged against its size before the first component.
The caller then fits that component, and every later one, along the largest
direction left in X (NIPALS with X as its own target, which converges to the leading
principal component of the deflated X). Its Y-loading is 0 (#1074), so the
predictions, R²Y and a cross-validated PRESS stay as they were at ``a`` components,
every fitted frame keeps the number of components asked for, and ``P'W`` stays
invertible. The caller asks once: deflating X along a direction of X keeps ``X'Y``
about zero.
Parameters
----------
estimator : str
Name used in the warning.
a : int
Number of components fitted so far.
x_blocks : list[np.ndarray]
The deflated X (one array, or one per block).
y : np.ndarray
The deflated Y.
initial : tuple[float, float]
:func:`y_reach` of the data before the first component.
Returns
-------
bool
True when the fit should stop at ``a`` components.
Raises
------
NotEnoughVarianceError
Before the first component, if Y is uncorrelated with every column of X
(:func:`check_x_explains_y`): there is then no component for the later ones to follow.
"""
initial_cross, initial_ssy = initial
if a == 0:
check_x_explains_y(x_blocks, y)
return False
if _cross_size(x_blocks, y) > epsqrt * initial_cross:
return False
r2y = 1.0 - float(np.nansum(y**2)) / initial_ssy if initial_ssy > 0 else 0.0
warnings.warn(
f"{estimator}: from component {a + 1} on there is nothing left in Y that X can explain: what is left of Y "
f"is uncorrelated with what is left of X (R²Y stays at {r2y:.4g}). The data support {a} component(s); "
"each later one takes the largest direction left in X and leaves the predictions unchanged.",
SpecificationWarning,
stacklevel=user_stacklevel(),
)
return True
def informative_widths(blocks: dict[str, np.ndarray]) -> dict[str, int]:
"""Return, per preprocessed block, the number of columns whose values vary, at least 1 (#913).
A multiblock model divides each block by the square root of its width, so that blocks
of unequal width carry equal weight. A constant column carries no information: centring
and scaling turn it into a column of zeros. Counted in the width, it lowered the
weight of its block, and the fit changed although the data had not. A column with
fewer than two observed values is constant too.
Parameters
----------
blocks : dict[str, np.ndarray]
Each block after centring and scaling, ``(n_samples, K_b)``, NaN where missing.
Returns
-------
dict[str, int]
The count of varying columns per block, at least 1 so a block of constants
divides by 1.
"""
out: dict[str, int] = {}
for name, values in blocks.items():
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning) # a column with no observed value has no spread
spread = np.nanstd(values, axis=0)
out[name] = max(1, int(np.count_nonzero(spread > epsqrt)))
return out
def text_table(header: list[str], rows: list[list[str]]) -> list[str]:
"""Lay out a table as lines of text: each column right-aligned to its widest cell, header included (#1083).
Parameters
----------
header : list[str]
The column headers.
rows : list[list[str]]
The formatted cells of each row, one per column.
Returns
-------
list[str]
The header, a rule of dashes as long as it, then one line per row, each indented by two
spaces, with the cells separated by ``" | "``.
"""
widths = [max(len(cell) for cell in column) for column in zip(header, *rows, strict=True)]
lines = [
" | ".join(cell.rjust(width) for cell, width in zip(line, widths, strict=True)) for line in (header, *rows)
]
return [f" {line}" for line in [lines[0], "-" * len(lines[0]), *lines[1:]]]
def _rows_with_data(values: np.ndarray) -> np.ndarray:
"""Report, for each row, whether it has at least one observed cell."""
return np.any(~np.isnan(values), axis=1)
def unobserved_columns(block: pd.DataFrame, label: str) -> np.ndarray:
"""Return which columns of a training block have no observed value, and warn, naming them, if any do.
Such a column carries no information, so every estimator of the family leaves it out of
the model: it gets a zero loading or weight, and a value a new row has in it adds nothing
to that row's scores or SPE. PCA, PLS, MBPCA and MBPLS used to differ here, from a silent
fit to a refusal (#1098); TPLS's wording (#982) is now the family's.
Parameters
----------
block : pd.DataFrame
The training block, NaN where missing.
label : str
Where the block sits, as ``'X'``, ``'Y'``, ``'X["name"]'`` or ``'F["group"]'``.
Returns
-------
np.ndarray of bool, shape (n_features,)
True for a column with no observed value.
"""
unobserved = block.isna().to_numpy().all(axis=0)
if unobserved.any():
named = ", ".join(repr(col) for col in block.columns[unobserved])
warnings.warn( # a UserWarning, warn's default category
f"Block {label} has columns with no observed values, which carry no information and are "
f"excluded from the model: {named}. Drop them, or check the data for a sensor that was never logged.",
stacklevel=user_stacklevel(),
)
return unobserved
def _equal_values(block: pd.DataFrame) -> np.ndarray:
"""Return, per column, whether its observed values are all equal, as :class:`MCUVScaler` reads a column.
True also for a column with one observed value, and so for every column of a block of one row.
False for a column with no observed value: its NaN maximum and minimum are not equal.
"""
return (block.max() == block.min()).to_numpy()
def check_x_varies(blocks: list[pd.DataFrame], *, explains_y: bool) -> None:
"""Refuse an X in which nothing varies: every column of every block is constant or has no observed value (#1144).
Such an X tells no row from another, so there is no direction for a component to take. MBPLS
refused it; PLS and OPLS refused it at their first component, with another message; PCA and
MBPCA fitted it, with scores that told no row apart. The family now refuses it in MBPLS's
words, before any warning names its columns. A block with no variance beside one that varies
is fitted, with zero weight, and :func:`constant_columns` names it.
Parameters
----------
blocks : list[pd.DataFrame]
The blocks of X as the caller gave them, NaN where missing.
explains_y : bool
Whether X explains a Y, as in PLS, or is described by its components, as in PCA.
Raises
------
NotEnoughVarianceError
If no column of any block varies.
"""
if any((~_equal_values(block) & block.notna().to_numpy().any(axis=0)).any() for block in blocks):
return
where = "every column of every block" if len(blocks) > 1 else "every column"
so = "there is nothing to explain Y by" if explains_y else "there is nothing for a component to describe"
raise NotEnoughVarianceError(f"There is no variance in X: {where} is constant or has no observed value, so {so}.")
def check_y_varies(y: np.ndarray, *, label: str = "Y", explained_by: str = "X") -> None:
"""Refuse a response with nothing to explain: every column is constant or has no observed value (#1149).
Every weight of such a fit is 0, so a new row would have no score. PLS refused it at its
first component, as a Y with no variance *left*; MBPLS, OPLS and TPLS refused it in these
words, which are now the family's. ``y`` is the response as the fit centres it, where such
a column is 0 throughout, and a missing cell counts as 0. A model that does not centre
refuses a Y of zeros the same way.
Parameters
----------
y : np.ndarray of shape (n_samples, n_targets)
The response block over the rows the fit uses, as the fit centres and scales it.
label : str
Where the block sits, as ``'Y'`` or ``'Y["name"]'``.
explained_by : str
The blocks that explain it, as ``'X'`` or ``'F, D and Z'``.
Raises
------
NotEnoughVarianceError
If no column of ``y`` varies.
"""
if not np.any(np.nan_to_num(y)):
raise NotEnoughVarianceError(
f"There is no variance in {label}: every response is constant or has no observed value, so there is "
f"nothing for {explained_by} to explain."
)
def constant_columns(block: pd.DataFrame, label: str, estimator: object) -> np.ndarray:
"""Return which columns of a training block have no variance, and warn, naming them, if any do (#1114).
A column whose observed values are all equal, or which has one, carries no information about
how the rows differ. TPLS named it; PCA, PLS, MBPCA and MBPLS fitted it without a word. Its
wording is now the family's. A column with no observed value is not counted here:
:func:`unobserved_columns` names it, for its own cause. A block of one row has no column that
varies.
A block in which no column varies, beside blocks that do, gets zero weight in the model: its
scores, weights and contributions are 0. It is named as a block, in one sentence, rather than
column by column (#1144). An X in which nothing varies at all is refused before this is called
(:func:`check_x_varies`), and so is a Y.
Only a fit the caller asked for warns. A fit that another routine of this package runs, on a
cross-validation fold or on a matrix it derived (ASCA's effect matrices, FPCA's basis
coefficients), would name a column the caller never gave it, or one that varies in the
caller's data.
Parameters
----------
block : pd.DataFrame
The training block, NaN where missing.
label : str
Where the block sits, as ``'X'``, ``'Y'``, ``'X["name"]'`` or ``'Z["group"]'``.
estimator : object
The estimator being fitted.
Returns
-------
np.ndarray of bool, shape (n_features,)
True for a column with no variance.
"""
constant = _equal_values(block)
if constant.any() and _fit_asked_for_by_the_caller(estimator):
if (constant | block.isna().to_numpy().all(axis=0)).all():
message = (
f"Block {label} has no variance: every column is constant or has no observed value, so it gets "
"zero weight in the model. Drop the block, or check the data for a stuck sensor."
)
else:
named = ", ".join(repr(col) for col in block.columns[constant])
message = (
f"Block {label} has columns with no variance, which carry no information: {named}. Drop them, "
"or check the data for a stuck sensor."
)
warnings.warn(message, stacklevel=user_stacklevel()) # a UserWarning, warn's default category
return constant
def _fit_asked_for_by_the_caller(estimator: object) -> bool:
"""Return whether the caller's code, not another routine of this package, is fitting ``estimator``.
The outermost frame of this package under the caller's is the routine the caller called. It
is a method of ``estimator`` itself (its ``fit``, past an inherited ``fit`` and scikit-learn's
wrapper) when the caller fits it; it is a cross-validation, a resampling or another estimator
when that routine fits it on the caller's behalf.
"""
frame = inspect.currentframe()
outermost = None
while frame is not None:
package = str(frame.f_globals.get("__name__", "")).partition(".")[0]
if package == "process_improve":
outermost = frame
elif package not in ("sklearn", "mutmut"): # scikit-learn's wrapper of fit; mutmut's trampoline
break
frame = frame.f_back
return outermost is not None and outermost.f_locals.get("self") is estimator
def check_response_block(Y: pd.DataFrame, n_samples: int) -> None:
"""Refuse a response block X's rows cannot pair with, or one with an infinite value, as X is refused (#1114).
``validate_data`` checks X alone. The next to meet Y was a scaler, which reported an infinite
response as "Input X contains infinity", and PRM's robust centring, which met a Y of another
length as an ``IndexError``.
Parameters
----------
Y : pd.DataFrame
The response block, NaN where missing.
n_samples : int
The number of rows of X.
Raises
------
ValueError
If ``Y`` has a row count other than ``n_samples``, or an infinite value.
"""
if Y.shape[0] != n_samples:
raise ValueError(
f"The X and Y arrays must have the same number of rows: X has {n_samples} and Y has {Y.shape[0]}."
)
check_array(Y, dtype="numeric", ensure_all_finite="allow-nan", ensure_min_features=0, input_name="Y")
[docs]
class BlockSet(dict):
"""A ``dict[str, pd.DataFrame]`` of equal-height blocks that can also be sliced by row (#193).
:meth:`MBPCA.fit <process_improve.multivariate.methods.MBPCA.fit>` and
:meth:`MBPLS.fit <process_improve.multivariate.methods.MBPLS.fit>` take a
plain ``dict[str, pd.DataFrame]``, which is convenient to build and
impossible to resample: a dict has no notion of "row 7 of every block". Any
resampling or cross-validation pass needs exactly that. ``Resampler``, for
instance, asks its data only for ``len(x)`` and ``x[indices]``.
TPLS already has :class:`~process_improve.multivariate.methods.DataFrameDict`
for this, but it is hardwired to the ``Z``/``F``/``Y`` block names and to a
nested ``dict[str, dict[str, DataFrame]]`` layout, so the multi-block models
could not borrow it. ``BlockSet`` is the flat equivalent: a real ``dict``
subclass, so anything that already accepts the plain dict keeps working,
plus row indexing.
.. warning::
``len(blocks)`` is the number of **rows**, not the number of blocks. That
is surprising for a dict, and it is deliberate: it is the convention
``DataFrameDict`` already set, and it is what the resampling code means
by the length of a dataset. Use ``len(blocks.keys())`` to count blocks.
Parameters
----------
blocks : dict[str, pd.DataFrame]
One entry per block. Every block must be a DataFrame with the same
number of rows; widths may differ.
Raises
------
ValueError
If ``blocks`` is empty, or the blocks disagree on their row count.
TypeError
If any value is not a DataFrame.
Examples
--------
>>> blocks = BlockSet({"a": df_a, "b": df_b}) # doctest: +SKIP
>>> len(blocks) # rows, not blocks # doctest: +SKIP
40
>>> blocks[[0, 1, 2]].keys() # a 3-row BlockSet # doctest: +SKIP
dict_keys(['a', 'b'])
"""
def __init__(self, blocks: dict[str, pd.DataFrame]):
if not isinstance(blocks, dict):
raise TypeError(f"blocks must be a dict of DataFrames, one per block; got {type(blocks).__name__}.")
if not blocks:
raise ValueError("At least one block is required.")
# Two passes, so the "is it a frame?" message lives in one place: the row
# count has to come from a block already known to be a DataFrame, and
# checking the first one separately would mean writing that message twice.
for name, block in blocks.items():
if not isinstance(block, pd.DataFrame):
raise TypeError(f"Block {name!r} must be a pandas DataFrame; got {type(block).__name__}.")
first_name, first = next(iter(blocks.items()))
n_samples = first.shape[0]
for name, block in blocks.items():
if block.shape[0] != n_samples:
raise ValueError(
f"Every block must have the same number of rows ({n_samples}, from block {first_name!r}). "
f"Block {name!r} has {block.shape[0]}."
)
super().__init__(blocks)
self.n_samples = int(n_samples)
self.shape = (self.n_samples, len(blocks))
def __len__(self) -> int:
"""Return the number of rows; see the warning in the class docstring."""
return self.n_samples
def __getitem__(self, lookup: str | int | list | np.ndarray) -> pd.DataFrame | BlockSet:
"""Look a block up by name, or slice every block to the same rows."""
if isinstance(lookup, str):
return cast("pd.DataFrame", super().__getitem__(lookup))
return BlockSet({name: _row_slice(block, lookup, name) for name, block in self.items()})
def __eq__(self, other: object) -> bool:
"""Value-based equality over the held blocks.
The inherited ``dict.__eq__`` compares the values with ``==``, which for
two DataFrames returns an element-wise frame; Python then asks that frame
for its truth value and pandas raises ``ValueError``. Equality therefore
appeared to work only when the two operands were the *same object* (and
the identity short-circuit fired) and blew up otherwise. Comparing with
:meth:`pandas.DataFrame.equals` makes the answer depend on the content,
as it must for a class that also carries ``n_samples`` and ``shape``
(CodeQL ``py/missing-equals``).
Any mapping is accepted on the other side, not just a ``BlockSet``, so
that comparing against the plain ``dict[str, DataFrame]`` the caller
started from answers instead of raising. That stays symmetric: Python
tries the subclass's ``__eq__`` first, so ``plain_dict == block_set``
reaches this method too.
"""
if self is other:
return True
if not isinstance(other, dict):
return NotImplemented
return self.keys() == other.keys() and all(_block_equal(block, other[name]) for name, block in self.items())
def __ne__(self, other: object) -> bool:
"""Negation of :meth:`__eq__`.
Defined explicitly because the C-level ``dict.__ne__`` would otherwise
bypass the Python ``__eq__`` above and compare the raw dict values again.
"""
result = self.__eq__(other)
return result if result is NotImplemented else not result
# Blocks are mutable frames, so a ``BlockSet`` is unhashable just as ``dict``
# is; make that explicit now that ``__eq__`` is defined.
__hash__ = None # type: ignore[assignment] # reason: intentionally unhashable, mirrors dict
def _block_equal(mine: pd.DataFrame, theirs: object) -> bool:
"""Compare one block against the other side's value for the same name.
``theirs`` is a DataFrame whenever the other side is a :class:`BlockSet`, but
an arbitrary object when it is a plain dict, and ``mine == theirs`` would then
hand pandas' element-wise result to :func:`bool`.
"""
return isinstance(theirs, pd.DataFrame) and mine.equals(theirs)
def _row_slice(block: pd.DataFrame, lookup: int | list | np.ndarray, name: str) -> pd.DataFrame:
"""Take rows from one block, keeping the result two-dimensional.
Every branch normalises the lookup to a *list* of row positions and the
function has a single exit, so a scalar can never leak out as a Series: one
row must still arrive at ``fit`` as a one-row frame.
"""
if isinstance(lookup, int | np.integer):
rows = [int(lookup)]
elif isinstance(lookup, np.ndarray):
rows = lookup.tolist()
elif isinstance(lookup, list):
rows = [int(index) for index in lookup]
else:
raise TypeError(
f"Row lookup must be an int, a list of ints, or an ndarray; "
f"got {type(lookup).__name__} while slicing block {name!r}."
)
return block.iloc[rows]
