diff --git a/CHANGE_LOG.txt b/CHANGE_LOG.txt index 3850dc9..ac25c3b 100644 --- a/CHANGE_LOG.txt +++ b/CHANGE_LOG.txt @@ -1,7 +1,12 @@ +Change Log 0.2.3 +----------------- + +- Added GMNet, based on the original implementation https://github.com/AICGijon/gmnet + Change Log 0.2.2 ----------------- -- Added HistNetQ, based on the original implementation https://github.com/pglez84/histnetq +- Added HistNetQ, based on the original implementation https://github.com/pglez82/histnetq - Minor fixes diff --git a/README.md b/README.md index 212d6c8..fdb80d8 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # QuaPy -## version 0.2.2 +## version 0.2.3 QuaPy is an open source framework for quantification (a.k.a. supervised prevalence estimation, or learning to quantify) written in Python. @@ -15,7 +15,7 @@ for facilitating the analysis and interpretation of the experimental results. ### Last updates: -* Version 0.2.2 is released! major changes can be consulted [here](CHANGE_LOG.txt). +* Version 0.2.3 is released! major changes can be consulted [here](CHANGE_LOG.txt). * The developer API documentation is available [here](https://hlt-isti.github.io/QuaPy/index.html) ### Installation diff --git a/TODO.txt b/TODO.txt index 17604ab..fd720cc 100644 --- a/TODO.txt +++ b/TODO.txt @@ -19,4 +19,8 @@ scale each value by per-class thresholds, i.e., [0.33*0.1, 0.33*1, 0.33*1]/sum. - [TODO] add Friedman's method and DeBias - [TODO] check ignore warning stuff check https://docs.python.org/3/library/warnings.html#temporarily-suppressing-warnings -- [TODO] nmd and md are not selectable from qp.evaluation.evaluate as a string \ No newline at end of file +- [TODO] nmd and md are not selectable from qp.evaluation.evaluate as a string +- [TODO] add https://www.kaggle.com/datasets/Cornell-University/arxiv dataset? in the paper "Online Adaptation to + Label Distribution Shift" they say it provides a natural label shift over time +- [TODO] add Bhattacharyya distance (https://en.wikipedia.org/wiki/Bhattacharyya_distance) + (which is actually not a proper distance) \ No newline at end of file diff --git a/docs/source/manuals.rst b/docs/source/manuals.rst index 97e0628..f81ab7c 100644 --- a/docs/source/manuals.rst +++ b/docs/source/manuals.rst @@ -7,6 +7,7 @@ Manuals manuals/datasets manuals/evaluation + manuals/label-shift-adaptation manuals/methods manuals/model-selection manuals/plotting diff --git a/docs/source/manuals/label-shift-adaptation.md b/docs/source/manuals/label-shift-adaptation.md new file mode 100644 index 0000000..f71994f --- /dev/null +++ b/docs/source/manuals/label-shift-adaptation.md @@ -0,0 +1,108 @@ +# Label Shift Adaptation + +Quantification methods estimate the class prevalence values of a test sample, but a prevalence estimate is +not, by itself, a classifier that has been corrected for the shift it describes. In some scenarios you may +want the latter: a classifier whose decisions account for the (estimated) change in class priors between +the training and the test distributions, rather than a single prevalence vector. + +QuaPy provides two building blocks for this: + +* `qp.method.aggregative.ImportanceWeightQuantifier`, an interface implemented by the quantifiers that + natively compute a vector of importance weights as part of estimating the target prevalence. +* `qp.classification.labelshift.LabelShiftedClassifier`, a classifier wrapper that uses any quantifier + (whether or not it implements the interface above) to adapt a classifier to a specific, (potentially) + shifted batch of test instances. + +## The `ImportanceWeightQuantifier` interface + +Some quantifiers (currently, `RLLS`, `BBSEhard`, and `BBSEsoft`; see the +{ref}`Regularized Learning under Label Shift (RLLS) ` +and {ref}`Black Box Shift Estimation (BBSE) ` sections of +the methods manual) work by first estimating a vector of importance weights + +:::{math} +w_y = \frac{Q(y)}{P(y)} +::: + +with :math:`P` and :math:`Q` the training and target class distributions, and only then rescaling the +training prevalence by these weights to obtain the target prevalence estimate: +:math:`\hat{Q}(y) \propto w_y \cdot P(y)`. For these methods, the weight vector is not an afterthought +computed from the prevalence estimate; it is the primitive quantity from which the estimate itself is +derived. + +`ImportanceWeightQuantifier` exposes this quantity directly, in addition to the usual `predict` method +inherited from any quantifier: + +* `get_importance_weights(instances)` returns the weight vector $w$ estimated for a given batch of + (unlabelled) target instances. +* `quantify_and_weigh(instances)` returns a tuple `(prevalence, weights)`, computed from a single pass of + classifier predictions over the instances (avoiding the double classification that calling `predict` and + `get_importance_weights` separately would incur). + +Neither method mutates any internal state of the quantifier, so both are safe to call concurrently, on the +same fitted instance, for different batches of target instances. + +```python +import quapy as qp +from quapy.method.aggregative import BBSEhard +from sklearn.linear_model import LogisticRegression + +train, test = qp.datasets.fetch_UCIBinaryDataset('haberman').train_test + +model = BBSEhard(LogisticRegression(max_iter=2000), val_split=5) +model.fit(*train.Xy) + +prevalence, weights = model.quantify_and_weigh(test.X) +``` + +## Adapting a classifier: `LabelShiftedClassifier` + +Rescaling the classifier's posterior probabilities after the fact (as, e.g., `EMQ` does) is one way of +accounting for an estimated prevalence shift. Another common practice is instead to retrain the classifier +with a reweighted loss, via scikit-learn's `class_weight` constructor parameter; this changes the decision +function of the classifier itself, rather than only rescaling its outputs, which can matter for classifiers +whose decision boundary is not a simple function of the posterior probabilities alone. + +`LabelShiftedClassifier` automates this second strategy. Given a base classifier `h` and any `quapy` +quantifier `q`, it fits `q` on the labelled training data (keeping a copy of it), and, for a given batch of +(unlabelled) target instances, estimates the importance weights: directly, via `get_importance_weights`, if +`q` is an `ImportanceWeightQuantifier`; or otherwise by dividing `q`'s estimated target prevalence by the +training prevalence. It then retrains a *fresh copy* of `h`, with `class_weight` set according to these +weights, on the original training data, and hands off classification of the target instances to this +adapted classifier. Neither the classifier instance `h` passed at construction time, nor `q`'s own internal +classifier (if any), are ever mutated. + +```python +import quapy as qp +from quapy.classification.labelshift import LabelShiftedClassifier +from quapy.method.aggregative import ACC +from sklearn.linear_model import LogisticRegression + +train, test = qp.datasets.fetch_UCIBinaryDataset('haberman').train_test + +h = LogisticRegression(max_iter=2000) +q = ACC(LogisticRegression(max_iter=2000), val_split=5) + +adapter = LabelShiftedClassifier(h, q) +adapter.fit(*train.Xy) + +# label predictions and posterior probabilities from a classifier adapted to test's estimated distribution +y_pred = adapter.predict(test.X) +y_proba = adapter.predict_proba(test.X) + +# or, to obtain the adapted classifier itself +adapted_h = adapter.get_classifier(test.X) +``` + +`h` is required to accept a `class_weight` parameter in its constructor, as `LogisticRegression` and many +other scikit-learn classifiers do; `q` can be any (fitted or unfitted) `quapy` quantifier, and need not be an +`ImportanceWeightQuantifier`. + +Since the adaptation is specific to the batch of target instances given, `predict` and `predict_proba` +retrain `h` on every call. If you plan to issue several predictions against the same (suspected) shifted +distribution, call `get_classifier` once and reuse the classifier instance it returns, rather than calling +`predict`/`predict_proba` repeatedly on the same or related samples. + +A small positive floor, controlled by the `weight_epsilon` constructor parameter (default `1e-4`), is +applied to the estimated weights before they are used as `class_weight`; this keeps a class with training +support from being assigned a weight of exactly `0`, which would silently drop it from the retraining loss. diff --git a/docs/source/manuals/methods.md b/docs/source/manuals/methods.md index b721957..dcab373 100644 --- a/docs/source/manuals/methods.md +++ b/docs/source/manuals/methods.md @@ -323,6 +323,60 @@ model.fit(*train.Xy) estim_prevalence = model.predict(test.X) ``` +RLLS computes its importance weights directly (it implements the `ImportanceWeightQuantifier` interface); see +the {ref}`Label Shift Adaptation manual ` for how to +access these weights, or use them to adapt a classifier itself rather than only estimating prevalence. + +### Black Box Shift Estimation (BBSE) + +`BBSEhard` and `BBSEsoft` are available at `qp.method.aggregative.BBSEhard` and +`qp.method.aggregative.BBSEsoft`, respectively, and implement the Black Box Shift Estimator +proposed in: + +_Lipton, Z., Wang, Y. X., & Smola, A. (2018, July). Detecting and correcting for label shift +with black box predictors. In International conference on machine learning +(pp. 3122-3130). PMLR._ ([link to paper](https://proceedings.mlr.press/v80/lipton18a.html)) + +BBSE is similar in spirit to ACC and PACC in that it exploits the label-shift invariance of +`P(hat{Y}|Y)` to correct for the change in class prevalence between the training and the test +distributions. However, while ACC solves the linear system `q = Mp` (with `M` the matrix of +class-conditional misclassification rates and `p` the sought prevalence vector), BBSE instead +solves `q = Cw` for the importance-weight vector `w`, with `w_i = Q(i)/P(i)` the ratio between +the target and the source class priors, and `C` the joint-distribution matrix +`C_ij = P(hat{Y}=i, Y=j)` estimated on a validation split. The target prevalence estimate is +then recovered as `Q(y) = w_y * P(y)`. + +`BBSEhard` estimates `C` from crisp classifier predictions (i.e., a standard confusion matrix, +normalized to sum to 1), while `BBSEsoft` estimates it from the classifier's posterior +probabilities instead, in the same spirit in which PACC generalizes ACC. + +```python +import quapy as qp +from quapy.method.aggregative import BBSEhard, BBSEsoft +from sklearn.linear_model import LogisticRegression + +train, test = qp.datasets.fetch_UCIBinaryDataset('haberman').train_test + +model = BBSEhard(LogisticRegression(max_iter=2000), val_split=5) +model.fit(*train.Xy) +estim_prevalence = model.predict(test.X) + +# or, using posterior probabilities instead of crisp counts: +model = BBSEsoft(LogisticRegression(max_iter=2000), val_split=5) +model.fit(*train.Xy) +estim_prevalence = model.predict(test.X) +``` + +As with ACC and RLLS, both variants require validation predictions and therefore expect +`val_split` to be set whenever `fit_classifier=True`. They also accept the same `solver` +(`"minimize"`, `"exact-raise"`, `"exact-cc"`) and `norm` (`"clip"`, `"mapsimplex"`, +`"condsoftmax"`) arguments discussed above for ACC/PACC. + +Like RLLS, both `BBSEhard` and `BBSEsoft` implement the `ImportanceWeightQuantifier` interface; see the +{ref}`Label Shift Adaptation manual ` for how to +access these weights directly, or use them to adapt a classifier itself rather than only estimating +prevalence. + ### Distribution Matching Distribution Matching (DM) methods search for the mixture parameter (the sought class prevalence values) @@ -853,7 +907,8 @@ estim_prevalence = model.predict(dataset.test.X) ### HistNetQ -QuaPy offers an implementation of HistNetQ, a deep learning model based on a differentiable +QuaPy offers an implementation of HistNetQ (based on the [original implementation](https://github.com/pglez82/histnetq)), +a deep learning model based on a differentiable histogram representation, presented in: [_Pérez-Mon, O., Moreo, A., del Coz, J.J., & González, P. (2025). @@ -880,6 +935,47 @@ HistNetQ can alternatively be trained directly from a protocol that already prov samples (e.g., when only bag-level prevalence values are available), via the `fit_from_samples` method; see the API documentation for further details. +### GMNet + +QuaPy offers an implementation of GMNet (based on the [original implementation](https://github.com/AICGijon/gmnet)), +a deep learning model that represents each instance of a +bag by its likelihood under one or more learned mixtures of Gaussians, presented in: + +[_Pérez-Mon, O., del Coz, J.J., & González, P. (2026). +Quantification via Gaussian latent space representations. +Neural Networks._](https://www.sciencedirect.com/science/article/pii/S0893608026003473) + +This model requires `torch` and `geotorch` to be installed (`geotorch` is used to keep the Gaussian +layers' covariance matrices positive-definite during training). Like HistNetQ, GMNet is trained +end-to-end on samples ("bags") of known prevalence rather than on individually labeled instances, +and requires no classifier, only an optional feature extraction module. + +```python +import quapy as qp +from quapy.method.meta import GMNet + +dataset = qp.datasets.fetch_UCIBinaryDataset('haberman') + +model = GMNet(bag_size=100, device='cpu') +model.fit(*dataset.training.Xy) +estim_prevalence = model.predict(dataset.test.X) +``` + +GMNet supports stacking multiple "GM branches" (via `n_gm_layers`, `num_gaussians`, and +`gaussian_dimensions`), optionally regularized with a CKA (Centered Kernel Alignment) term that +encourages the branches to learn dissimilar latent representations, set through +`cka_regularization`: + +```python +model = GMNet( + n_gm_layers=2, num_gaussians=(4, 4), gaussian_dimensions=(8, 8), + cka_regularization=0.1, bag_size=100, device='cpu' +) +``` + +Like HistNetQ, GMNet can alternatively be trained directly from a protocol that already provides the +training samples, via the `fit_from_samples` method; see the API documentation for further details. + ## Quantifiers with Uncertainty Quantification diff --git a/docs/source/quapy.classification.rst b/docs/source/quapy.classification.rst index cfc7d9b..c07b72e 100644 --- a/docs/source/quapy.classification.rst +++ b/docs/source/quapy.classification.rst @@ -12,6 +12,14 @@ quapy.classification.calibration module :undoc-members: :show-inheritance: +quapy.classification.labelshift module +--------------------------------------- + +.. automodule:: quapy.classification.labelshift + :members: + :undoc-members: + :show-inheritance: + quapy.classification.methods module ----------------------------------- diff --git a/examples/22.histnet_gmnet_bagtraining.py b/examples/22.histnet_gmnet_bagtraining.py new file mode 100644 index 0000000..985c138 --- /dev/null +++ b/examples/22.histnet_gmnet_bagtraining.py @@ -0,0 +1,68 @@ +import quapy as qp +from quapy.data.datasets import fetch_image_embeddings +from quapy.method.meta import HistNetQ, GMNet +from quapy.protocol import UPP + +# This example showcases HistNetQ and GMNet, two neural quantifiers that -- unlike the +# classify-then-aggregate methods seen in most other examples (e.g., EMQ, PACC) -- are trained +# end-to-end on whole samples ("bags") of known prevalence, instead of on individually labelled +# instances. Both require `torch` to be installed (GMNet additionally requires `geotorch`); see +# quapy.method._neural_bags.BagTrainedQuantifier, the base class they share, and the "Methods" manual +# for further details. +# +# In particular, we exercise `fit_from_samples`, the entry point meant for the situation in which +# individually labelled training instances are NOT available at all, and all one has access to is a +# pool of pre-built bags of known prevalence -- e.g., as released for the LeQua challenges, or as +# collected by a third party. We simulate such a pool here using CIFAR10 (see also example 20), by +# drawing a fixed set of bags from its (otherwise fully instance-labelled) training collection with +# the UPP protocol, and then training HistNetQ and GMNet as if that pool of bags -- and not the +# underlying instance-level labels -- were all we had access to. +# +# The settings below (number of bags, training epochs, patience) are deliberately modest so that the +# example runs in a reasonable time on a CPU; scale them up for serious experimentation. + +if __name__ == '__main__': + + qp.environ['SAMPLE_SIZE'] = 500 + + # The datasets available in quapy do not consist of raw image files, but are instead + # pre-generated embeddings (see the manuals for further information); 'features' corresponds to + # the penultimate-layer representations of a resnet18, which is what HistNetQ and GMNet expect to + # receive, since both operate on continuous instance representations. + print('fetching cifar10 embeddings') + train, test = fetch_image_embeddings(dataset_name='cifar10', embedding='features').train_test + print('training:', train) + print('test:', test) + + # Let us pretend that, instead of having access to the individually labelled instances of `train`, + # we are only given a pool of 1000 pre-built bags of 500 instances each, with heterogeneous class + # prevalence values (sampled uniformly at random from the simplex). This is exactly what + # fit_from_samples expects: an AbstractProtocol that already yields (sample, prevalence) pairs. + # UPP is deterministic by default (fixed random_state=0), so this pool is the same every time it + # is iterated, as befits a fixed, pre-built collection of bags. + print('simulating a pool of 1000 pre-built training bags (as if released by a third party)') + bag_pool = UPP(train, repeats=1000, random_state=0) + + # Neither method requires a classifier, nor (in this case) a feature extraction module: since the + # CIFAR10 embeddings are already vectorized, the default identity module is used, and both + # networks operate directly on the 512-dimensional resnet18 features + models = { + 'HistNetQ': HistNetQ(n_bins=8, bag_size=qp.environ['SAMPLE_SIZE'], train_epochs=50, patience=10, verbose=1), + 'GMNet': GMNet(n_gm_layers=8, num_gaussians=10, gaussian_dimensions=5, hidden_size_fe=(50,), + bag_size=qp.environ['SAMPLE_SIZE'], train_epochs=50, patience=10, verbose=1), + } + + for name, model in models.items(): + print(f'\nfitting {name} from the pre-built bag pool via fit_from_samples ' + f'(no instance-level labels are used)') + model.fit_from_samples(bag_pool) + + # we now evaluate both quantifiers on cifar10's test set, using a fresh artificial-prevalence + # protocol (note that, at prediction time, HistNetQ and GMNet behave just like any other + # quantifier: the bag-based training is an internal detail of these two particular methods) + test_prot = UPP(test, repeats=1000) + + for name, model in models.items(): + report = qp.evaluation.evaluation_report(model, protocol=test_prot, error_metrics=['mae']) + print(f'\n{name}:') + print(report.mean(numeric_only=True)) diff --git a/quapy/__init__.py b/quapy/__init__.py index 4ed6726..de4575c 100644 --- a/quapy/__init__.py +++ b/quapy/__init__.py @@ -17,7 +17,7 @@ try: except ImportError: plot = None -__version__ = '0.2.2' +__version__ = '0.2.3' def _default_cls(): diff --git a/quapy/classification/labelshift.py b/quapy/classification/labelshift.py new file mode 100644 index 0000000..6225f1e --- /dev/null +++ b/quapy/classification/labelshift.py @@ -0,0 +1,110 @@ +import numpy as np +from sklearn.base import BaseEstimator, clone + +import quapy.functional as F + + +class LabelShiftedClassifier(BaseEstimator): + """ + Adapts a classifier to a shifted (unlabelled) target distribution by means of a `quapy` quantifier. + + Given a base classifier `h` and a quantifier `q`, this wrapper: (i) at `fit` time, fits `q` on the labelled + training data and keeps a copy of this data; (ii) at prediction time, for a given batch of (unlabelled) + target instances, estimates the importance weights :math:`w_y=Q(y)/P(y)` that would explain the shift + between the training distribution and the target one (directly, if `q` is an + :class:`quapy.method.aggregative.ImportanceWeightQuantifier`, or otherwise by dividing `q`'s prevalence + estimate by the training prevalence); and (iii) retrains a fresh copy of `h`, reweighted according to `w` + via `class_weight`, on the original training data, before delegating the classification of the target + instances to this adapted classifier. + + Note that a fresh copy of `h` is retrained for every batch of target instances passed to `predict`, + `predict_proba`, or `get_classifier`, since the adaptation is specific to that batch's estimated shift; `q` + and its own internal classifier (if any) are otherwise left untouched. `h` itself is required to accept a + `class_weight` parameter in its constructor (e.g., `sklearn.linear_model.LogisticRegression`). + + :param classifier: a scikit-learn classifier that accepts a `class_weight` parameter in its constructor + :param quantifier: a `quapy` quantifier, used to estimate the target prevalence (or, if it is an + :class:`quapy.method.aggregative.ImportanceWeightQuantifier`, the importance weights directly) + :param weight_epsilon: float, a small positive floor applied to the estimated importance weights (default + 1e-4), so that a class with training support is never assigned a `class_weight` of exactly 0 (which + would silently drop it from the retraining loss) + """ + + def __init__(self, classifier, quantifier, weight_epsilon=1e-4): + self.classifier = classifier + self.quantifier = quantifier + self.weight_epsilon = weight_epsilon + + def fit(self, X, y): + """ + Fits the internal quantifier on the training data, and stores the training data for the subsequent, + per-target-batch, classifier adaptation. + + :param X: array-like of shape `(n_samples, n_features)` with the training instances + :param y: array-like of shape `(n_samples,)` with the class labels + :return: self + """ + self.quantifier.fit(X, y) + self.classes_ = self.quantifier.classes_ + self.train_prevalence_ = F.prevalence_from_labels(y, classes=self.classes_) + self.X_, self.y_ = X, y + return self + + def get_importance_weights(self, X): + """ + Estimates the vector of importance weights :math:`w_y=Q(y)/P(y)` that would explain the shift between + the training distribution and the distribution of the given target instances. This method does not + mutate any internal state, so it is safe to call concurrently for different batches of target + instances. + + :param X: array-like of shape `(n_samples, n_features)` with the (unlabelled) target instances + :return: np.ndarray of shape `(n_classes,)` + """ + from quapy.method.aggregative import ImportanceWeightQuantifier + if isinstance(self.quantifier, ImportanceWeightQuantifier): + weights = self.quantifier.get_importance_weights(X) + else: + test_prevalence = self.quantifier.predict(X) + weights = test_prevalence / np.maximum(self.train_prevalence_, self.weight_epsilon) + return np.clip(weights, self.weight_epsilon, None) + + def get_classifier(self, X): + """ + Returns a fresh classifier instance, sharing `classifier`'s hyperparameters, retrained on the original + training data with `class_weight` set according to the importance weights estimated for the given + target instances. The classifier passed at construction time, and the one used internally by the + quantifier, are left untouched. + + :param X: array-like of shape `(n_samples, n_features)` with the (unlabelled) target instances + :return: a fitted scikit-learn classifier + """ + weights = self.get_importance_weights(X) + class_weight = dict(zip(self.classes_, weights)) + adapted = clone(self.classifier) + adapted.set_params(class_weight=class_weight) + adapted.fit(self.X_, self.y_) + return adapted + + def predict(self, X): + """ + Adapts the classifier to the target instances in `X`, and returns its label predictions for `X`. + Adaptation implies retraining. If you intend to generate inferences over many samples from the same + suspected distribution, then call the `get_classifier` method once, and use the returned classifier + over such samples. + + :param X: array-like of shape `(n_samples, n_features)` with the (unlabelled) target instances + :return: array-like of shape `(n_samples,)` with the class label predictions + """ + return self.get_classifier(X).predict(X) + + def predict_proba(self, X): + """ + Adapts the classifier to the target instances in `X`, and returns its posterior probabilities for `X`. + Adaptation implies retraining. If you intend to generate inferences over many samples from the same + suspected distribution, then call the `get_classifier` method once, and use the returned classifier + over such samples. + + :param X: array-like of shape `(n_samples, n_features)` with the (unlabelled) target instances + :return: array-like of shape `(n_samples, n_classes)` with posterior probabilities + """ + return self.get_classifier(X).predict_proba(X) diff --git a/quapy/functional.py b/quapy/functional.py index 392309e..84504ee 100644 --- a/quapy/functional.py +++ b/quapy/functional.py @@ -462,7 +462,7 @@ def argmin_prevalence(loss: Callable, raise NotImplementedError() -def optim_minimize(loss: Callable, n_classes: int, return_loss=False): +def optim_minimize(loss: Callable, n_classes: int, x0='uniform', bounds='simplex', constraints='simplex', return_loss=False): """ Searches for the optimal prevalence values, i.e., an `n_classes`-dimensional vector of the (`n_classes`-1)-simplex that yields the smallest lost. This optimization is carried out by means of a constrained search using scipy's @@ -470,19 +470,32 @@ def optim_minimize(loss: Callable, n_classes: int, return_loss=False): :param loss: (callable) the function to minimize :param n_classes: (int) the number of classes, i.e., the dimensionality of the prevalence vector + :param x0: initial solution; if the string 'uniform' is passed (default) then the initial solution is the + uniform distribution; otherwise, a valid object must be provided + :param bounds: the bounds of the search space; if the string `simplex` is passed (default) then the bounds + of a simplex of appropriate dimension is instantiated; otherwise, a valid tuple must be provided + :param constraints: the constraints of valid solutions; if the string `simplex` is passed (default) then the + constraint that a point must lie on the simplex is assumed; otherwise, valid constrains must be provided + (see scipy.optimize) :param return_loss: bool, if True, returns also the value of the loss (default is False). :return: (ndarray) the best prevalence vector found or a tuple which also contains the value of the loss if return_loss=True """ from scipy import optimize - # the initial point is set as the uniform distribution - uniform_distribution = uniform_prevalence(n_classes=n_classes) + if isinstance(x0, str) and x0=='uniform': + # the initial point is set as the uniform distribution + x0 = uniform_prevalence(n_classes=n_classes) - # solutions are bounded to those contained in the unit-simplex - bounds = tuple((0, 1) for _ in range(n_classes)) # values in [0,1] - constraints = ({'type': 'eq', 'fun': lambda x: 1 - sum(x)}) # values summing up to 1 - r = optimize.minimize(loss, x0=uniform_distribution, method='SLSQP', bounds=bounds, constraints=constraints) + if isinstance(bounds, str) and bounds=='simplex': + # solutions are bounded to those contained in the unit-simplex + bounds = tuple((0, 1) for _ in range(n_classes)) # values in [0,1] + + if isinstance(constraints, str) and constraints=='simplex': + # the point must lie on the simplex + constraints = ({'type': 'eq', 'fun': lambda x: 1 - sum(x)}) # values summing up to 1 + + r = optimize.minimize(loss, x0=x0, method='SLSQP', bounds=bounds, constraints=constraints) if return_loss: return r.x, r.fun @@ -623,84 +636,6 @@ def solve_adjustment_binary(prevalence_estim: ArrayLike, tpr: float, fpr: float, return adjusted -def solve_adjustment( - class_conditional_rates: np.ndarray, - unadjusted_counts: np.ndarray, - method: Literal["inversion", "invariant-ratio"], - solver: Literal["exact", "minimize", "exact-raise", "exact-cc"]) -> np.ndarray: - """ - Function that tries to solve for :math:`p` the equation :math:`q = M p`, where :math:`q` is the vector of - `unadjusted counts` (as estimated, e.g., via classify and count) with :math:`q_i` an estimate of - :math:`P(\hat{Y}=y_i)`, and where :math:`M` is the matrix of `class-conditional rates` with :math:`M_{ij}` an - estimate of :math:`P(\hat{Y}=y_i|Y=y_j)`. - - :param class_conditional_rates: array of shape `(n_classes, n_classes,)` with entry `(i,j)` being the estimate - of :math:`P(\hat{Y}=y_i|Y=y_j)`, that is, the probability that an instance that belongs to class :math:`y_j` - ends up being classified as belonging to class :math:`y_i` - - :param unadjusted_counts: array of shape `(n_classes,)` containing the unadjusted prevalence values (e.g., as - estimated by CC or PCC) - - :param str method: indicates the adjustment method to be used. Valid options are: - - * `inversion`: tries to solve the equation :math:`q = M p` as :math:`p = M^{-1} q` where - :math:`M^{-1}` is the matrix inversion of :math:`M`. This inversion may not exist in - degenerated cases. - * `invariant-ratio`: invariant ratio estimator of `Vaz et al. 2018 `_, - which replaces the last equation in :math:`M` with the normalization condition (i.e., that the sum of - all prevalence values must equal 1). - - :param str solver: the method to use for solving the system of linear equations. Valid options are: - - * `exact-raise`: tries to solve the system using matrix inversion. Raises an error if the matrix has rank - strictly lower than `n_classes`. - * `exact-cc`: if the matrix is not full rank, returns :math:`q` (i.e., the unadjusted counts) as the estimates - * `exact`: deprecated, defaults to 'exact-cc' (will be removed in future versions) - * `minimize`: minimizes a loss, so the solution always exists - """ - if solver == "exact": - warnings.warn( - "The 'exact' solver is deprecated. Use 'exact-raise' or 'exact-cc'", DeprecationWarning, stacklevel=2) - solver = "exact-cc" - - A = np.asarray(class_conditional_rates, dtype=float) - B = np.asarray(unadjusted_counts, dtype=float) - - if method == "inversion": - pass # We leave A and B unchanged - elif method == "invariant-ratio": - # Change the last equation to replace it with the normalization condition; - # copy first so this does not mutate the caller's arrays (np.asarray above - # returns the same object, not a copy, when the input is already float64) - A = A.copy() - B = B.copy() - A[-1, :] = 1.0 - B[-1] = 1.0 - else: - raise ValueError(f"unknown {method=}") - - if solver == "minimize": - def loss(prev): - return np.linalg.norm(A @ prev - B) - return optim_minimize(loss, n_classes=A.shape[0]) - elif solver in ["exact-raise", "exact-cc"]: - # Solvers based on matrix inversion, so we use try/except block - try: - return np.linalg.solve(A, B) - except np.linalg.LinAlgError: - # The matrix is not invertible. - # Depending on the solver, we either raise an error - # or return the classifier predictions without adjustment - if solver == "exact-raise": - raise - elif solver == "exact-cc": - return unadjusted_counts - else: - raise ValueError(f"Solver {solver} not known.") - else: - raise ValueError(f'unknown {solver=}') - - # ------------------------------------------------------------------------------------------ # Transformations from Compositional analysis # ------------------------------------------------------------------------------------------ diff --git a/quapy/method/__init__.py b/quapy/method/__init__.py index 0bcab12..d295382 100644 --- a/quapy/method/__init__.py +++ b/quapy/method/__init__.py @@ -15,6 +15,8 @@ AGGREGATIVE_METHODS = { aggregative.ACC, aggregative.PCC, aggregative.PACC, + aggregative.BBSEhard, + aggregative.BBSEsoft, aggregative.RLLS, aggregative.EMQ, aggregative.HDy, @@ -30,6 +32,7 @@ AGGREGATIVE_METHODS = { aggregative.KDEyML, aggregative.KDEyCS, aggregative.KDEyHD, + aggregative.LEIP, # aggregative.OneVsAllAggregative, confidence.BayesianCC, _bayesian.BayesianKDEy, @@ -54,12 +57,15 @@ MULTICLASS_METHODS = { aggregative.ACC, aggregative.PCC, aggregative.PACC, + aggregative.BBSEhard, + aggregative.BBSEsoft, aggregative.RLLS, aggregative.EMQ, aggregative.EDy, aggregative.KDEyML, aggregative.KDEyCS, aggregative.KDEyHD, + aggregative.LEIP, confidence.BayesianCC, _bayesian.BayesianKDEy, _bayesian.BayesianMAPLS, @@ -75,7 +81,8 @@ NON_AGGREGATIVE_METHODS = { META_METHODS = { meta.Ensemble, meta.QuaNet, - meta.HistNetQ + meta.HistNetQ, + meta.GMNet } QUANTIFICATION_METHODS = AGGREGATIVE_METHODS | NON_AGGREGATIVE_METHODS | META_METHODS diff --git a/quapy/method/_gmnet.py b/quapy/method/_gmnet.py new file mode 100644 index 0000000..6222705 --- /dev/null +++ b/quapy/method/_gmnet.py @@ -0,0 +1,328 @@ +""" +GMNet implementation. + +Ported from the reference implementation at https://github.com/AICGijon/gmnet (the `GMNet`/ +`DLQuantification` classes in that repo), adapted to QuaPy's own protocol-based sample generation +(replacing that repo's custom, `quantificationlib`-backed bag generators), and reusing the shared +bag-based training loop already factored out for :class:`quapy.method._histnet.HistNetQ` (see +:class:`quapy.method._neural_bags.BagTrainedQuantifier`). + +The overall architecture is: one or more "GM branches" -- each a small per-branch feature extractor +followed by a layer of Gaussian likelihoods (a :class:`_GMLayer`) evaluated at every instance of a bag +-- concatenated and mean-pooled over the bag, followed by the shared quantification MLP head. Like +HistNetQ (and QuaNet), GMNet is trained end-to-end by minimizing a quantification loss over samples +("bags") of known prevalence, rather than over individually labeled instances. + +Two deliberate deviations from the reference implementation, both required for the model to satisfy +QuaPy's `predict(X)` contract (i.e., to be usable on a real test collection of arbitrary size, as +opposed to only on bags resampled at the fixed `bag_size` used for training): + +* the original `GMNet_Module` reshapes each branch's per-instance likelihoods around a *fixed*, + constructor-time `bag_size` (via `torch.nn.Unflatten(0, (-1, bag_size))`), which only works when + every forward pass is fed bags of exactly that size. Here, the reshape is instead computed from the + actual input shape at forward time (see :class:`_GMBranch`), which is equivalent when the bag size + matches but also supports bags (or, at prediction time, whole test samples) of any other size. +* the forward hooks used by the original code to capture each branch's pre-Gaussian latent activations + (for the CKA regularization term) are replaced by simply storing that activation as an attribute + during `forward` (see :attr:`_GMBranch.latent_activation`), since branches are now implemented with a + plain `forward` method rather than an opaque `torch.nn.Sequential`. +""" +import numpy as np +import scipy.spatial.distance +import torch +import torch.nn as nn + +import geotorch + +from quapy.method._neural_bags import BagTrainedQuantifier +from quapy.protocol import UPP + + +def _cka(latent_activations): + """Feature-space linear CKA (Centered Kernel Alignment), averaged over every pair of latent + activations, following the `CKARegularization` class in the reference implementation. Used to + encourage the Gaussian components learned by different GM branches to capture complementary + (dissimilar) aspects of the instances. + + :param latent_activations: a list of tensors, one per GM branch, all of shape (n_instances, dim_i) + (dim_i may differ across branches) + """ + cka_sum = 0. + n_pairs = 0 + for i in range(len(latent_activations)): + for j in range(i + 1, len(latent_activations)): + x = latent_activations[i] + y = latent_activations[j] + x = x - torch.mean(x, dim=0, keepdim=True) + y = y - torch.mean(y, dim=0, keepdim=True) + dot_product_similarity = torch.norm(torch.matmul(x.t(), y)) ** 2 + normalization_x = torch.norm(torch.matmul(x.t(), x)) + normalization_y = torch.norm(torch.matmul(y.t(), y)) + cka_sum = cka_sum + dot_product_similarity / (normalization_x * normalization_y) + n_pairs += 1 + return cka_sum / n_pairs + + +class _GMLayer(nn.Module): + """A layer of `num_gaussians` (unnormalized) Gaussian likelihoods, evaluated at every instance of + a bag. `centers` and `covariance` are learned; `covariance` is constrained to stay positive-definite + throughout training via `geotorch.positive_definite`. + """ + + def __init__(self, n_features, num_gaussians): + super().__init__() + self.n_features = n_features + self.num_gaussians = num_gaussians + self.centers = nn.Parameter(torch.rand(num_gaussians, n_features)) + self.covariance = nn.Parameter(torch.eye(n_features).repeat(num_gaussians, 1, 1)) + geotorch.positive_definite(self, "covariance") + + # initialize the centers' covariance from the (squared, halved) nearest-neighbor distance + # between the randomly initialized centers, so that gaussians start with a sensible spread + centers = self.centers.detach().cpu().numpy() + distances = scipy.spatial.distance.cdist(centers, centers) + np.fill_diagonal(distances, np.inf) + cov = (np.mean(np.min(distances, axis=1)) / 2) ** 2 + self.covariance = torch.eye(n_features).repeat(num_gaussians, 1, 1) * cov + + def forward(self, x): + # x: (batch_size, bag_size, n_features) + centers = self.centers.unsqueeze(0).unsqueeze(0) # (1, 1, num_gaussians, n_features) + diff = x.unsqueeze(2) - centers # (batch_size, bag_size, num_gaussians, n_features) + + cov_inv = torch.inverse(self.covariance) + det_cov = torch.linalg.det(self.covariance) + + mahalanobis = torch.einsum('...i,...ij,...j->...', diff, cov_inv.unsqueeze(0).unsqueeze(0), diff) + normalization_term = torch.log((2 * torch.pi) ** self.n_features * det_cov).unsqueeze(0).unsqueeze(0) + log_probs = -0.5 * (mahalanobis + normalization_term) + return torch.exp(log_probs) # (batch_size, bag_size, num_gaussians) + + +class _GMBranch(nn.Module): + """One GM branch: an optional small MLP mapping the (already feature-extracted) instances into a + `gaussian_dimensions`-sized latent space, followed by a Sigmoid, a :class:`_GMLayer`, and a + BatchNorm applied instance-wise (i.e., over the merged batch*bag_size dimension, matching the + reference implementation). + """ + + def __init__(self, input_size, num_gaussians, gaussian_dimensions, hidden_size_fe, dropout_fe): + super().__init__() + self.pre = nn.Sequential() + prev_size = input_size + latent_size = gaussian_dimensions if gaussian_dimensions is not None else input_size + if gaussian_dimensions is not None: + for j, layer_size in enumerate(hidden_size_fe or ()): + self.pre.add_module(f'hidden_{j}', nn.Linear(prev_size, layer_size)) + self.pre.add_module(f'leakyrelu_{j}', nn.LeakyReLU()) + self.pre.add_module(f'dropout_{j}', nn.Dropout(dropout_fe)) + prev_size = layer_size + self.pre.add_module('latent_linear', nn.Linear(prev_size, gaussian_dimensions)) + self.pre.add_module('sigmoid', nn.Sigmoid()) + + self.gm_layer = _GMLayer(n_features=latent_size, num_gaussians=num_gaussians) + self.batch_norm = nn.BatchNorm1d(num_features=num_gaussians) + self.output_size = num_gaussians + self.latent_activation = None # populated on every forward(), read by GMNet's CKA regularization + + def forward(self, x): + # x: (batch_size, bag_size, input_size) + batch_size, bag_size = x.shape[0], x.shape[1] + latent = self.pre(x) + self.latent_activation = latent.reshape(-1, latent.shape[-1]) + likelihoods = self.gm_layer(latent) # (batch_size, bag_size, num_gaussians) + flat = self.batch_norm(likelihoods.reshape(batch_size * bag_size, -1)) + return flat.reshape(batch_size, bag_size, -1) + + +class _GMNetModule(nn.Module): + """The quantification module for GMNet: one or more :class:`_GMBranch` instances, each producing a + per-instance representation that is concatenated across branches and mean-pooled over the bag, as + required by :class:`quapy.method._neural_bags.BagTrainedQuantifier`. + """ + + def __init__(self, input_size, num_gaussians, n_gm_layers, gaussian_dimensions, hidden_size_fe=None, + dropout_fe=0., cka_regularization=0.): + super().__init__() + if len(num_gaussians) != n_gm_layers: + raise ValueError('num_gaussians should be a tuple of the same size as n_gm_layers') + if len(gaussian_dimensions) != n_gm_layers: + raise ValueError('gaussian_dimensions should be a tuple of the same size as n_gm_layers') + + self.n_gm_layers = n_gm_layers + self.cka_regularization = cka_regularization + self.branches = nn.ModuleList([ + _GMBranch(input_size, num_gaussians[i], gaussian_dimensions[i], hidden_size_fe, dropout_fe) + for i in range(n_gm_layers) + ]) + self.output_size = sum(num_gaussians) + + def forward(self, x): + outputs = [branch(x) for branch in self.branches] + return torch.mean(torch.cat(outputs, dim=-1), dim=1) + + def apply_regularization(self): + """Whether the CKA regularization term should be added to the training loss: requires at least + two GM branches (CKA is a pairwise measure) and a nonzero `cka_regularization` weight.""" + return self.n_gm_layers > 1 and self.cka_regularization != 0 + + def regularization_term(self): + latent_activations = [branch.latent_activation for branch in self.branches] + return self.cka_regularization * _cka(latent_activations) + + +class GMNet(BagTrainedQuantifier): + """ + Implementation of `GMNet `_, a neural network for quantification + that represents each instance of a bag by its likelihood under one or more learned mixtures of + Gaussians, mean-pools these representations over the bag, and predicts the class prevalence from the + result, trained end-to-end by minimizing a quantification loss over many samples ("bags") of known + prevalence. + + Like :class:`quapy.method._histnet.HistNetQ` and :class:`quapy.method.meta.QuaNet`, GMNet does not + follow the classify-then-aggregate pattern of :class:`quapy.method.aggregative.AggregativeQuantifier`; + it is instead trained and evaluated end-to-end on whole bags (see + :class:`quapy.method._neural_bags.BagTrainedQuantifier` for the shared training/prediction logic, + including the two entry points, :meth:`fit` and :meth:`fit_from_samples`). + + :param feature_extraction_module: a `torch.nn.Module` exposing an `output_size` attribute, used to + embed each instance before it is passed to every GM branch. If None (default), an identity + module is used, i.e., the instances in `X` are assumed to already be in their final numeric + representation. + :param n_gm_layers: number of GM branches (default 1). + :param num_gaussians: number of gaussians per branch: either a single int (used for every branch) or + a tuple/list of `n_gm_layers` ints (default 4). + :param gaussian_dimensions: dimensionality of the latent space in which each branch's gaussians live: + either a single int/None (used for every branch) or a tuple/list of `n_gm_layers` int/None + values. If None for a given branch, that branch's gaussians operate directly on the + feature-extracted instances, with no extra per-branch projection (default None). + :param hidden_size_fe: sizes of the hidden layers of the small per-branch MLP that maps the + feature-extracted instances into the latent space (only used when `gaussian_dimensions` is not + None for the corresponding branch); default None (no hidden layers, i.e., a single linear + projection). + :param dropout_fe: dropout applied after each of the `hidden_size_fe` layers (default 0). + :param cka_regularization: weight of the CKA regularization term encouraging the different branches' + latent representations to be dissimilar; only applied when `n_gm_layers > 1` (default 0, i.e., + disabled). + :param linear_sizes: tuple of ints with the sizes of the linear layers used in the shared + quantification head, after the GM branches (default empty, i.e., only the final classification + layer is used). + :param dropout: dropout applied after each of the `linear_sizes` layers (default 0). + :param output_function: either 'softmax' or 'normalize' (L1); both yield a valid prevalence vector + (default 'softmax'). + :param bag_size: number of instances per training/validation bag (default 500). + :param n_bags_train: number of bags generated per training epoch (default 500). + :param n_bags_val: number of bags generated per validation epoch (default 500). + :param train_epochs: maximum number of training epochs (default 200). + :param patience: number of epochs without improvement in validation loss before early-stopping + (default 20). + :param start_lr: initial learning rate (default 1e-3). + :param end_lr: once the learning rate decays below this value, training stops (default 1e-6). + :param lr_factor: factor by which the learning rate is reduced after `patience` epochs without + improvement (default 0.1). + :param weight_decay: L2 regularization (default 0). + :param quant_loss: the quantification loss to minimize (default `torch.nn.L1Loss()`), called as + `quant_loss(true_prevalences, predicted_prevalences)`. + :param batch_size: number of bags per gradient update (default 16). + :param protocol: the :class:`quapy.protocol.AbstractStochasticSeededProtocol` subclass used by + :meth:`fit` to resample bags from the given labelled collection (default + :class:`quapy.protocol.UPP`, which draws bags with prevalence sampled uniformly at random from + the simplex). + :param protocol_params: dict of extra keyword arguments passed to `protocol` (besides `data`, + `sample_size`, `repeats`, and `random_state`, which are set internally); default None. + :param val_split: float in (0,1), the proportion of the collection given to :meth:`fit` that is held + out (via stratified sampling) for validation and early stopping (default 0.4). + :param device: `'cpu'` or `'cuda'` (default 'cpu'). + :param random_state: seed used for the train/validation split and for the (fixed) validation + sampling sequence, as well as for the random initialization of the GM branches (default 0). + :param checkpointdir: directory where the best model found during training is stored (default + '../checkpoint'). + :param checkpointname: name of the checkpoint file; if None (default), a random name is generated. + :param verbose: verbosity level; if >0, shows a progress bar with the current losses (default 0). + """ + + def __init__(self, + feature_extraction_module=None, + n_gm_layers=1, + num_gaussians=4, + gaussian_dimensions=None, + hidden_size_fe=None, + dropout_fe=0., + cka_regularization=0., + linear_sizes=(), + dropout=0., + output_function='softmax', + bag_size=500, + n_bags_train=500, + n_bags_val=500, + train_epochs=200, + patience=20, + start_lr=1e-3, + end_lr=1e-6, + lr_factor=0.1, + weight_decay=0., + quant_loss=None, + batch_size=16, + protocol=UPP, + protocol_params=None, + val_split=0.4, + device='cpu', + random_state=0, + checkpointdir='../checkpoint', + checkpointname=None, + verbose=0): + super().__init__( + feature_extraction_module=feature_extraction_module, + linear_sizes=linear_sizes, + dropout=dropout, + output_function=output_function, + bag_size=bag_size, + n_bags_train=n_bags_train, + n_bags_val=n_bags_val, + train_epochs=train_epochs, + patience=patience, + start_lr=start_lr, + end_lr=end_lr, + lr_factor=lr_factor, + weight_decay=weight_decay, + quant_loss=quant_loss, + batch_size=batch_size, + protocol=protocol, + protocol_params=protocol_params, + val_split=val_split, + device=device, + random_state=random_state, + checkpointdir=checkpointdir, + checkpointname=checkpointname, + verbose=verbose, + ) + self.n_gm_layers = n_gm_layers + self.num_gaussians = num_gaussians if isinstance(num_gaussians, (tuple, list)) \ + else [num_gaussians] * n_gm_layers + self.gaussian_dimensions = gaussian_dimensions if isinstance(gaussian_dimensions, (tuple, list)) \ + else [gaussian_dimensions] * n_gm_layers + self.hidden_size_fe = hidden_size_fe + self.dropout_fe = dropout_fe + self.cka_regularization = cka_regularization + + @property + def _checkpoint_prefix(self): + return 'GMNet' + + def _build_quantmodule(self, n_features): + torch.manual_seed(self.random_state) + return _GMNetModule( + input_size=n_features, + num_gaussians=self.num_gaussians, + n_gm_layers=self.n_gm_layers, + gaussian_dimensions=self.gaussian_dimensions, + hidden_size_fe=self.hidden_size_fe, + dropout_fe=self.dropout_fe, + cka_regularization=self.cka_regularization, + ) + + def _extra_loss(self): + quantmodule = self.model.quantmodule + if quantmodule.apply_regularization(): + return quantmodule.regularization_term() + return 0. diff --git a/quapy/method/_histnet.py b/quapy/method/_histnet.py index ac4cc40..d718464 100644 --- a/quapy/method/_histnet.py +++ b/quapy/method/_histnet.py @@ -15,33 +15,15 @@ arXiv preprint arXiv:2012.06311 (2020). The overall architecture is: feature_extraction -> Sigmoid -> histogram layer -> small MLP -> softmax, trained by minimizing a quantification loss over samples ("bags") of known prevalence, rather than -over individually labeled instances (in the spirit of QuaNet, see method/_quanet.py). +over individually labeled instances (in the spirit of QuaNet, see method/_quanet.py). The bag-based +training loop itself (bag generation, early stopping, LR scheduling, checkpointing, prediction) is +shared with :class:`quapy.method._gmnet.GMNet` via :class:`quapy.method._neural_bags.BagTrainedQuantifier`. """ -import copy -import os -import random - -import numpy as np import torch import torch.nn as nn -from tqdm import tqdm -from quapy.data import LabelledCollection -from quapy.method.base import BaseQuantifier -from quapy.protocol import AbstractProtocol, UPP -from quapy.util import EarlyStop - - -class _IdentityFeatureExtractionModule(nn.Module): - """Used when no feature extraction module is provided: instances are assumed to already be in - their final numeric representation.""" - - def __init__(self, input_size): - super().__init__() - self.output_size = input_size - - def forward(self, x): - return x +from quapy.method._neural_bags import BagTrainedQuantifier +from quapy.protocol import UPP class _HardHistogramLayer(nn.Module): @@ -99,86 +81,24 @@ class _HardHistogramLayer(nn.Module): return result -class _HistNetModule(nn.Module): - """The full HistNetQ network: feature extraction, histogram, and the quantification MLP.""" +class _SigmoidHistogram(nn.Module): + """The quantification module for HistNetQ: squashes the (already feature-extracted) instances + through a Sigmoid and builds a differentiable histogram of them, as required by + :class:`quapy.method._neural_bags.BagTrainedQuantifier`.""" - def __init__(self, feature_extraction_module, n_classes, n_bins=8, quantiles=False, linear_sizes=(), - dropout=0., output_function='softmax'): + def __init__(self, n_features, n_bins=8, quantiles=False): super().__init__() - self.feature_extraction_module = feature_extraction_module self.sigmoid = nn.Sigmoid() - self.histogram = _HardHistogramLayer( - n_features=feature_extraction_module.output_size, n_bins=n_bins, quantiles=quantiles - ) + self.histogram = _HardHistogramLayer(n_features=n_features, n_bins=n_bins, quantiles=quantiles) + self.output_size = self.histogram.output_size - self.output_function = output_function - self.output_module = nn.Sequential() - prev_size = self.histogram.output_size - for i, linear_size in enumerate(linear_sizes): - self.output_module.add_module(f'linear_{i}', nn.Linear(prev_size, linear_size)) - self.output_module.add_module(f'leakyrelu_{i}', nn.LeakyReLU()) - self.output_module.add_module(f'dropout_{i}', nn.Dropout(dropout)) - prev_size = linear_size - self.output_module.add_module('last_linear', nn.Linear(prev_size, n_classes)) - if output_function == 'softmax': - self.output_module.add_module('softmax', nn.Softmax(dim=1)) - elif output_function == 'normalize': - self.output_module.add_module('relu', nn.ReLU()) - else: - raise ValueError(f"unknown {output_function=}; valid ones are 'softmax', 'normalize'") - - def forward(self, bag): - # bag: (batch_size, bag_size, n_features) - features = self.feature_extraction_module(bag) - features = self.sigmoid(features) - histogram = self.histogram(features) - out = self.output_module(histogram) - if self.output_function == 'normalize': - out = nn.functional.normalize(out, p=1, dim=1) - return out + def forward(self, input): + return self.histogram(self.sigmoid(input)) -def _to_tensor(x, device): - if torch.is_tensor(x): - return x.to(device=device, dtype=torch.float32) - if hasattr(x, 'toarray'): # scipy sparse - x = x.toarray() - return torch.as_tensor(np.asarray(x), dtype=torch.float32, device=device) - - -def _stack_bags(bags, device): +class HistNetQ(BagTrainedQuantifier): """ - :param bags: an iterable of (X_bag, prevalence) pairs, all X_bag with the same number of instances - :return: a pair of tensors (X, P) of shape (n_bags, bag_size, n_features) and (n_bags, n_classes) - """ - Xs, ps = zip(*bags) - X = torch.stack([_to_tensor(x, device) for x in Xs]) - P = torch.stack([_to_tensor(p, device) for p in ps]) - return X, P - - -def _mix_two_bags(bag_a, bag_b, bag_size, rng): - """Synthesizes a new bag of size `bag_size` by mixing two given bags with a random ratio, following - the "mixer" idea from the original HistNetQ repo (`UnlabeledMixerBagGenerator`): useful when the - only available training material is a modest number of pre-built samples (e.g., LeQua's dev - samples) and one wants extra intermediate-prevalence bags without access to instance-level labels. - """ - Xa, pa = bag_a - Xb, pb = bag_b - m = rng.random() - na = round(m * bag_size) - nb = bag_size - na - idx_a = rng.choices(range(len(Xa)), k=na) if na > 0 else [] - idx_b = rng.choices(range(len(Xb)), k=nb) if nb > 0 else [] - Xa, Xb = np.asarray(Xa), np.asarray(Xb) - X_mixed = np.concatenate([Xa[idx_a], Xb[idx_b]], axis=0) - p_mixed = m * np.asarray(pa, dtype=float) + (1 - m) * np.asarray(pb, dtype=float) - return X_mixed, p_mixed - - -class HistNetQ(BaseQuantifier): - """ - Implementation of `HistNetQ `_, a neural network for + Implementation of `HistNetQ `_, a neural network for quantification that learns a differentiable histogram-based representation of a sample, trained end-to-end by minimizing a quantification loss over many samples ("bags") of known prevalence. The method was proposed in `Pérez-Mon, O., Moreo, A., Coz, JJ del, & González, P. (2025). @@ -193,14 +113,9 @@ class HistNetQ(BaseQuantifier): end-to-end on whole samples rather than on individually labeled instances, following a symmetric problem setting (learning from bags, predicting on bags). - Training data can be provided in two ways: - - * via :meth:`fit`, from a plain labelled collection (`X`, `y`): training/validation bags are then - generated by resampling from it using a QuaPy sampling protocol (:class:`quapy.protocol.UPP` by - default). - * via :meth:`fit_from_samples`, from a :class:`quapy.protocol.AbstractProtocol` that already yields - the training bags (e.g., :class:`quapy.data._lequa.SamplesFromDir` for LeQua-style pre-built - samples), optionally enriched with synthetic bags mixed from the given ones. + Training data can be provided in two ways: via :meth:`fit`, from a plain labelled collection; or via + :meth:`fit_from_samples`, from a protocol that already yields the training bags. See + :class:`quapy.method._neural_bags.BagTrainedQuantifier` for details. :param feature_extraction_module: a `torch.nn.Module` exposing an `output_size` attribute, used to embed each instance before computing the histogram (e.g., a small MLP for tabular data, a CNN @@ -270,213 +185,37 @@ class HistNetQ(BaseQuantifier): checkpointdir='../checkpoint', checkpointname=None, verbose=0): - self.feature_extraction_module = feature_extraction_module + super().__init__( + feature_extraction_module=feature_extraction_module, + linear_sizes=linear_sizes, + dropout=dropout, + output_function=output_function, + bag_size=bag_size, + n_bags_train=n_bags_train, + n_bags_val=n_bags_val, + train_epochs=train_epochs, + patience=patience, + start_lr=start_lr, + end_lr=end_lr, + lr_factor=lr_factor, + weight_decay=weight_decay, + quant_loss=quant_loss, + batch_size=batch_size, + protocol=protocol, + protocol_params=protocol_params, + val_split=val_split, + device=device, + random_state=random_state, + checkpointdir=checkpointdir, + checkpointname=checkpointname, + verbose=verbose, + ) self.n_bins = n_bins self.quantiles = quantiles - self.linear_sizes = linear_sizes - self.dropout = dropout - self.output_function = output_function - self.bag_size = bag_size - self.n_bags_train = n_bags_train - self.n_bags_val = n_bags_val - self.train_epochs = train_epochs - self.patience = patience - self.start_lr = start_lr - self.end_lr = end_lr - self.lr_factor = lr_factor - self.weight_decay = weight_decay - self.quant_loss = quant_loss if quant_loss is not None else torch.nn.L1Loss() - self.batch_size = batch_size - self.protocol = protocol - self.protocol_params = protocol_params - self.val_split = val_split - self.device = torch.device(device) - self.random_state = random_state - if checkpointname is None: - local_random = random.Random() - random_code = '-'.join(str(local_random.randint(0, 1000000)) for _ in range(5)) - checkpointname = 'HistNetQ-' + random_code - self.checkpointdir = checkpointdir - self.checkpoint = os.path.join(checkpointdir, checkpointname) - self.verbose = verbose - self._classes_ = None @property - def classes_(self): - return self._classes_ + def _checkpoint_prefix(self): + return 'HistNetQ' - def fit(self, X, y): - """ - Trains HistNetQ from a plain labelled collection, generating training and validation bags by - resampling from it via `self.protocol` (a fresh random sequence of bags every epoch for - training, and a fixed, reproducible sequence for validation). - - :param X: the training instances - :param y: the labels of X - :return: self - """ - data = LabelledCollection(X, y) - self._classes_ = data.classes_ - train_data, val_data = data.split_stratified(train_prop=1 - self.val_split, random_state=self.random_state) - - protocol_params = self.protocol_params or {} - - def train_bags(): - sampler = self.protocol( - train_data, sample_size=self.bag_size, repeats=self.n_bags_train, random_state=None, - **protocol_params - ) - return sampler() - - def val_bags(): - sampler = self.protocol( - val_data, sample_size=self.bag_size, repeats=self.n_bags_val, random_state=self.random_state, - **protocol_params - ) - return sampler() - - n_features = train_data.instances.shape[1] - self._fit_loop(train_bags, val_bags, n_features, n_bags_train=self.n_bags_train, n_bags_val=self.n_bags_val) - return self - - def fit_from_samples(self, protocol: AbstractProtocol, val_protocol: AbstractProtocol = None, - mix_bags=False, mix_bags_proportion=0.5): - """ - Trains HistNetQ from a protocol that already yields the training bags (e.g., - :class:`quapy.data._lequa.SamplesFromDir`, for LeQua-style pre-built samples), instead of - resampling from a labelled collection. This is the entry point to use whenever only bags of - known prevalence are available (no instance-level labels). - - :param protocol: an :class:`AbstractProtocol` yielding `(sample, prevalence)` pairs; consumed - once and kept in memory (expected to be of modest size, as is typical of pre-built sample - collections). - :param val_protocol: an optional, separate protocol providing the validation bags; if None, a - `val_split` fraction of the bags returned by `protocol` is held out instead. - :param mix_bags: if True, in addition to the bags returned by `protocol`, synthesize extra bags - each epoch by mixing random pairs of the given bags with a random ratio (a substitute for - the original repo's `UnlabeledMixerBagGenerator`, useful to broaden the coverage of - prevalence values beyond what the given bags exhibit). - :param mix_bags_proportion: proportion (relative to the number of base training bags) of extra - mixed bags to generate per epoch when `mix_bags=True` (default 0.5). - :return: self - """ - assert isinstance(protocol, AbstractProtocol), 'protocol must be an instance of AbstractProtocol' - base_bags = list(protocol()) - n_classes = len(np.asarray(base_bags[0][1])) - self._classes_ = np.arange(n_classes) - - if val_protocol is not None: - val_bags_list = list(val_protocol()) - else: - n_val = max(1, int(len(base_bags) * self.val_split)) - val_bags_list = base_bags[:n_val] - base_bags = base_bags[n_val:] - - rng = random.Random(self.random_state) - n_mixed = round(len(base_bags) * mix_bags_proportion) if mix_bags else 0 - - def train_bags(): - bags = list(base_bags) - if n_mixed > 0: - for _ in range(n_mixed): - a, b = rng.choice(base_bags), rng.choice(base_bags) - bags.append(_mix_two_bags(a, b, self.bag_size, rng)) - rng.shuffle(bags) - return bags - - def val_bags(): - return val_bags_list - - n_features = np.asarray(base_bags[0][0]).shape[1] - self._fit_loop( - train_bags, val_bags, n_features, - n_bags_train=len(base_bags) + n_mixed, n_bags_val=len(val_bags_list) - ) - return self - - def _fit_loop(self, train_bags_fn, val_bags_fn, n_features, n_bags_train, n_bags_val): - os.makedirs(self.checkpointdir, exist_ok=True) - n_classes = len(self._classes_) - - fe = self.feature_extraction_module - if fe is None: - fe = _IdentityFeatureExtractionModule(n_features) - self.histnet = _HistNetModule( - fe, n_classes, n_bins=self.n_bins, quantiles=self.quantiles, linear_sizes=self.linear_sizes, - dropout=self.dropout, output_function=self.output_function - ).to(self.device) - - optim = torch.optim.Adam(self.histnet.parameters(), lr=self.start_lr, weight_decay=self.weight_decay) - scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optim, patience=self.patience, factor=self.lr_factor) - early_stop = EarlyStop(self.patience, lower_is_better=True) - - best_state = copy.deepcopy(self.histnet.state_dict()) - for epoch in range(self.train_epochs): - tr_loss = self._run_epoch(train_bags_fn(), n_bags_train, optim, train=True, epoch=epoch) - va_loss = self._run_epoch(val_bags_fn(), n_bags_val, optim=None, train=False, epoch=epoch) - - early_stop(va_loss, epoch) - if early_stop.IMPROVED: - best_state = copy.deepcopy(self.histnet.state_dict()) - torch.save(best_state, self.checkpoint) - elif early_stop.STOP: - if self.verbose > 0: - print(f'[HistNetQ] training ended by patience exhausted at epoch {epoch}; ' - f'restoring best model from epoch {early_stop.best_epoch}') - break - - scheduler.step(va_loss) - if optim.param_groups[0]['lr'] < self.end_lr: - if self.verbose > 0: - print(f'[HistNetQ] early stopping in epoch {epoch} (learning rate below end_lr)') - break - - self.histnet.load_state_dict(best_state) - - def _run_epoch(self, bags, n_bags, optim, train, epoch): - self.histnet.train(mode=train) - losses = [] - pbar = tqdm(bags, total=n_bags, disable=self.verbose == 0) - batch = [] - - def process_batch(batch): - X, P = _stack_bags(batch, self.device) - if train: - optim.zero_grad() - P_hat = self.histnet.forward(X) - loss = self.quant_loss(P, P_hat) - loss.backward() - optim.step() - else: - with torch.no_grad(): - P_hat = self.histnet.forward(X) - loss = self.quant_loss(P, P_hat) - return loss.item() - - for bag in pbar: - batch.append(bag) - if len(batch) == self.batch_size: - losses.append(process_batch(batch)) - batch = [] - pbar.set_description( - f'[HistNetQ] epoch={epoch} {"train" if train else "val"}-loss={np.mean(losses):.5f}' - ) - if batch: - losses.append(process_batch(batch)) - - return np.mean(losses) if losses else float('inf') - - def predict(self, X): - """ - Generates a class prevalence estimate for the sample `X`, via a single forward pass of the - trained network (the histogram layer aggregates over however many instances are given, so `X` - need not match the `bag_size` used during training). - - :param X: the test instances - :return: `np.ndarray` of shape `(n_classes,)` with the class prevalence estimates - """ - self.histnet.eval() - with torch.no_grad(): - X_t = _to_tensor(X, self.device).unsqueeze(0) - prevalence = self.histnet.forward(X_t) - return prevalence.cpu().numpy().flatten() + def _build_quantmodule(self, n_features): + return _SigmoidHistogram(n_features=n_features, n_bins=self.n_bins, quantiles=self.quantiles) diff --git a/quapy/method/_liep.py b/quapy/method/_liep.py new file mode 100644 index 0000000..3c28203 --- /dev/null +++ b/quapy/method/_liep.py @@ -0,0 +1,227 @@ +import warnings + +import numpy as np +from sklearn.base import BaseEstimator + +import quapy.functional as F +from quapy.method._helper import _labels_to_indices +from quapy.method.aggregative import AggregativeSoftQuantifier + + +def _normalize(v, eps=1e-12): + v = np.maximum(np.asarray(v, dtype=float), eps) + return v / v.sum() + + +def _bayes_prior_update(probs, new_prior, source_prior, eps=1e-12): + """ + Applies the standard label-shift prior correction: + + p_new(y|x) \\propto p_old(y|x) * new_prior(y) / source_prior(y) + + :param probs: array of shape `(n_instances, n_classes)` with the (uncorrected) posterior probabilities + :param new_prior: array of shape `(n_classes,)`, the current estimate of the target prior + :param source_prior: array of shape `(n_classes,)`, the source (training) prior + :param eps: numerical stabilizer + :return: array of shape `(n_instances, n_classes)` with the corrected posterior probabilities + """ + weights = _normalize(new_prior, eps) / np.maximum(_normalize(source_prior, eps), eps) + updated = probs * weights[np.newaxis, :] + updated /= np.maximum(updated.sum(axis=1, keepdims=True), eps) + return updated + + +def _confusion_statistic(y_true, y_pred, n_classes, mode='recall', eps=1e-12): + """ + Computes the minimum diagonal statistic of the confusion matrix, used to choose the fraction of target + instances to be retained in the high-confidence set (LEIP, Section 4.1). + + :param y_true: array of true label indices in `[0, n_classes)` + :param y_pred: array of predicted label indices in `[0, n_classes)` + :param n_classes: number of classes + :param mode: "recall" (per-true-class recall, i.e., diag(C) / row sums) or "precision" (per-predicted-class + precision, i.e., diag(C) / column sums) + :param eps: numerical stabilizer + :return: float, the minimum, across classes with at least one valid instance, of the chosen statistic + """ + C = np.zeros((n_classes, n_classes), dtype=float) + for yt, yp in zip(y_true, y_pred): + C[yt, yp] += 1. + + denom = C.sum(axis=1) if mode == 'recall' else C.sum(axis=0) + diag = np.diag(C) + valid = denom > 0 + if not np.any(valid): + raise ValueError('no valid classes found in the confusion matrix') + + return float(np.min(diag[valid] / np.maximum(denom[valid], eps))) + + +class LEIP(AggregativeSoftQuantifier): + """ + `Label Shift Estimation With Incremental Prior update + `_ (LEIP). + + Zhang, Y., Batista, G., & Kanhere, S.S. (2025). Label Shift Estimation With Incremental Prior Update. + In Proceedings of the 2025 SIAM International Conference on Data Mining (SDM), pp. 134-142. + + LEIP is a variant of :class:`EMQ` (aka SLD) that estimates the target prior without resorting to an + iterative Expectation-Maximization procedure. The target instances are first split into a high-confidence + set (those instances whose maximum posterior probability is at or above a threshold `tau`) and a + low-confidence set. The high-confidence set is used to obtain an initial classify-and-count estimate of + the target prior. The low-confidence instances are then relabelled one at a time, in order of decreasing + confidence, applying at each step the standard label-shift Bayes correction (see the module-level function + :func:`_bayes_prior_update`) using the running target-prior estimate; each newly assigned pseudo-label + updates this running estimate before the next instance is processed. Finally, a single Bayes correction is + applied to the whole target set using the resulting prior, and a hard-label classify-and-count pass on the + corrected posteriors yields the returned prevalence estimate. + + :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be + the one indicated in `qp.environ['DEFAULT_CLS']` + + :param fit_classifier: whether to train the classifier (default is True). Set to False if the + given classifier has already been trained. + + :param val_split: specifies the data used, when `tau=None`, for estimating the fraction of target + instances to retain in the high-confidence set (see `threshold_mode`). This specification can be made + as a float in (0, 1) indicating the proportion of stratified held-out validation set to be extracted + from the training set; or as an integer (default 5), indicating that the predictions are to be + generated in a `k`-fold cross-validation manner (with this integer indicating the value for `k`); or + as a tuple `(X, y)` defining the specific set of data to use for validation. This parameter is ignored + when a fixed `tau` is provided. + + :param tau: float or None (default). A fixed confidence threshold used to split the target instances into + the high- and low-confidence sets. If None, the threshold is instead derived, independently for each + target sample, as the quantile of that sample's own posterior-confidence distribution corresponding to + a retention fraction estimated from `val_split` (see `threshold_mode`). + + :param threshold_mode: either "recall" (default) or "precision", the statistic computed on the validation + confusion matrix that determines the fraction of target instances to be retained in the high-confidence + set (see the module-level function :func:`_confusion_statistic`). Only used when `tau=None`. + + :param count_smoothing: float, additive smoothing applied to the pseudo-label counts used to compute the + running target-prior estimate (default 0.0, matching the original paper; a small positive value, e.g., + 1e-8, can instead be used for extra numerical robustness). + """ + + THRESHOLD_MODES = ['recall', 'precision'] + + def __init__(self, classifier: BaseEstimator = None, fit_classifier=True, val_split=5, tau=None, + threshold_mode='recall', count_smoothing=0.0): + + assert threshold_mode in LEIP.THRESHOLD_MODES, \ + f'invalid {threshold_mode=}; valid ones are {LEIP.THRESHOLD_MODES}' + assert tau is None or isinstance(tau, (int, float)), f'invalid {tau=}; must be None or a number' + assert isinstance(count_smoothing, (int, float)) and count_smoothing >= 0, \ + f'invalid {count_smoothing=}; must be a non-negative number' + + super().__init__(classifier, fit_classifier, val_split) + self.tau = tau + self.threshold_mode = threshold_mode + self.count_smoothing = count_smoothing + + def _check_init_parameters(self): + if self.tau is None: + if self.val_split is None: + raise ValueError( + 'LEIP requires validation predictions to estimate the retention threshold when tau=None; ' + 'please set val_split to an integer, float, or validation tuple, or provide a fixed tau.' + ) + elif self.val_split is not None: + warnings.warn(f'a fixed {self.tau=} was set; {self.val_split=} will be ignored, since LEIP only ' + f'needs validation data to estimate a retention threshold when tau is not fixed.') + + def classify(self, X): + """ + Provides the posterior probabilities for the given instances. + + :param X: array-like of shape `(n_instances, n_dimensions,)` + :return: np.ndarray of shape `(n_instances, n_classes,)` with posterior probabilities + """ + return self.classifier.predict_proba(X) + + def classifier_fit_predict(self, X, y): + classif_predictions = super().classifier_fit_predict(X, y) + self.train_prevalence = F.prevalence_from_labels(y, classes=self.classes_) + return classif_predictions + + def aggregation_fit(self, classif_predictions, labels): + """ + Trains the aggregation function of LEIP. When a fixed `tau` was not provided at construction time, + this estimates, from the validation predictions, the fraction of target instances that should be + retained in the high-confidence set (see :func:`_confusion_statistic`). + + :param classif_predictions: array-like with the posterior probabilities of the validation instances + :param labels: array-like with the true labels associated to each classifier prediction + """ + if self.tau is None: + n_classes = len(self.classes_) + y_val_idx = _labels_to_indices(labels, self.classes_) + y_val_pred_idx = classif_predictions.argmax(axis=1) + self.retain_fraction_ = _confusion_statistic( + y_val_idx, y_val_pred_idx, n_classes, mode=self.threshold_mode + ) + else: + self.retain_fraction_ = None + + def aggregate(self, classif_posteriors): + prevalence, _ = self._leip(classif_posteriors) + return prevalence + + def predict_proba(self, instances): + """ + Returns the posterior probabilities updated by the final Bayes correction applied by LEIP. + + :param instances: np.ndarray of shape `(n_instances, n_dimensions)` + :return: np.ndarray of shape `(n_instances, n_classes)` + """ + classif_posteriors = self.classify(instances) + _, corrected_posteriors = self._leip(classif_posteriors) + return corrected_posteriors + + def _leip(self, target_probs): + n_classes = len(self.classes_) + source_prior = self.train_prevalence + + target_conf = target_probs.max(axis=1) + target_top = target_probs.argmax(axis=1) + + if self.tau is not None: + tau = self.tau + else: + retain_fraction = float(np.clip(self.retain_fraction_, 0., 1.)) + if retain_fraction <= 0: + tau = np.inf + elif retain_fraction >= 1: + tau = -np.inf + else: + tau = np.quantile(target_conf, 1. - retain_fraction) + self.tau_ = tau + + # high-confidence set A: an initial classify-and-count estimate of the target prior + A_mask = target_conf >= tau + counts = np.full(n_classes, self.count_smoothing, dtype=float) + if np.any(A_mask): + counts += np.bincount(target_top[A_mask], minlength=n_classes) + current_prior = counts / counts.sum() + else: + current_prior = source_prior.copy() + + # low-confidence set B, relabelled incrementally in order of decreasing confidence + B_indices = np.where(~A_mask)[0] + B_indices = B_indices[np.argsort(-target_conf[B_indices])] + + for idx in B_indices: + corrected = _bayes_prior_update(target_probs[idx:idx + 1], current_prior, source_prior)[0] + pseudo_label = int(np.argmax(corrected)) + counts[pseudo_label] += 1. + current_prior = counts / counts.sum() + + self.intermediate_prior_ = current_prior.copy() + + # final full-batch Bayes correction and hard-label classify-and-count pass + corrected_posteriors = _bayes_prior_update(target_probs, current_prior, source_prior) + final_labels = corrected_posteriors.argmax(axis=1) + prevalence = F.prevalence_from_labels(final_labels, classes=np.arange(n_classes)) + + return prevalence, corrected_posteriors diff --git a/quapy/method/_neural_bags.py b/quapy/method/_neural_bags.py new file mode 100644 index 0000000..f32af71 --- /dev/null +++ b/quapy/method/_neural_bags.py @@ -0,0 +1,443 @@ +""" +Shared machinery for QuaPy's "bag-trained" neural quantifiers, i.e., methods that -- like +:class:`quapy.method.meta.QuaNet` -- do not follow the classify-then-aggregate pattern of +:class:`quapy.method.aggregative.AggregativeQuantifier`, but are instead trained and evaluated +end-to-end on whole samples ("bags") of known prevalence. + +:class:`quapy.method._histnet.HistNetQ` and :class:`quapy.method._gmnet.GMNet` share the same overall +architecture (feature_extraction -> quantification module -> small MLP head -> softmax/normalize) and +the same bag-based training protocol (bag generation via a QuaPy sampling protocol, early stopping, LR +scheduling, checkpointing). This module factors that common part out into :class:`BagTrainedQuantifier`; +concrete subclasses only need to supply the quantification module placed between the feature extractor +and the shared head (see :meth:`BagTrainedQuantifier._build_quantmodule`). +""" +import copy +import os +import random +from abc import ABC, abstractmethod + +import numpy as np +import torch +import torch.nn as nn +from tqdm import tqdm + +from quapy.data import LabelledCollection +from quapy.method.base import BaseQuantifier +from quapy.protocol import AbstractProtocol, UPP +from quapy.util import EarlyStop + + +class IdentityFeatureExtractionModule(nn.Module): + """Used when no feature extraction module is provided: instances are assumed to already be in + their final numeric representation.""" + + def __init__(self, input_size): + super().__init__() + self.output_size = input_size + + def forward(self, x): + return x + + +def to_tensor(x, device): + if torch.is_tensor(x): + return x.to(device=device, dtype=torch.float32) + if hasattr(x, 'toarray'): # scipy sparse + x = x.toarray() + return torch.as_tensor(np.asarray(x), dtype=torch.float32, device=device) + + +def stack_bags(bags, device): + """ + :param bags: an iterable of (X_bag, prevalence) pairs, all X_bag with the same number of instances + :return: a pair of tensors (X, P) of shape (n_bags, bag_size, n_features) and (n_bags, n_classes) + """ + Xs, ps = zip(*bags) + X = torch.stack([to_tensor(x, device) for x in Xs]) + P = torch.stack([to_tensor(p, device) for p in ps]) + return X, P + + +def mix_two_bags(bag_a, bag_b, bag_size, rng): + """Synthesizes a new bag of size `bag_size` by mixing two given bags with a random ratio, following + the "mixer" idea from the original HistNetQ repo (`UnlabeledMixerBagGenerator`): useful when the + only available training material is a modest number of pre-built samples (e.g., LeQua's dev + samples) and one wants extra intermediate-prevalence bags without access to instance-level labels. + """ + Xa, pa = bag_a + Xb, pb = bag_b + m = rng.random() + na = round(m * bag_size) + nb = bag_size - na + idx_a = rng.choices(range(len(Xa)), k=na) if na > 0 else [] + idx_b = rng.choices(range(len(Xb)), k=nb) if nb > 0 else [] + Xa, Xb = np.asarray(Xa), np.asarray(Xb) + X_mixed = np.concatenate([Xa[idx_a], Xb[idx_b]], axis=0) + p_mixed = m * np.asarray(pa, dtype=float) + (1 - m) * np.asarray(pb, dtype=float) + return X_mixed, p_mixed + + +def build_output_head(input_size, n_classes, linear_sizes, dropout, output_function): + """Builds the small MLP + output activation shared by every bag-trained quantifier's head: a stack + of (Linear, LeakyReLU, Dropout) blocks sized by `linear_sizes`, followed by a final Linear to + `n_classes` and either a softmax or an L1-normalization (applied in :class:`BagNetworkModule`), both + yielding a valid prevalence vector.""" + output_module = nn.Sequential() + prev_size = input_size + for i, linear_size in enumerate(linear_sizes): + output_module.add_module(f'linear_{i}', nn.Linear(prev_size, linear_size)) + output_module.add_module(f'leakyrelu_{i}', nn.LeakyReLU()) + output_module.add_module(f'dropout_{i}', nn.Dropout(dropout)) + prev_size = linear_size + output_module.add_module('last_linear', nn.Linear(prev_size, n_classes)) + if output_function == 'softmax': + output_module.add_module('softmax', nn.Softmax(dim=1)) + elif output_function == 'normalize': + output_module.add_module('relu', nn.ReLU()) + else: + raise ValueError(f"unknown {output_function=}; valid ones are 'softmax', 'normalize'") + return output_module + + +class BagNetworkModule(nn.Module): + """The full network shared by every bag-trained quantifier: feature extraction, a pluggable + quantification module (mapping a bag of instance-level features to a single per-bag + representation), and the shared MLP head. + + :param quantmodule: a `torch.nn.Module` exposing an `output_size` attribute, mapping a tensor of + shape (batch_size, bag_size, n_features) to one of shape (batch_size, quantmodule.output_size). + """ + + def __init__(self, feature_extraction_module, quantmodule, n_classes, linear_sizes, dropout, output_function): + super().__init__() + self.feature_extraction_module = feature_extraction_module + self.quantmodule = quantmodule + self.output_function = output_function + self.output_module = build_output_head( + quantmodule.output_size, n_classes, linear_sizes, dropout, output_function + ) + + def forward(self, bag): + # bag: (batch_size, bag_size, n_features) + features = self.feature_extraction_module(bag) + representation = self.quantmodule(features) + out = self.output_module(representation) + if self.output_function == 'normalize': + out = nn.functional.normalize(out, p=1, dim=1) + return out + + +class BagTrainedQuantifier(BaseQuantifier, ABC): + """ + Base class for QuaPy's neural quantifiers trained end-to-end on whole samples ("bags") of known + prevalence, rather than following the classify-then-aggregate pattern of + :class:`quapy.method.aggregative.AggregativeQuantifier` (in the spirit of + :class:`quapy.method.meta.QuaNet`). Concrete subclasses only need to provide the quantification + module placed between the feature extractor and the shared MLP head (see + :meth:`_build_quantmodule`) and a checkpoint-name prefix (see :attr:`_checkpoint_prefix`); bag + generation, the training/validation loop, early stopping, LR scheduling, checkpointing, and + prediction are all shared. + + Training data can be provided in two ways: + + * via :meth:`fit`, from a plain labelled collection (`X`, `y`): training/validation bags are then + generated by resampling from it using a QuaPy sampling protocol (:class:`quapy.protocol.UPP` by + default). + * via :meth:`fit_from_samples`, from a :class:`quapy.protocol.AbstractProtocol` that already yields + the training bags (e.g., :class:`quapy.data._lequa.SamplesFromDir` for LeQua-style pre-built + samples), optionally enriched with synthetic bags mixed from the given ones. + + :param feature_extraction_module: a `torch.nn.Module` exposing an `output_size` attribute, used to + embed each instance before the quantification module (e.g., a small MLP for tabular data, a CNN + for images). If None (default), an identity module is used, i.e., the instances in `X` are + assumed to already be in their final numeric representation. + :param linear_sizes: tuple of ints with the sizes of the linear layers used in the shared head + (default empty, i.e., only the final classification layer is used). + :param dropout: dropout applied after each of the `linear_sizes` layers (default 0). + :param output_function: either 'softmax' or 'normalize' (L1); both yield a valid prevalence vector + (default 'softmax'). + :param bag_size: number of instances per training/validation bag (default 500). + :param n_bags_train: number of bags generated per training epoch (default 500). + :param n_bags_val: number of bags generated per validation epoch (default 500). + :param train_epochs: maximum number of training epochs (default 200). + :param patience: number of epochs without improvement in validation loss before early-stopping + (default 20). + :param start_lr: initial learning rate (default 1e-3). + :param end_lr: once the learning rate decays below this value, training stops (default 1e-6). + :param lr_factor: factor by which the learning rate is reduced after `patience` epochs without + improvement (default 0.1). + :param weight_decay: L2 regularization (default 0). + :param quant_loss: the quantification loss to minimize (default `torch.nn.L1Loss()`), called as + `quant_loss(true_prevalences, predicted_prevalences)`. + :param batch_size: number of bags per gradient update (default 16). + :param protocol: the :class:`quapy.protocol.AbstractStochasticSeededProtocol` subclass used by + :meth:`fit` to resample bags from the given labelled collection (default + :class:`quapy.protocol.UPP`, which draws bags with prevalence sampled uniformly at random from + the simplex). + :param protocol_params: dict of extra keyword arguments passed to `protocol` (besides `data`, + `sample_size`, `repeats`, and `random_state`, which are set internally); default None. + :param val_split: float in (0,1), the proportion of the collection given to :meth:`fit` that is held + out (via stratified sampling) for validation and early stopping (default 0.4). + :param device: `'cpu'` or `'cuda'` (default 'cpu'). + :param random_state: seed used for the train/validation split and for the (fixed) validation + sampling sequence (default 0). + :param checkpointdir: directory where the best model found during training is stored (default + '../checkpoint'). + :param checkpointname: name of the checkpoint file; if None (default), a random name prefixed by + :attr:`_checkpoint_prefix` is generated. + :param verbose: verbosity level; if >0, shows a progress bar with the current losses (default 0). + """ + + def __init__(self, + feature_extraction_module=None, + linear_sizes=(), + dropout=0., + output_function='softmax', + bag_size=500, + n_bags_train=500, + n_bags_val=500, + train_epochs=200, + patience=20, + start_lr=1e-3, + end_lr=1e-6, + lr_factor=0.1, + weight_decay=0., + quant_loss=None, + batch_size=16, + protocol=UPP, + protocol_params=None, + val_split=0.4, + device='cpu', + random_state=0, + checkpointdir='../checkpoint', + checkpointname=None, + verbose=0): + self.feature_extraction_module = feature_extraction_module + self.linear_sizes = linear_sizes + self.dropout = dropout + self.output_function = output_function + self.bag_size = bag_size + self.n_bags_train = n_bags_train + self.n_bags_val = n_bags_val + self.train_epochs = train_epochs + self.patience = patience + self.start_lr = start_lr + self.end_lr = end_lr + self.lr_factor = lr_factor + self.weight_decay = weight_decay + self.quant_loss = quant_loss if quant_loss is not None else torch.nn.L1Loss() + self.batch_size = batch_size + self.protocol = protocol + self.protocol_params = protocol_params + self.val_split = val_split + self.device = torch.device(device) + self.random_state = random_state + if checkpointname is None: + local_random = random.Random() + random_code = '-'.join(str(local_random.randint(0, 1000000)) for _ in range(5)) + checkpointname = f'{self._checkpoint_prefix}-{random_code}' + self.checkpointdir = checkpointdir + self.checkpoint = os.path.join(checkpointdir, checkpointname) + self.verbose = verbose + self._classes_ = None + self.model = None + + @property + def classes_(self): + return self._classes_ + + @property + @abstractmethod + def _checkpoint_prefix(self): + """Short name used as the default checkpoint filename prefix (e.g. 'HistNetQ', 'GMNet').""" + ... + + @abstractmethod + def _build_quantmodule(self, n_features): + """Builds the module placed between the (already feature-extracted) instances and the shared + MLP head. Must expose an `output_size` attribute and accept input of shape + (batch_size, bag_size, n_features), returning one of shape (batch_size, output_size).""" + ... + + def _extra_loss(self): + """Optional additional term added to the quantification loss during training (e.g., GMNet's CKA + regularization across GM layers). Returns 0 by default.""" + return 0. + + def fit(self, X, y): + """ + Trains the quantifier from a plain labelled collection, generating training and validation bags + by resampling from it via `self.protocol` (a fresh random sequence of bags every epoch for + training, and a fixed, reproducible sequence for validation). + + :param X: the training instances + :param y: the labels of X + :return: self + """ + data = LabelledCollection(X, y) + self._classes_ = data.classes_ + train_data, val_data = data.split_stratified(train_prop=1 - self.val_split, random_state=self.random_state) + + protocol_params = self.protocol_params or {} + + def train_bags(): + sampler = self.protocol( + train_data, sample_size=self.bag_size, repeats=self.n_bags_train, random_state=None, + **protocol_params + ) + return sampler() + + def val_bags(): + sampler = self.protocol( + val_data, sample_size=self.bag_size, repeats=self.n_bags_val, random_state=self.random_state, + **protocol_params + ) + return sampler() + + n_features = train_data.instances.shape[1] + self._fit_loop(train_bags, val_bags, n_features, n_bags_train=self.n_bags_train, n_bags_val=self.n_bags_val) + return self + + def fit_from_samples(self, protocol: AbstractProtocol, val_protocol: AbstractProtocol = None, + mix_bags=False, mix_bags_proportion=0.5): + """ + Trains the quantifier from a protocol that already yields the training bags (e.g., + :class:`quapy.data._lequa.SamplesFromDir`, for LeQua-style pre-built samples), instead of + resampling from a labelled collection. This is the entry point to use whenever only bags of + known prevalence are available (no instance-level labels). + + :param protocol: an :class:`AbstractProtocol` yielding `(sample, prevalence)` pairs; consumed + once and kept in memory (expected to be of modest size, as is typical of pre-built sample + collections). + :param val_protocol: an optional, separate protocol providing the validation bags; if None, a + `val_split` fraction of the bags returned by `protocol` is held out instead. + :param mix_bags: if True, in addition to the bags returned by `protocol`, synthesize extra bags + each epoch by mixing random pairs of the given bags with a random ratio (a substitute for + the original HistNetQ repo's `UnlabeledMixerBagGenerator`, useful to broaden the coverage of + prevalence values beyond what the given bags exhibit). + :param mix_bags_proportion: proportion (relative to the number of base training bags) of extra + mixed bags to generate per epoch when `mix_bags=True` (default 0.5). + :return: self + """ + assert isinstance(protocol, AbstractProtocol), 'protocol must be an instance of AbstractProtocol' + base_bags = list(protocol()) + n_classes = len(np.asarray(base_bags[0][1])) + self._classes_ = np.arange(n_classes) + + if val_protocol is not None: + val_bags_list = list(val_protocol()) + else: + n_val = max(1, int(len(base_bags) * self.val_split)) + val_bags_list = base_bags[:n_val] + base_bags = base_bags[n_val:] + + rng = random.Random(self.random_state) + n_mixed = round(len(base_bags) * mix_bags_proportion) if mix_bags else 0 + + def train_bags(): + bags = list(base_bags) + if n_mixed > 0: + for _ in range(n_mixed): + a, b = rng.choice(base_bags), rng.choice(base_bags) + bags.append(mix_two_bags(a, b, self.bag_size, rng)) + rng.shuffle(bags) + return bags + + def val_bags(): + return val_bags_list + + n_features = np.asarray(base_bags[0][0]).shape[1] + self._fit_loop( + train_bags, val_bags, n_features, + n_bags_train=len(base_bags) + n_mixed, n_bags_val=len(val_bags_list) + ) + return self + + def _fit_loop(self, train_bags_fn, val_bags_fn, n_features, n_bags_train, n_bags_val): + os.makedirs(self.checkpointdir, exist_ok=True) + n_classes = len(self._classes_) + + fe = self.feature_extraction_module + if fe is None: + fe = IdentityFeatureExtractionModule(n_features) + quantmodule = self._build_quantmodule(fe.output_size) + self.model = BagNetworkModule( + fe, quantmodule, n_classes, linear_sizes=self.linear_sizes, dropout=self.dropout, + output_function=self.output_function + ).to(self.device) + + optim = torch.optim.Adam(self.model.parameters(), lr=self.start_lr, weight_decay=self.weight_decay) + scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optim, patience=self.patience, factor=self.lr_factor) + early_stop = EarlyStop(self.patience, lower_is_better=True) + + best_state = copy.deepcopy(self.model.state_dict()) + for epoch in range(self.train_epochs): + self._run_epoch(train_bags_fn(), n_bags_train, optim, train=True, epoch=epoch) + va_loss = self._run_epoch(val_bags_fn(), n_bags_val, optim=None, train=False, epoch=epoch) + + early_stop(va_loss, epoch) + if early_stop.IMPROVED: + best_state = copy.deepcopy(self.model.state_dict()) + torch.save(best_state, self.checkpoint) + elif early_stop.STOP: + if self.verbose > 0: + print(f'[{self._checkpoint_prefix}] training ended by patience exhausted at epoch {epoch}; ' + f'restoring best model from epoch {early_stop.best_epoch}') + break + + scheduler.step(va_loss) + if optim.param_groups[0]['lr'] < self.end_lr: + if self.verbose > 0: + print(f'[{self._checkpoint_prefix}] early stopping in epoch {epoch} (learning rate below end_lr)') + break + + self.model.load_state_dict(best_state) + + def _run_epoch(self, bags, n_bags, optim, train, epoch): + self.model.train(mode=train) + losses = [] + pbar = tqdm(bags, total=n_bags, disable=self.verbose == 0) + batch = [] + + def process_batch(batch): + X, P = stack_bags(batch, self.device) + if train: + optim.zero_grad() + P_hat = self.model.forward(X) + loss = self.quant_loss(P, P_hat) + self._extra_loss() + loss.backward() + optim.step() + else: + with torch.no_grad(): + P_hat = self.model.forward(X) + loss = self.quant_loss(P, P_hat) + return loss.item() + + for bag in pbar: + batch.append(bag) + if len(batch) == self.batch_size: + losses.append(process_batch(batch)) + batch = [] + pbar.set_description( + f'[{self._checkpoint_prefix}] epoch={epoch} {"train" if train else "val"}-' + f'loss={np.mean(losses):.5f}' + ) + if batch: + losses.append(process_batch(batch)) + + return np.mean(losses) if losses else float('inf') + + def predict(self, X): + """ + Generates a class prevalence estimate for the sample `X`, via a single forward pass of the + trained network (the quantification module aggregates over however many instances are given, so + `X` need not match the `bag_size` used during training). + + :param X: the test instances + :return: `np.ndarray` of shape `(n_classes,)` with the class prevalence estimates + """ + self.model.eval() + with torch.no_grad(): + X_t = to_tensor(X, self.device).unsqueeze(0) + prevalence = self.model.forward(X_t) + return prevalence.cpu().numpy().flatten() diff --git a/quapy/method/aggregative.py b/quapy/method/aggregative.py index f678546..6dcd025 100644 --- a/quapy/method/aggregative.py +++ b/quapy/method/aggregative.py @@ -368,6 +368,65 @@ class BinaryAggregativeQuantifier(AggregativeQuantifier, BinaryQuantifier): return super().fit(X, y) +class ImportanceWeightQuantifier(AggregativeQuantifier, ABC): + """ + Abstract mixin for aggregative quantifiers that estimate the target prevalence by first computing a vector + of importance weights :math:`w_y=Q(y)/P(y)` (with :math:`P` and :math:`Q` the training and target + distributions), and then rescaling the training prevalence by these weights, i.e., + :math:`\\hat{p}(y) \\propto w_y \\cdot P(y)`. + + Subclasses must set a fitted attribute `self.train_prevalence_` (typically in :meth:`aggregation_fit`) and + must implement :meth:`_weights_from_predictions`. This class provides a template implementation of + :meth:`aggregate`, together with :meth:`get_importance_weights` and :meth:`quantify_and_weigh`, none of + which mutate any internal state, so that they are all safe to call concurrently (e.g., from different + threads) on the same fitted instance for different batches of target instances. + """ + + @abstractmethod + def _weights_from_predictions(self, classif_predictions) -> np.ndarray: + """ + Computes the vector of importance weights from the classifier predictions on a batch of (unlabelled) + target instances. + + :param classif_predictions: array-like with the classifier predictions (crisp or soft, depending on + the subclass) for the target instances + :return: np.ndarray of shape `(n_classes,)` + """ + ... + + def _prevalence_from_weights(self, weights: np.ndarray) -> np.ndarray: + estimate = self.train_prevalence_ * weights + return F.normalize_prevalence(estimate, method=self.norm) + + def aggregate(self, classif_predictions): + weights = self._weights_from_predictions(classif_predictions) + return self._prevalence_from_weights(weights) + + def get_importance_weights(self, instances) -> np.ndarray: + """ + Estimates the vector of importance weights :math:`w_y=Q(y)/P(y)` for the given (unlabelled) target + instances. + + :param instances: array-like of shape `(n_instances, n_dimensions)`, the target instances + :return: np.ndarray of shape `(n_classes,)` + """ + classif_predictions = self.classify(instances) + return self._weights_from_predictions(classif_predictions) + + def quantify_and_weigh(self, instances): + """ + Jointly returns the estimated target prevalence and the importance weights used to obtain it, computed + from a single pass of classifier predictions over the given (unlabelled) target instances. + + :param instances: array-like of shape `(n_instances, n_dimensions)`, the target instances + :return: a tuple `(prevalence, weights)`, both np.ndarray of shape `(n_classes,)` + """ + classif_predictions = self.classify(instances) + weights = self._weights_from_predictions(classif_predictions) + prevalence = self._prevalence_from_weights(weights) + return prevalence, weights + + # Methods # ------------------------------------ class CC(AggregativeCrispQuantifier): @@ -456,8 +515,10 @@ class ACC(AggregativeCrispQuantifier): :param str method: adjustment method to be used: - * 'inversion': matrix inversion method based on the matrix equality :math:`P(C)=P(C|Y)P(Y)`, - which tries to invert :math:`P(C|Y)` matrix. + * 'inversion': matrix inversion method. Based on the matrix equality :math:`q=M p`, with + :math:`q` the prevalence vector estimated by CC, :math:`M` the matrix with entries :math:`i,j` representing + :math:`P(\\hat{Y}=i|Y=j)`, and :math:`p` the sought class prevalence vector, the matrix inversion + tries to solve for :math:`p=M^{-1} q` * 'invariant-ratio': invariant ratio estimator of `Vaz et al. 2018 `_, which replaces the last equation with the normalization condition. @@ -467,8 +528,7 @@ class ACC(AggregativeCrispQuantifier): strictly less than `n_classes`. * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) - * 'exact': deprecated, defaults to 'exact-cc' - * 'minimize': minimizes the L2 norm of :math:`|Ax-B|`. This one generally works better, and is the + * 'minimize': minimizes the squared L2 norm of :math:`|Ax-B|`. This one generally works better, and is the default parameter. More details about this can be consulted in `Bunse, M. "On Multi-Class Extensions of Adjusted Classify and Count", on proceedings of the 2nd International Workshop on Learning to Quantify: Methods and Applications (LQ 2022), ECML/PKDD 2022, Grenoble (France) @@ -489,7 +549,7 @@ class ACC(AggregativeCrispQuantifier): classifier: BaseEstimator = None, fit_classifier = True, val_split = 5, - solver: Literal['minimize', 'exact', 'exact-raise', 'exact-cc'] = 'minimize', + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', method: Literal['inversion', 'invariant-ratio'] = 'inversion', norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', n_jobs=None, @@ -500,7 +560,7 @@ class ACC(AggregativeCrispQuantifier): self.method = method self.norm = norm - SOLVERS = ['exact', 'minimize', 'exact-raise', 'exact-cc'] + SOLVERS = ['minimize', 'exact-raise', 'exact-cc'] METHODS = ['inversion', 'invariant-ratio'] NORMALIZATIONS = ['clip', 'mapsimplex', 'condsoftmax', None] @@ -551,8 +611,8 @@ class ACC(AggregativeCrispQuantifier): @classmethod def getPteCondEstim(cls, classes, y, y_): """ - Estimate the matrix with entry (i,j) being the estimate of P(hat_yi|yj), that is, the probability that a - document that belongs to yj ends up being classified as belonging to yi + Estimate the matrix with entry (i,j) being the estimate of P(hat_yi|yj), that is, the probability that an + instance that belongs to class j ends up being classified as belonging to class i :param classes: array-like with the class names :param y: array-like with the true labels @@ -570,15 +630,234 @@ class ACC(AggregativeCrispQuantifier): return conf def aggregate(self, classif_predictions): - prevs_estim = self.cc.aggregate(classif_predictions) - estimate = F.solve_adjustment( + prevs_estim_cc = self.cc.aggregate(classif_predictions) + estimate = ACC.solve_adjustment( class_conditional_rates=self.Pte_cond_estim_, - unadjusted_counts=prevs_estim, + unadjusted_counts=prevs_estim_cc, solver=self.solver, method=self.method, ) return F.normalize_prevalence(estimate, method=self.norm) + @classmethod + def solve_adjustment(cls, + class_conditional_rates: np.ndarray, + unadjusted_counts: np.ndarray, + method: Literal["inversion", "invariant-ratio"], + solver: Literal["minimize", "exact-raise", "exact-cc"]) -> np.ndarray: + """ + Function that tries to solve for :math:`p` the equation :math:`q = M p`, where :math:`q` is the vector of + `unadjusted counts` (as estimated, e.g., via classify and count) with :math:`q_i` an estimate of + :math:`P(\\hat{Y}=y_i)`, and where :math:`M` is the matrix of `class-conditional rates` with :math:`M_{ij}` an + estimate of :math:`P(\\hat{Y}=y_i|Y=y_j)`. + + :param class_conditional_rates: array of shape `(n_classes, n_classes,)` with entry `(i,j)` being the estimate + of :math:`P(\\hat{Y}=y_i|Y=y_j)`, that is, the probability that an instance ends up being classified as + belonging to class :math:`y_i` given it actually belonged to class :math:`y_j` + + :param unadjusted_counts: array of shape `(n_classes,)` containing the unadjusted prevalence values (e.g., as + estimated by CC or PCC) + + :param str method: indicates the adjustment method to be used. Valid options are: + + * `inversion`: tries to solve the equation :math:`q = M p` as :math:`p = M^{-1} q` where + :math:`M^{-1}` is the matrix inversion of :math:`M`. This inversion may not exist in + degenerated cases. + * `invariant-ratio`: invariant ratio estimator of `Vaz et al. 2018 `_, + which replaces the last equation in :math:`M` with the normalization condition (i.e., that the sum of + all prevalence values must equal 1). + + :param str solver: the method to use for solving the system of linear equations. Valid options are: + + * `exact-raise`: tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly lower than `n_classes`. + * `exact-cc`: if the matrix is not full rank, returns :math:`q` (i.e., the unadjusted counts) as the estimates + * `minimize`: minimizes a loss, so the solution always exists + """ + + A = np.asarray(class_conditional_rates, dtype=float).copy() + B = np.asarray(unadjusted_counts, dtype=float).copy() + + if method == "inversion": + pass # leave A and B unchanged + elif method == "invariant-ratio": + # Change the last equation to replace it with the normalization condition; + # copy first so this does not mutate the caller's arrays (np.asarray above + # returns the same object, not a copy, when the input is already float64) + A[-1, :] = 1.0 + B[-1] = 1.0 + else: + raise ValueError(f"unknown {method=}") + + if solver in ["exact-raise", "exact-cc"]: + try: + return np.linalg.solve(A, B) + except np.linalg.LinAlgError: + if solver=='exact-cc': + return unadjusted_counts + else: + raise + elif solver == "minimize": + def loss(prev): + return np.linalg.norm(A @ prev - B) + return F.optim_minimize(loss, n_classes=A.shape[0], return_loss=False) + else: + raise ValueError(f"unknown {solver=}") + + +class BBSEhard(ImportanceWeightQuantifier, AggregativeCrispQuantifier): + """ + `Black Box Shift Estimator` (BBSE) hard aims at finding the importance weights :math:`w_i=Q(i)/P(i)`, + with :math:`P` and :math:`Q` the training and test distributions. + BBSE is similar in spirit to ACC, but it solves the problem :math:`q=C w`, with :math:`q` the prevalence vector + estimated by CC, :math:`C` the matrix with entries :math:`i,j` representing :math:`P(\\hat{Y}=i,Y=j)`, and + :math:`w` the sought vector of importance weights. The `hard` variant estimates these quantities using + crisp counts. + + BBSE was proposed in + `Lipton, Z., Wang, Y. X., & Smola, A. (2018, July). + Detecting and correcting for label shift with black box predictors. + In International conference on machine learning (pp. 3122-3130). PMLR. + `_. + + :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be + the one indicated in `qp.environ['DEFAULT_CLS']` + + :param fit_classifier: whether to train the learner (default is True). Set to False if the + learner has been trained outside the quantifier. + + :param val_split: specifies the data used for generating classifier predictions. This specification + can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to + be extracted from the training set; or as an integer (default 5), indicating that the predictions + are to be generated in a `k`-fold cross-validation manner (with this integer indicating the value + for `k`); or as a tuple (X,y) defining the specific set of data to use for validation. + + :param str solver: indicates the method to use for solving the system of linear equations. Valid options are: + + * 'exact-raise': tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly less than `n_classes`. + * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to + no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) + * 'exact': deprecated, defaults to 'exact-cc' + * 'minimize': minimizes the squared L2 norm of :math:`|Ax-B|`. This one generally works better, and is the + default parameter. More details about this can be consulted in + `Tachet des Combes, R., Zhao, H., Wang, Y. X., & Gordon, G. J. (2020). + Domain adaptation with conditional distribution matching and generalized label shift. + Advances in Neural Information Processing Systems, 33, 19276-19289. + `_. + + :param str norm: the method to use for normalization. + + * `clip`, the values are clipped to the range [0,1] and then L1-normalized. + * `mapsimplex` projects vectors onto the probability simplex. This implementation relies on + `Mathieu Blondel's projection_simplex_sort `_ + * `condsoftmax`, applies a softmax normalization only to prevalence vectors that lie outside the simplex + + :param n_jobs: number of parallel workers + """ + + def __init__( + self, + classifier: BaseEstimator = None, + fit_classifier = True, + val_split = 5, + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', + norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', + n_jobs=None, + ): + super().__init__(classifier, fit_classifier, val_split) + self.n_jobs = qp._get_njobs(n_jobs) + self.solver = solver + self.norm = norm + + def _check_init_parameters(self): + if self.solver not in ACC.SOLVERS: + raise ValueError(f"unknown solver; valid ones are {ACC.SOLVERS}") + if self.norm not in ACC.NORMALIZATIONS: + raise ValueError(f"unknown normalization; valid ones are {ACC.NORMALIZATIONS}") + + def aggregation_fit(self, classif_predictions, labels): + """ + Estimates the misclassification rates. + :param classif_predictions: array-like with the predicted labels + :param labels: array-like with the true labels associated to each predicted label + """ + true_labels = labels + pred_labels = classif_predictions + self.cc = CC(self.classifier, fit_classifier=False) + self.confusion = BBSEhard.getConfusionJointProb(self.classifier.classes_, true_labels, pred_labels) + self.train_prevalence_ = F.prevalence_from_labels(labels, classes=self.classifier.classes_) + + @classmethod + def getConfusionJointProb(cls, classes, y, y_): + """ + Estimate the matrix with entry (i,j) being the estimate of :math:`P(\\hat{Y}=i,Y=j)` + + :param classes: array-like with the class names + :param y: array-like with the true labels + :param y_: array-like with the estimated labels + :return: np.ndarray + """ + conf = confusion_matrix(y, y_, labels=classes).T + joint_probs = conf / conf.sum() + return joint_probs + + def _weights_from_predictions(self, classif_predictions): + prevs_estim_cc = self.cc.aggregate(classif_predictions) + weights = BBSEhard.solve_importance_weights( + joint_probs=self.confusion, + unadjusted_counts=prevs_estim_cc, + solver=self.solver, + ) + return np.clip(weights, 0.0, None) + + @classmethod + def solve_importance_weights(cls, + joint_probs: np.ndarray, + unadjusted_counts: np.ndarray, + solver: Literal["minimize", "exact-raise", "exact-cc"]) -> np.ndarray: + """ + Function that tries to solve for :math:`p` the equation :math:`q = C w`, where :math:`q` is the vector of + `unadjusted counts` (as estimated, e.g., via classify and count) with :math:`q_i` an estimate of + :math:`P(\\hat{Y}=y_i)`, and where :math:`C` is the confusion matrix distribution with :math:`C_{ij}` an + estimate of :math:`P(\\hat{Y}=y_i,Y=y_j)`. + + :param joint_probs: array of shape `(n_classes, n_classes,)` with entry `(i,j)` being the estimate + of :math:`P(\\hat{Y}=y_i,Y=y_j)` + + :param unadjusted_counts: array of shape `(n_classes,)` containing the unadjusted prevalence values (e.g., as + estimated by CC or PCC) + + :param str solver: the method to use for solving the system of linear equations. Valid options are: + + * `exact-raise`: tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly lower than `n_classes`. + * `exact-cc`: if the matrix is not full rank, returns a vector of ones as the weights + * `minimize`: minimizes a loss, so the solution always exists + """ + A = np.asarray(joint_probs, dtype=float) + B = np.asarray(unadjusted_counts, dtype=float) + all_ones = np.full_like(B, fill_value=1., dtype=float) + + if solver in ["exact-raise", "exact-cc"]: + try: + return np.linalg.solve(A, B) + except np.linalg.LinAlgError: + if solver=='exact-cc': + return all_ones + else: + raise + elif solver == "minimize": + def loss(prev): + return np.linalg.norm(A @ prev - B) + + n_dims = len(all_ones) + bounds = [(0, np.inf)] * n_dims + return F.optim_minimize(loss, n_classes=A.shape[0], x0=all_ones, bounds=bounds, constraints=()) + else: + raise ValueError(f"unknown {solver=}") + + class PACC(AggregativeSoftQuantifier): """ @@ -615,7 +894,6 @@ class PACC(AggregativeSoftQuantifier): Raises an error if the matrix has rank strictly less than `n_classes`. * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) - * 'exact': deprecated, defaults to 'exact-cc' * 'minimize': minimizes the L2 norm of :math:`|Ax-B|`. This one generally works better, and is the default parameter. More details about this can be consulted in `Bunse, M. "On Multi-Class Extensions of Adjusted Classify and Count", on proceedings of the 2nd International Workshop on Learning to @@ -637,7 +915,7 @@ class PACC(AggregativeSoftQuantifier): classifier: BaseEstimator = None, fit_classifier=True, val_split=5, - solver: Literal['minimize', 'exact', 'exact-raise', 'exact-cc'] = 'minimize', + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', method: Literal['inversion', 'invariant-ratio'] = 'inversion', norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', n_jobs=None @@ -671,7 +949,7 @@ class PACC(AggregativeSoftQuantifier): def aggregate(self, classif_posteriors): prevs_estim = self.pcc.aggregate(classif_posteriors) - estimate = F.solve_adjustment( + estimate = ACC.solve_adjustment( class_conditional_rates=self.Pte_cond_estim_, unadjusted_counts=prevs_estim, solver=self.solver, @@ -693,7 +971,120 @@ class PACC(AggregativeSoftQuantifier): return confusion.T -class RLLS(AggregativeSoftQuantifier): +class BBSEsoft(ImportanceWeightQuantifier, AggregativeSoftQuantifier): + """ + `Black Box Shift Estimator` (BBSE) soft, the probabilistic variant of :class:`BBSEhard` that relies on the + posterior probabilities returned by a probabilistic classifier, instead of on crisp counts, to estimate the + joint distribution :math:`P(\\hat{Y}=i,Y=j)`. As in :class:`BBSEhard`, the sought importance weights + :math:`w_i=Q(i)/P(i)` (with :math:`P` and :math:`Q` the training and test distributions) are obtained by + solving :math:`q=C w`, with :math:`q` the (now probabilistic) prevalence vector estimated by PCC and + :math:`C` the joint-probability matrix. + + BBSE was proposed in + `Lipton, Z., Wang, Y. X., & Smola, A. (2018, July). + Detecting and correcting for label shift with black box predictors. + In International conference on machine learning (pp. 3122-3130). PMLR. + `_. + + :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be + the one indicated in `qp.environ['DEFAULT_CLS']` + + :param fit_classifier: whether to train the learner (default is True). Set to False if the + learner has been trained outside the quantifier. + + :param val_split: specifies the data used for generating classifier predictions. This specification + can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to + be extracted from the training set; or as an integer (default 5), indicating that the predictions + are to be generated in a `k`-fold cross-validation manner (with this integer indicating the value + for `k`); or as a tuple (X,y) defining the specific set of data to use for validation. + + :param str solver: indicates the method to use for solving the system of linear equations. Valid options are: + + * 'exact-raise': tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly less than `n_classes`. + * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to + no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) + * 'minimize': minimizes the squared L2 norm of :math:`|Ax-B|`. This one generally works better, and is the + default parameter. More details about this can be consulted in + `Tachet des Combes, R., Zhao, H., Wang, Y. X., & Gordon, G. J. (2020). + Domain adaptation with conditional distribution matching and generalized label shift. + Advances in Neural Information Processing Systems, 33, 19276-19289. + `_. + + :param str norm: the method to use for normalization. + + * `clip`, the values are clipped to the range [0,1] and then L1-normalized. + * `mapsimplex` projects vectors onto the probability simplex. This implementation relies on + `Mathieu Blondel's projection_simplex_sort `_ + * `condsoftmax`, applies a softmax normalization only to prevalence vectors that lie outside the simplex + + :param n_jobs: number of parallel workers + """ + + def __init__( + self, + classifier: BaseEstimator = None, + fit_classifier = True, + val_split = 5, + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', + norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', + n_jobs=None, + ): + super().__init__(classifier, fit_classifier, val_split) + self.n_jobs = qp._get_njobs(n_jobs) + self.solver = solver + self.norm = norm + + def _check_init_parameters(self): + if self.solver not in ACC.SOLVERS: + raise ValueError(f"unknown solver; valid ones are {ACC.SOLVERS}") + if self.norm not in ACC.NORMALIZATIONS: + raise ValueError(f"unknown normalization; valid ones are {ACC.NORMALIZATIONS}") + + def aggregation_fit(self, classif_predictions, labels): + """ + Estimates the joint distribution P(hat_Y,Y), using posterior probabilities in place of crisp counts. + + :param classif_predictions: array-like with posterior probabilities + :param labels: array-like with the true labels associated to each vector of posterior probabilities + """ + posteriors = classif_predictions + true_labels = labels + self.pcc = PCC(self.classifier, fit_classifier=False) + self.confusion = BBSEsoft.getConfusionJointProb(self.classifier.classes_, true_labels, posteriors) + self.train_prevalence_ = F.prevalence_from_labels(labels, classes=self.classifier.classes_) + + @classmethod + def getConfusionJointProb(cls, classes, y, posteriors): + """ + Estimate the matrix with entry (i,j) being the estimate of :math:`P(\\hat{Y}=i,Y=j)`, using the + posterior probabilities of the instances belonging to class :math:`j` in place of their hard predictions. + + :param classes: array-like with the class names + :param y: array-like with the true labels + :param posteriors: array-like of shape `(n_instances, n_classes,)` with posterior probabilities + :return: np.ndarray + """ + n_classes = len(classes) + joint_probs = np.zeros((n_classes, n_classes), dtype=float) + for j, class_ in enumerate(classes): + idx = y == class_ + if idx.any(): + joint_probs[:, j] = posteriors[idx].sum(axis=0) + joint_probs /= joint_probs.sum() + return joint_probs + + def _weights_from_predictions(self, classif_posteriors): + prevs_estim_pcc = self.pcc.aggregate(classif_posteriors) + weights = BBSEhard.solve_importance_weights( + joint_probs=self.confusion, + unadjusted_counts=prevs_estim_pcc, + solver=self.solver, + ) + return np.clip(weights, 0.0, None) + + +class RLLS(ImportanceWeightQuantifier, AggregativeSoftQuantifier): """ `Regularized Learning for Domain Adaptation under Label Shifts `_, used here as an aggregative @@ -750,7 +1141,6 @@ class RLLS(AggregativeSoftQuantifier): self.delta = delta self.clip_weights = clip_weights self.norm = norm - self.last_w_ = None def _check_init_parameters(self): _get_cvxpy() @@ -781,18 +1171,15 @@ class RLLS(AggregativeSoftQuantifier): self.pz_ = _rlls_predicted_marginal(classif_predictions, mode=self.mode) self.rho_ = _rlls_compute_3deltaC(len(self.classes_), len(labels), self.delta) - def aggregate(self, classif_posteriors): + def _weights_from_predictions(self, classif_posteriors): qz = _rlls_predicted_marginal(classif_posteriors, mode=self.mode) - w = _rlls_compute_weights( + return _rlls_compute_weights( self.C_zy_, qz, self.pz_, rho=self.alpha * self.rho_, clip=self.clip_weights, ) - self.last_w_ = w - estimate = self.train_prevalence_ * w - return F.normalize_prevalence(estimate, method=self.norm) class EMQ(AggregativeSoftQuantifier): @@ -1831,6 +2218,9 @@ KDEyML = _kdey.KDEyML KDEyHD = _kdey.KDEyHD KDEyCS = _kdey.KDEyCS +from . import _liep +LEIP = _liep.LEIP + # --------------------------------------------------------------- # aliases diff --git a/quapy/method/meta.py b/quapy/method/meta.py index 7d15b7c..58c26dd 100644 --- a/quapy/method/meta.py +++ b/quapy/method/meta.py @@ -37,6 +37,17 @@ if _histnet: else: HistNetQ = "HistNetQ is not available due to missing torch package" +try: + from . import _gmnet +except ModuleNotFoundError: + _gmnet = None + + +if _gmnet: + GMNet = _gmnet.GMNet +else: + GMNet = "GMNet is not available due to missing torch and/or geotorch packages" + class MedianEstimator(BinaryQuantifier): """ diff --git a/quapy/tests/test_labelshift.py b/quapy/tests/test_labelshift.py new file mode 100644 index 0000000..adaae03 --- /dev/null +++ b/quapy/tests/test_labelshift.py @@ -0,0 +1,84 @@ +import unittest + +import numpy as np +from sklearn.linear_model import LogisticRegression +from sklearn.exceptions import NotFittedError +from sklearn.utils.validation import check_is_fitted + +from quapy.classification.labelshift import LabelShiftedClassifier +from quapy.method.aggregative import ACC, BBSEhard, BBSEsoft, ImportanceWeightQuantifier, RLLS +from quapy.functional import check_prevalence_vector +from quapy.tests._synthetic import make_dataset + + +class TestLabelShift(unittest.TestCase): + + dataset = make_dataset(n_train=200, n_test=100, n_classes=3, n_features=10, random_state=7, name='labelshift') + + def test_importance_weight_quantifiers(self): + Xtr, ytr = self.dataset.training.Xy + Xte = self.dataset.test.X + for cls in [BBSEhard, BBSEsoft, RLLS]: + if cls is RLLS: + try: + import cvxpy # noqa: F401 + except ImportError: + continue + + q = cls(LogisticRegression(max_iter=2000), val_split=3) + q.fit(Xtr, ytr) + + weights = q.get_importance_weights(Xte) + prevalence, weights2 = q.quantify_and_weigh(Xte) + prevalence_direct = q.predict(Xte) + + self.assertTrue(check_prevalence_vector(prevalence)) + np.testing.assert_allclose(weights, weights2) + np.testing.assert_allclose(prevalence, prevalence_direct) + # no leftover mutable per-call state (safety against races under concurrent calls) + self.assertFalse(hasattr(q, 'weights')) + self.assertFalse(hasattr(q, 'last_w_')) + + def test_label_shifted_classifier_generic_quantifier(self): + Xtr, ytr = self.dataset.training.Xy + Xte = self.dataset.test.X + + base_classifier = LogisticRegression(max_iter=2000) + q = ACC(LogisticRegression(max_iter=2000), val_split=3) + self.assertNotIsInstance(q, ImportanceWeightQuantifier) + + lsc = LabelShiftedClassifier(base_classifier, q) + lsc.fit(Xtr, ytr) + + weights = lsc.get_importance_weights(Xte) + self.assertEqual(weights.shape, (self.dataset.training.n_classes,)) + self.assertTrue(np.all(weights >= lsc.weight_epsilon)) + + preds = lsc.predict(Xte) + probs = lsc.predict_proba(Xte) + self.assertEqual(len(preds), len(Xte)) + self.assertEqual(probs.shape, (len(Xte), self.dataset.training.n_classes)) + np.testing.assert_allclose(probs.sum(axis=1), 1.0, rtol=1e-5) + + # the classifier instance passed at construction time must be left untouched + with self.assertRaises(NotFittedError): + check_is_fitted(base_classifier) + + def test_label_shifted_classifier_importance_weight_quantifier(self): + Xtr, ytr = self.dataset.training.Xy + Xte = self.dataset.test.X + + q = BBSEhard(LogisticRegression(max_iter=2000), val_split=3) + self.assertIsInstance(q, ImportanceWeightQuantifier) + + lsc = LabelShiftedClassifier(LogisticRegression(max_iter=2000), q) + lsc.fit(Xtr, ytr) + + weights = lsc.get_importance_weights(Xte) + adapted = lsc.get_classifier(Xte) + check_is_fitted(adapted) # should not raise + self.assertEqual(set(adapted.get_params()['class_weight'].keys()), set(lsc.classes_)) + + +if __name__ == '__main__': + unittest.main() diff --git a/quapy/tests/test_methods.py b/quapy/tests/test_methods.py index dfd3b12..0e41ff5 100644 --- a/quapy/tests/test_methods.py +++ b/quapy/tests/test_methods.py @@ -8,7 +8,7 @@ from sklearn.linear_model import LogisticRegression from quapy.method import AGGREGATIVE_METHODS, BINARY_METHODS, NON_AGGREGATIVE_METHODS from quapy.method.non_aggregative import DMx, EDx, HDx -from quapy.method.aggregative import ACC, DMy, EDy, KDEyCS, RLLS +from quapy.method.aggregative import ACC, BBSEhard, BBSEsoft, DMy, EDy, KDEyCS, LEIP, RLLS from quapy.method.meta import Ensemble from quapy.functional import check_prevalence_vector from quapy.tests._synthetic import make_dataset @@ -21,6 +21,7 @@ OPTIONAL_AGGREGATIVE_METHODS = { 'PQ', 'RLLS', 'EDy', + 'LEIP', } OPTIONAL_NON_AGGREGATIVE_METHODS = { @@ -126,6 +127,38 @@ class TestMethods(unittest.TestCase): estim_prevalences2 = model2.predict(dataset.test.X) self.assertTrue(check_prevalence_vector(estim_prevalences2)) + def test_gmnet(self): + try: + import torch + import geotorch + except ModuleNotFoundError: + print('the torch and/or geotorch packages are not installed; skipping unit test for GMNet') + return + + from quapy.method.meta import GMNet + from quapy.protocol import UPP + + for dataset in TestMethods.datasets: + # single GM layer, no CKA regularization + model = GMNet( + bag_size=20, n_bags_train=10, n_bags_val=5, train_epochs=2, patience=1, batch_size=2, + device='cpu', checkpointdir='./checkpoint_test_gmnet' + ) + model.fit(*dataset.training.Xy) + estim_prevalences = model.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + + # multiple GM layers + CKA regularization, and fit_from_samples + given_samples = UPP(dataset.training, sample_size=20, repeats=8, random_state=1) + val_samples = UPP(dataset.training, sample_size=20, repeats=4, random_state=2) + model2 = GMNet( + n_gm_layers=2, num_gaussians=3, gaussian_dimensions=4, cka_regularization=0.1, + bag_size=20, train_epochs=2, patience=1, batch_size=2, device='cpu', + checkpointdir='./checkpoint_test_gmnet' + ) + model2.fit_from_samples(given_samples, val_protocol=val_samples, mix_bags=True) + estim_prevalences2 = model2.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences2)) def test_composable(self): try: @@ -178,6 +211,29 @@ class TestMethods(unittest.TestCase): self.assertTrue(check_prevalence_vector(estim_prevalences)) + def test_leip(self): + dataset = TestMethods.tiny_dataset_multiclass + q = LEIP(LogisticRegression(max_iter=2000), val_split=3) + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + + def test_leip_fixed_tau(self): + dataset = TestMethods.tiny_dataset_binary + q = LEIP(LogisticRegression(max_iter=2000), val_split=None, tau=0.6) + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + + def test_bbse(self): + dataset = TestMethods.tiny_dataset_multiclass + for cls in (BBSEhard, BBSEsoft): + for solver in ('minimize', 'exact-raise', 'exact-cc'): + q = cls(LogisticRegression(max_iter=2000), val_split=3, solver=solver) + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + def test_edy(self): try: import quadprog # noqa: F401 diff --git a/setup.py b/setup.py index 22c424f..3ff18da 100644 --- a/setup.py +++ b/setup.py @@ -131,7 +131,7 @@ setup( # projects. extras_require={ # Optional 'bayes': ['jax', 'jaxlib', 'numpyro', 'pystan', 'setuptools<82'], - 'neural': ['torch'], + 'neural': ['torch', 'geotorch'], 'tests': ['certifi'], 'docs' : ['pydata-sphinx-theme', 'myst-parser', 'sphinx-design'], },