Fix broken imports left over from the BBSE/LEIP commit

quapy/classification/__init__.py imported a labelshift module that doesn't exist
anywhere in the repo (and nothing else references it), breaking `import quapy`
entirely; the import is removed.

quapy/method/aggregative.py imported a nonexistent `_liep` module instead of the
actual file, _liep_draft.py, which itself never got a LEIP class added (only
helper functions such as leip()). Guards the import so LEIP degrades to an
"not available" placeholder, like the other optional neural methods, instead of
crashing the whole package. test_leip/test_leip_fixed_tau still fail as a result;
finishing LEIP is left for a follow-up.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

CHANGE_LOG.txt

@@ -1,3 +1,10 @@
+Change Log 0.2.2
+-----------------
+
+- Added HistNetQ, based on the original implementation https://github.com/pglez84/histnetq
+
+- Minor fixes
+
 Change Log 0.2.1
 -----------------
 

TODO.txt

@@ -1,14 +1,7 @@
-Adapt examples; remaining: example 4-onwards
-not working: 15 (qunfold)
-
-Solve the warnings issue; right now there is a warning ignore in method/__init__.py:
-
 Add 'platt' to calib options in EMQ?
 
 Allow n_prevpoints in APP to be specified by a user-defined grid?
 
-Update READMEs, wiki, & examples for new fit-predict interface
-
 Add the fix suggested by Alexander:
 
 For a more general application, I would maybe first establish a per-class threshold value of plausible prevalence
@@ -20,17 +13,14 @@ scale each value by per-class thresholds, i.e., [0.33*0.1, 0.33*1, 0.33*1]/sum.
 - This affects LabelledCollection
 - This functionality should be accessible via sampling protocols and evaluation functions
 
-Solve the pre-trained classifier issues. An example is the coptic-codes script I did, which needed a mock_lr to
-work for having access to classes_; think also the case in which the precomputed outputs are already generated
-as in the unifying problems code.
-
-
-- [TODO] document confidence in manuals
 - [TODO] Test the return_type="index" in protocols and finish the "distributing_samples.py" example
 - [TODO] add ensemble methods SC-MQ, MC-SQ, MC-MQ
-- [TODO] add HistNetQ
 - [TODO] add CDE-iteration and Bayes-CDE methods
 - [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
+- [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)

docs/source/manuals.rst

@@ -7,6 +7,7 @@ Manuals
 
    manuals/datasets
    manuals/evaluation
+   manuals/label-shift-adaptation
    manuals/methods
    manuals/model-selection
    manuals/plotting

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 <manuals/label-shift-adaptation:Label Shift Adaptation>` 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 <manuals/label-shift-adaptation:Label Shift Adaptation>` 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)
@@ -851,7 +905,36 @@ model.fit(*dataset.training.Xy)
 estim_prevalence = model.predict(dataset.test.X)
 ```
 
-(confidence-regions-for-class-prevalence-estimation)=
+### HistNetQ
+
+QuaPy offers an implementation of 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).
+Quantification using permutation-invariant networks based on histograms.
+Neural Computing and Applications, 37(5), 3505-3520._](https://doi.org/10.1007/s00521-024-10721-1)
+
+This model requires `torch` to be installed. Like QuaNet, HistNetQ is trained end-to-end on
+samples ("bags") of known prevalence rather than on individually labeled instances; unlike QuaNet,
+it requires no classifier at all, only an optional feature extraction module (a plain identity
+module is used by default, for already-vectorized data).
+
+```python
+import quapy as qp
+from quapy.method.meta import HistNetQ
+
+dataset = qp.datasets.fetch_UCIBinaryDataset('haberman')
+
+model = HistNetQ(bag_size=100, device='cpu')
+model.fit(*dataset.training.Xy)
+estim_prevalence = model.predict(dataset.test.X)
+```
+
+HistNetQ can alternatively be trained directly from a protocol that already provides the training
+samples (e.g., when only bag-level prevalence values are available), via the `fit_from_samples`
+method; see the API documentation for further details.
+
+
 ## Quantifiers with Uncertainty Quantification
 
 _(New in v0.2.0!)_ Some quantification methods go beyond providing a single point estimate of class prevalence values and also produce confidence regions, which characterize the uncertainty around the point estimate. In QuaPy, two such families are currently implemented: bootstrap methods and Bayesian methods.

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
 -----------------------------------
 

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 <https://jmlr.org/papers/v20/18-456.html>`_,
-          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
 # ------------------------------------------------------------------------------------------

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

quapy/method/_gmnet.py

@@ -0,0 +1,328 @@
+"""
+GMNet implementation.
+
+Ported from the reference implementation at https://github.com/pglez84/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 <https://github.com/pglez84/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.

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,84 +81,22 @@ 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):
-    """
-    :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):
+class HistNetQ(BagTrainedQuantifier):
     """
     Implementation of `HistNetQ <https://github.com/pglez84/histnetq>`_, a neural network for
     quantification that learns a differentiable histogram-based representation of a sample, trained
@@ -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)

quapy/method/_liep_draft.py

@@ -0,0 +1,261 @@
+import numpy as np
+
+
+def normalize(v, eps=1e-12):
+    v = np.asarray(v, dtype=float)
+    v = np.maximum(v, 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) ∝ p_old(y|x) * new_prior(y) / source_prior(y)
+
+    Parameters
+    ----------
+    probs : array, shape (n_samples, n_classes) or (n_classes,)
+        Probabilistic outputs of the classifier.
+    new_prior : array, shape (n_classes,)
+        Current estimate of the target prior.
+    source_prior : array, shape (n_classes,)
+        Source/training prior.
+    """
+    probs = np.asarray(probs, dtype=float)
+    one_dim = probs.ndim == 1
+    if one_dim:
+        probs = probs[None, :]
+
+    new_prior = normalize(new_prior, eps=eps)
+    source_prior = normalize(source_prior, eps=eps)
+
+    weights = new_prior / np.maximum(source_prior, eps)
+    updated = probs * weights[None, :]
+    updated = updated / np.maximum(updated.sum(axis=1, keepdims=True), eps)
+
+    return updated[0] if one_dim else updated
+
+
+def confusion_statistic(y_true, y_pred, n_classes, mode="recall", eps=1e-12):
+    """
+    Computes the minimum diagonal statistic used to choose tau.
+
+    mode="recall":
+        diag(C) / row sums, i.e. per-true-class recall.
+
+    mode="precision":
+        diag(C) / column sums, i.e. per-predicted-class precision.
+
+    The paper calls this quantity 'minimum recall', but also describes
+    a column-normalized confusion matrix, which would correspond more
+    closely to precision. I expose both options.
+    """
+    y_true = np.asarray(y_true, dtype=int)
+    y_pred = np.asarray(y_pred, dtype=int)
+
+    C = np.zeros((n_classes, n_classes), dtype=float)
+    for yt, yp in zip(y_true, y_pred):
+        C[yt, yp] += 1.0
+
+    if mode == "recall":
+        denom = C.sum(axis=1)
+    elif mode == "precision":
+        denom = C.sum(axis=0)
+    else:
+        raise ValueError("mode must be 'recall' or 'precision'.")
+
+    diag = np.diag(C)
+    valid = denom > 0
+
+    if not np.any(valid):
+        raise ValueError("No valid classes found in the confusion matrix.")
+
+    scores = diag[valid] / np.maximum(denom[valid], eps)
+    return float(np.min(scores))
+
+
+def choose_tau_from_validation(
+    target_probs,
+    val_probs,
+    y_val,
+    mode="recall",
+    retain_fraction=None,
+):
+    """
+    Chooses tau following the spirit of Section 4.1.
+
+    The paper says tau is selected from the top n percentile of target
+    confidences, with n = min-recall * 100. To make the selected set A
+    contain approximately min-recall fraction of the target samples, we set:
+
+        tau = quantile(max_probs, 1 - min_recall)
+
+    so that about min_recall of the target points satisfy max_prob >= tau.
+
+    If retain_fraction is provided, it overrides the validation-derived value.
+    """
+    target_probs = np.asarray(target_probs, dtype=float)
+    val_probs = np.asarray(val_probs, dtype=float)
+    y_val = np.asarray(y_val, dtype=int)
+
+    n_classes = target_probs.shape[1]
+    y_val_pred = val_probs.argmax(axis=1)
+
+    if retain_fraction is None:
+        retain_fraction = confusion_statistic(
+            y_true=y_val,
+            y_pred=y_val_pred,
+            n_classes=n_classes,
+            mode=mode,
+        )
+
+    retain_fraction = float(np.clip(retain_fraction, 0.0, 1.0))
+
+    target_conf = target_probs.max(axis=1)
+
+    if retain_fraction <= 0:
+        tau = np.inf
+    elif retain_fraction >= 1:
+        tau = -np.inf
+    else:
+        tau = np.quantile(target_conf, 1.0 - retain_fraction)
+
+    return tau, retain_fraction
+
+
+def leip(
+    target_probs,
+    source_prior,
+    tau=None,
+    val_probs=None,
+    y_val=None,
+    threshold_mode="recall",
+    count_smoothing=0.0,
+    eps=1e-12,
+    return_details=False,
+):
+    """
+    LEIP: Label shift Estimation with Incremental Prior update.
+
+    Parameters
+    ----------
+    target_probs : array, shape (n_target, n_classes)
+        Probabilistic classifier outputs on the target/test set.
+    source_prior : array, shape (n_classes,)
+        Source class prior p_s(y).
+    tau : float or None
+        Confidence threshold. If None, it is estimated using validation data.
+    val_probs : array, shape (n_val, n_classes), optional
+        Validation probabilistic outputs, required if tau is None.
+    y_val : array, shape (n_val,), optional
+        Validation labels, required if tau is None.
+    threshold_mode : {"recall", "precision"}
+        Statistic used to choose tau from validation data.
+    count_smoothing : float
+        Optional additive smoothing for pseudo-label counts. Set to 0.0
+        for a closer implementation of the paper; use a small value such
+        as 1e-8 for extra numerical robustness.
+    eps : float
+        Numerical stabilizer.
+    return_details : bool
+        If True, returns diagnostic information.
+
+    Returns
+    -------
+    estimated_prior : array, shape (n_classes,)
+        Estimated target class distribution.
+    details : dict, optional
+        Returned only if return_details=True.
+    """
+    target_probs = np.asarray(target_probs, dtype=float)
+    if target_probs.ndim != 2:
+        raise ValueError("target_probs must have shape (n_samples, n_classes).")
+
+    n_target, n_classes = target_probs.shape
+    source_prior = normalize(source_prior, eps=eps)
+
+    if n_classes != len(source_prior):
+        raise ValueError("source_prior must have one entry per class.")
+
+    # Step 1: choose tau if needed
+    if tau is None:
+        if val_probs is None or y_val is None:
+            raise ValueError("val_probs and y_val are required when tau is None.")
+
+        tau, retain_fraction = choose_tau_from_validation(
+            target_probs=target_probs,
+            val_probs=val_probs,
+            y_val=y_val,
+            mode=threshold_mode,
+        )
+    else:
+        retain_fraction = None
+
+    target_conf = target_probs.max(axis=1)
+    target_top = target_probs.argmax(axis=1)
+
+    # Step 2: high-confidence set A
+    A_mask = target_conf >= tau
+    A_labels = target_top[A_mask]
+
+    counts = np.full(n_classes, count_smoothing, dtype=float)
+
+    if len(A_labels) > 0:
+        counts += np.bincount(A_labels, minlength=n_classes)
+        current_prior = counts / counts.sum()
+    else:
+        # Fallback if tau is too strict.
+        # One could also use classify-and-count over the full target set.
+        current_prior = source_prior.copy()
+
+    # Step 3: low-confidence set B, sorted by decreasing confidence
+    B_indices = np.where(~A_mask)[0]
+    B_indices = B_indices[np.argsort(-target_conf[B_indices])]
+
+    # Step 4: incremental pass over B
+    incremental_labels = []
+
+    for idx in B_indices:
+        corrected = bayes_prior_update(
+            probs=target_probs[idx],
+            new_prior=current_prior,
+            source_prior=source_prior,
+            eps=eps,
+        )
+        pseudo_label = int(np.argmax(corrected))
+        incremental_labels.append(pseudo_label)
+
+        counts[pseudo_label] += 1.0
+        current_prior = counts / counts.sum()
+
+    estimated_intermediate_prior = current_prior.copy()
+
+    # Step 5: final complete pass over all target instances
+    corrected_all = bayes_prior_update(
+        probs=target_probs,
+        new_prior=estimated_intermediate_prior,
+        source_prior=source_prior,
+        eps=eps,
+    )
+
+    final_labels = corrected_all.argmax(axis=1)
+    estimated_prior = np.bincount(final_labels, minlength=n_classes).astype(float)
+    estimated_prior /= estimated_prior.sum()
+
+    if not return_details:
+        return estimated_prior
+
+    details = {
+        "tau": tau,
+        "retain_fraction": retain_fraction,
+        "n_A": int(A_mask.sum()),
+        "n_B": int((~A_mask).sum()),
+        "A_mask": A_mask,
+        "intermediate_prior": estimated_intermediate_prior,
+        "final_labels": final_labels,
+        "corrected_probs": corrected_all,
+        "incremental_labels": np.asarray(incremental_labels, dtype=int),
+    }
+
+    return estimated_prior, details

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()

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 <https://jmlr.org/papers/v20/18-456.html>`_,
           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 <https://jmlr.org/papers/v20/18-456.html>`_,
+              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.
+    <https://proceedings.mlr.press/v80/lipton18a.html>`_.
+
+    :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.
+          <https://proceedings.neurips.cc/paper_files/paper/2020/hash/dfbfa7ddcfffeb581f50edcf9a0204bb-Abstract.html>`_.
+
+    :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 <https://gist.github.com/mblondel/6f3b7aaad90606b98f71>`_
+        * `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.
+    <https://proceedings.mlr.press/v80/lipton18a.html>`_.
+
+    :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.
+          <https://proceedings.neurips.cc/paper_files/paper/2020/hash/dfbfa7ddcfffeb581f50edcf9a0204bb-Abstract.html>`_.
+
+    :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 <https://gist.github.com/mblondel/6f3b7aaad90606b98f71>`_
+        * `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
     <https://arxiv.org/abs/1903.09734>`_, 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,12 @@ KDEyML = _kdey.KDEyML
 KDEyHD = _kdey.KDEyHD
 KDEyCS = _kdey.KDEyCS
 
+try:
+    from . import _liep_draft as _liep
+    LEIP = _liep.LEIP
+except AttributeError:
+    LEIP = "LEIP is not available (incomplete implementation in _liep_draft.py)"
+
 
 # ---------------------------------------------------------------
 # aliases

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):
     """

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

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'],
     },