diff --git a/CHANGE_LOG.txt b/CHANGE_LOG.txt index 836b0f9..a384d39 100644 --- a/CHANGE_LOG.txt +++ b/CHANGE_LOG.txt @@ -29,6 +29,9 @@ Change Log 0.2.1 - Deep code revision and improved codebase +- Added EDx/EDy from quantificationlib (thanks to Pablo and Juanjo!) + + Change Log 0.2.0 ----------------- diff --git a/docs/source/manuals.rst b/docs/source/manuals.rst index a426786..97e0628 100644 --- a/docs/source/manuals.rst +++ b/docs/source/manuals.rst @@ -7,7 +7,6 @@ Manuals manuals/datasets manuals/evaluation - manuals/explicit-loss-minimization manuals/methods manuals/model-selection manuals/plotting diff --git a/docs/source/manuals/explicit-loss-minimization.md b/docs/source/manuals/explicit-loss-minimization.md deleted file mode 100644 index f80c434..0000000 --- a/docs/source/manuals/explicit-loss-minimization.md +++ /dev/null @@ -1,26 +0,0 @@ -# Explicit Loss Minimization - -QuaPy makes available several Explicit Loss Minimization (ELM) methods, including -SVM(Q), SVM(KLD), SVM(NKLD), SVM(AE), or SVM(RAE). -These methods require to first download the -[svmperf](http://www.cs.cornell.edu/people/tj/svm_light/svm_perf.html) -package, apply the patch -[svm-perf-quantification-ext.patch](https://github.com/HLT-ISTI/QuaPy/blob/master/svm-perf-quantification-ext.patch), and compile the sources. -The script [prepare_svmperf.sh](https://github.com/HLT-ISTI/QuaPy/blob/master/prepare_svmperf.sh) does all the job. Simply run: - -``` -./prepare_svmperf.sh -``` - -The resulting directory `svm_perf_quantification/` contains the -patched version of _svmperf_ with quantification-oriented losses. - -The [svm-perf-quantification-ext.patch](https://github.com/HLT-ISTI/QuaPy/blob/master/prepare_svmperf.sh) is an extension of the patch made available by -[Esuli et al. 2015](https://dl.acm.org/doi/abs/10.1145/2700406?casa_token=8D2fHsGCVn0AAAAA:ZfThYOvrzWxMGfZYlQW_y8Cagg-o_l6X_PcF09mdETQ4Tu7jK98mxFbGSXp9ZSO14JkUIYuDGFG0) -that allows SVMperf to optimize for -the _Q_ measure as proposed by [Barranquero et al. 2015](https://www.sciencedirect.com/science/article/abs/pii/S003132031400291X) -and for the _KLD_ and _NKLD_ measures as proposed by [Esuli et al. 2015](https://dl.acm.org/doi/abs/10.1145/2700406?casa_token=8D2fHsGCVn0AAAAA:ZfThYOvrzWxMGfZYlQW_y8Cagg-o_l6X_PcF09mdETQ4Tu7jK98mxFbGSXp9ZSO14JkUIYuDGFG0). -This patch extends the above one by also allowing SVMperf to optimize for -_AE_ and _RAE_. -See the [](./methods) manual for more details and code examples. - diff --git a/docs/source/manuals/methods.md b/docs/source/manuals/methods.md index 7851362..4d7d75e 100644 --- a/docs/source/manuals/methods.md +++ b/docs/source/manuals/methods.md @@ -488,66 +488,37 @@ the last two methods (SVM(AE) and SVM(RAE)) have been implemented in QuaPy in order to make available ELM variants for what nowadays are considered the most well-behaved evaluation metrics in quantification. -In order to make these models work, you would need to run the script -`prepare_svmperf.sh` (distributed along with QuaPy) that -downloads `SVMperf`' source code, applies a patch that -implements the quantification oriented losses, and compiles the -sources. +#### Installing the SVMperf backend -If you want to add any custom loss, you would need to modify -the source code of `SVMperf` in order to implement it, and -assign a valid loss code to it. Then you must re-compile -the whole thing and instantiate the quantifier in QuaPy -as follows: +These methods rely on Joachim's [SVMperf](https://www.cs.cornell.edu/people/tj/svm_light/svm_perf.html), +patched with quantification-oriented losses. QuaPy provides the script +[`prepare_svmperf.sh`](https://github.com/HLT-ISTI/QuaPy/blob/master/prepare_svmperf.sh), +which downloads the original sources, applies the patch, and compiles the +resulting binary. In practice, this amounts to running: -```python -# you can either set the path to your custom svm_perf_quantification implementation -# in the environment variable, or as an argument to the constructor of ELM -qp.environ['SVMPERF_HOME'] = './path/to/svm_perf_quantification' - -# assign an alias to your custom loss and the id you have assigned to it -svmperf = qp.classification.svmperf.SVMperf -svmperf.valid_losses['mycustomloss'] = 28 - -# instantiate the ELM method indicating the loss -model = qp.method.aggregative.ELM(loss='mycustomloss') +```sh +./prepare_svmperf.sh ``` -All ELM are binary quantifiers since they rely on `SVMperf`, that -currently supports only binary classification. -ELM variants (any binary quantifier in general) can be extended -to operate in single-label scenarios trivially by adopting a -"one-vs-all" strategy (as, e.g., in -[_Gao, W. and Sebastiani, F. (2016). From classification to quantification in tweet sentiment -analysis. Social Network Analysis and Mining, 6(19):1–22_](https://link.springer.com/article/10.1007/s13278-016-0327-z)). -In QuaPy this is possible by using the `OneVsAll` class. - -There are two ways for instantiating this class, `OneVsAllGeneric` that works for -any quantifier, and `OneVsAllAggregative` that is optimized for aggregative quantifiers. -In general, you can simply use the `newOneVsAll` function and QuaPy will choose -the more convenient of the two. +This creates a directory `svm_perf_quantification/`. Once this is available, +you can point QuaPy to it with: ```python -import quapy as qp -from quapy.method.aggregative import SVMQ - -# load a single-label dataset (this one contains 3 classes) -train, test = qp.datasets.fetch_twitter('hcr', pickle=True).train_test - -# let qp know where svmperf is -qp.environ['SVMPERF_HOME'] = '../svm_perf_quantification' - -model = newOneVsAll(SVMQ(), n_jobs=-1) # run them on parallel -model.fit(*train.Xy) -estim_prevalence = model.predict(test.X) +qp.environ['SVMPERF_HOME'] = './svm_perf_quantification' ``` -Check the examples on [explicit loss minimization](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/17.explicit_loss_minimization.py) -and on [one versus all quantification](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/10.one_vs_all.py) for more details. -**Note** that the _one versus all_ approach is considered inappropriate under prior probability shift, though. - - +The patch extends the one originally released for +[Esuli and Sebastiani (2015)](https://dl.acm.org/doi/abs/10.1145/2700406) +and also covers the `Q`, `AE`, and `RAE` losses used by QuaPy's ELM wrappers. +All ELM methods are binary because `SVMperf` itself is binary. They can still +be wrapped in a one-vs-all scheme for single-label multiclass problems, though +this strategy is generally considered inappropriate under prior probability +shift. See the examples on +[explicit loss minimization](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/17.explicit_loss_minimization.py) +and on +[one versus all quantification](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/10.one_vs_all.py) +for minimal working code. ## Non-Aggregative Methods @@ -635,6 +606,18 @@ estim_prevalence_hdx = hdx.predict(test.X) Note that, unlike HDy, HDx requires no classifier whatsoever, since it operates directly on the covariates. +### Energy Distance x (EDx) + +QuaPy also provides `qp.method.non_aggregative.EDx`, which is the +feature-space counterpart of `EDy`: it keeps the same energy-distance +formulation and quadratic-program solver, but applies them directly to the raw +instances instead of first projecting them onto posterior probabilities through +a classifier. In this sense, `EDx` is to `EDy` what `DMx` is to `DMy`. + +`EDx` works for binary and multiclass problems, accepts the same `distance` +options as `EDy` (`'manhattan'`, `'euclidean'`, or a custom callable), and +requires the optional dependency `pip install quadprog`. + ### ReadMe `ReadMe` is a non-aggregative quantification method proposed by @@ -1026,7 +1009,8 @@ estim_prevalence, conf_region = model.predict_conf(test.X) #### BayesianKDEy (a Bayesian implementation of KDEyML) `BayesianKDEy`, available at `qp.method._bayesian.BayesianKDEy`, is a Bayesian -version of KDEy. Instead of solving for the single prevalence vector that +version of KDEy proposed by [Moreo et al. 2026](https://arxiv.org/abs/2607.04977). +Instead of solving for the single prevalence vector that minimizes a divergence between the test distribution and a KDE-based mixture model (as the KDEy variants above do), `BayesianKDEy` places a Dirichlet prior over the prevalence vector and samples its posterior via Markov Chain diff --git a/docs/source/manuals/protocols.md b/docs/source/manuals/protocols.md index aae89cb..7129f83 100644 --- a/docs/source/manuals/protocols.md +++ b/docs/source/manuals/protocols.md @@ -33,7 +33,7 @@ QuaPy provides implementations of most popular sample generation protocols used in literature. This is the subject of the following sections. -## Artificial-Prevalence Protocol +## APP: Artificial-Prevalence Protocol The "artificial-sampling protocol" (APP) proposed by [Forman (2005)](https://link.springer.com/chapter/10.1007/11564096_55) @@ -138,7 +138,7 @@ Each point corresponds to one sampled prevalence vector. As expected, the points lie on a regular grid over the simplex, ensuring systematic coverage of the prevalence space. -## Sampling from the unit-simplex, the Uniform-Prevalence Protocol (UPP) +## UPP: Sampling from the unit-simplex, the Uniform-Prevalence Protocol Generating all possible combinations from a grid of prevalence values (APP) in multiclass is cumbersome, and when the number of classes increases it rapidly @@ -177,6 +177,10 @@ regular lattice. Instead, it spreads samples over the simplex in a statistically uniform way, making it attractive when the number of classes is large and exhaustive grids become impractical. +*Note* that UPP is actually a different (modern) implementation of the Artificial Prevalence Protocol, +and is here given a different name simply to allow both implementations coexist in QuaPy. "UPP" is not a +proper accademic name, and practitioners should rather refer to it as APP. + ## Natural-Prevalence Protocol The "natural-prevalence protocol" (NPP) comes down to generating samples drawn diff --git a/quapy/method/__init__.py b/quapy/method/__init__.py index 1273a5a..ff6a769 100644 --- a/quapy/method/__init__.py +++ b/quapy/method/__init__.py @@ -63,11 +63,13 @@ MULTICLASS_METHODS = { confidence.BayesianCC, _bayesian.BayesianKDEy, _bayesian.BayesianMAPLS, + non_aggregative.EDx, } NON_AGGREGATIVE_METHODS = { non_aggregative.MaximumLikelihoodPrevalenceEstimation, - non_aggregative.DMx + non_aggregative.DMx, + non_aggregative.EDx } META_METHODS = { diff --git a/quapy/method/_edy.py b/quapy/method/_edy.py deleted file mode 100644 index 9d0f938..0000000 --- a/quapy/method/_edy.py +++ /dev/null @@ -1,269 +0,0 @@ -from typing import Callable, Union - -import numpy as np -from sklearn.base import BaseEstimator -from sklearn.metrics.pairwise import euclidean_distances, manhattan_distances - -import quapy as qp -import quapy.functional as F -from quapy.method._helper import _get_quadprog -from quapy.method.aggregative import AggregativeSoftQuantifier - - -class EDy(AggregativeSoftQuantifier): - """ - Energy Distance y (EDy), a posterior-space distribution-matching quantifier - based on energy distance. - - The method represents each class by the posterior-probability vectors - produced by a probabilistic classifier on validation data, and estimates the - test prevalence vector by matching the test posterior distribution against - the class-conditional validation distributions through an energy-distance - objective solved as a quadratic program. The method is therefore another - instance of the general mixture-matching view of quantification, but it - operates directly on posterior vectors rather than on histogram summaries. - - This implementation works for binary and multiclass single-label - quantification and relies on the optional ``quadprog`` dependency. It was - adapted to QuaPy's current aggregative API from the original implementation - available in `quantificationlib `_. - - The current implementation follows the energy-distance formulation discussed - in: - - * Alberto Castaño, Laura Morán-Fernández, Jaime Alonso, - Verónica Bolón-Canedo, Amparo Alonso-Betanzos, and Juan José del Coz. - *An analysis of quantification methods based on matching distributions*. - * Hideko Kawakubo, Marthinus Christoffel du Plessis, and Masashi Sugiyama - (2016). *Computationally efficient class-prior estimation under class - balance change using energy distance*. IEICE Transactions on Information - and Systems, 99(1):176-186. - - :param classifier: a scikit-learn ``BaseEstimator``, or ``None`` to use - ``qp.environ['DEFAULT_CLS']`` - :param fit_classifier: whether to train the learner (default ``True``). - Set to ``False`` if the learner has already been trained outside the - quantifier - :param val_split: specification of the data used for generating validation - posterior probabilities. This can be an integer (default ``5``) for - k-fold cross-validation, a float in ``(0, 1)`` for a held-out split, - or a tuple ``(X, y)`` with explicit validation data - :param distance: distance used to compare posterior vectors. Valid string - aliases are ``'manhattan'`` (default) and ``'euclidean'``; a custom - callable compatible with pairwise-distance signatures can also be used - :param n_jobs: number of parallel workers (default ``None``, meaning the - value is taken from the environment) - """ - - def __init__( - self, - classifier: BaseEstimator = None, - fit_classifier: bool = True, - val_split=5, - distance: Union[str, Callable] = 'manhattan', - n_jobs=None, - ): - super().__init__(classifier, fit_classifier, val_split) - self.distance = distance - self.n_jobs = qp._get_njobs(n_jobs) - self.train_n_cls_i_ = None - self.train_distrib_ = None - self.K_ = None - self.G_ = None - self.C_ = None - self.b_ = None - self.a_ = None - - def _check_init_parameters(self): - self.distance = self._resolve_distance_function(self.distance) - - @staticmethod - def _resolve_distance_function(distance): - if isinstance(distance, str): - if distance == 'manhattan': - return manhattan_distances - if distance == 'euclidean': - return euclidean_distances - raise ValueError( - f"unknown distance {distance!r}; valid aliases are 'manhattan' and 'euclidean'" - ) - if not hasattr(distance, '__call__'): - raise ValueError('distance must be a valid string alias or a callable function') - return distance - - def _is_pd(self, m): - """Check whether a symmetric matrix is positive definite. - - This helper is used before invoking ``quadprog`` because the quadratic - term of the optimization problem must be positive definite. - """ - return self._dpofa(m)[0] == 0 - - def _dpofa(self, m): - """Factor a symmetric positive definite matrix. - - This is a lightweight Python adaptation of the ``dpofa`` routine used by - ``quadprog``. Here it is mainly employed as a numerical check while - preparing the quadratic-program matrix. - """ - r = np.array(m, copy=True) - n = len(r) - for k in range(n): - s = 0.0 - if k >= 1: - for i in range(k): - t = r[i, k] - if i > 0: - t = t - np.sum(r[0:i, i] * r[0:i, k]) - t = t / r[i, i] - r[i, k] = t - s = s + t * t - s = r[k, k] - s - if s <= 0.0: - return k + 1, r - r[k, k] = np.sqrt(s) - return 0, r - - def _nearest_pd(self, A): - """Project a matrix onto the cone of positive-definite matrices. - - In some cases the matrix induced by the energy-distance objective is not - numerically positive definite, even though the underlying optimization - problem is well posed. In those cases we replace it with the nearest - positive-definite approximation before calling ``quadprog``. - """ - B = (A + A.T) / 2 - _, s, V = np.linalg.svd(B) - H = V.T @ np.diag(s) @ V - A2 = (B + H) / 2 - A3 = (A2 + A2.T) / 2 - - if self._is_pd(A3): - return A3 - - spacing = np.spacing(np.linalg.norm(A)) - identity_matrix = np.eye(A.shape[0]) - k = 1 - while not self._is_pd(A3): - mineig = np.min(np.real(np.linalg.eigvals(A3))) - A3 += identity_matrix * (-mineig * k ** 2 + spacing) - k += 1 - - return A3 - - def _compute_ed_param_train(self, distance_func, train_distrib, classes, n_cls_i): - """Pre-compute the training-side terms of the ED optimization problem. - - Given the class-conditional posterior clouds observed on validation - data, this routine computes the pairwise average distances between - classes and derives the matrices required by the quadratic program. - These terms depend only on the validation distribution and can therefore - be cached after ``aggregation_fit``. - """ - n_classes = len(classes) - K = np.zeros((n_classes, n_classes), dtype=float) - for i in range(n_classes): - K[i, i] = distance_func(train_distrib[classes[i]], train_distrib[classes[i]]).sum() - for j in range(i + 1, n_classes): - K[i, j] = distance_func(train_distrib[classes[i]], train_distrib[classes[j]]).sum() - K[j, i] = K[i, j] - - K = K / np.dot(n_cls_i, n_cls_i.T) - - B = np.zeros((n_classes - 1, n_classes - 1), dtype=float) - for i in range(n_classes - 1): - B[i, i] = -K[i, i] - K[-1, -1] + 2 * K[i, -1] - for j in range(n_classes - 1): - if j == i: - continue - B[i, j] = -K[i, j] - K[-1, -1] + K[i, -1] + K[j, -1] - - G = 2 * B - if not self._is_pd(G): - G = self._nearest_pd(G) - - C = -np.vstack([np.ones((1, n_classes - 1)), -np.eye(n_classes - 1)]).T - b = -np.array([1] + [0] * (n_classes - 1), dtype=float) - - return K, G, C, b - - def _compute_ed_param_test(self, distance_func, train_distrib, test_distrib, K, classes, n_cls_i): - """Compute the test-dependent linear term of the ED objective. - - Once the training-side matrices have been computed, each new test sample - only requires estimating the distances between its posterior cloud and - the class-conditional validation clouds. - """ - n_classes = len(classes) - Kt = np.zeros(n_classes, dtype=float) - for i in range(n_classes): - Kt[i] = distance_func(train_distrib[classes[i]], test_distrib).sum() - - Kt = Kt / (n_cls_i.squeeze() * float(len(test_distrib))) - return 2 * (-Kt[:-1] + K[:-1, -1] + Kt[-1] - K[-1, -1]) - - def _solve_ed(self, G, a, C, b): - """Solve the energy-distance quadratic program. - - The optimization is carried out over the first ``n_classes - 1`` - prevalences; the prevalence of the last class is recovered afterwards by - the simplex constraint. The resulting vector is finally normalized as a - precaution against small numerical deviations. - """ - quadprog = _get_quadprog() - sol = quadprog.solve_qp(G=G, a=a, C=C, b=b) - prevalences = sol[0] - prevalences = np.append(prevalences, 1 - prevalences.sum()) - return F.normalize_prevalence(prevalences, method='clip') - - def aggregation_fit(self, classif_predictions, labels): - """ - Estimate the class-conditional posterior distributions on validation - data and pre-compute the quadratic-program parameters that depend only - on the training side. - - In EDy, the validation posteriors are not discretized into histograms. - Instead, each class is represented by the cloud of posterior vectors - observed for that class, and these clouds are then compared through the - selected pairwise distance. - - :param classif_predictions: posterior probabilities returned by the - classifier on validation data - :param labels: true labels associated to each posterior vector - """ - posteriors = np.asarray(classif_predictions, dtype=float) - labels = np.asarray(labels) - - self.train_distrib_ = { - class_: posteriors[labels == class_] for class_ in self.classes_ - } - self.train_n_cls_i_ = np.asarray( - [[len(self.train_distrib_[class_])] for class_ in self.classes_], - dtype=float, - ) - - self.K_, self.G_, self.C_, self.b_ = self._compute_ed_param_train( - self.distance, - self.train_distrib_, - self.classes_, - self.train_n_cls_i_, - ) - return self - - def aggregate(self, posteriors: np.ndarray): - """Estimate the prevalence vector for a test sample. - - :param posteriors: posterior probabilities returned by the classifier - for the instances in the test sample - :return: a prevalence vector of shape ``(n_classes,)`` - """ - posteriors = np.asarray(posteriors, dtype=float) - self.a_ = self._compute_ed_param_test( - self.distance, - self.train_distrib_, - posteriors, - self.K_, - self.classes_, - self.train_n_cls_i_, - ) - return self._solve_ed(G=self.G_, a=self.a_, C=self.C_, b=self.b_) diff --git a/quapy/method/_energy.py b/quapy/method/_energy.py new file mode 100644 index 0000000..e0ca8ed --- /dev/null +++ b/quapy/method/_energy.py @@ -0,0 +1,144 @@ +from typing import Callable, Union + +import numpy as np +from sklearn.metrics.pairwise import euclidean_distances, manhattan_distances + +import quapy as qp +import quapy.functional as F +from quapy.method._helper import _get_quadprog + + +class _EnergyDistanceCore: + """Shared numerical core for energy-distance quantifiers.""" + + def _check_ed_init_parameters(self): + self.distance = self._resolve_distance_function(self.distance) + + @staticmethod + def _resolve_distance_function(distance): + if isinstance(distance, str): + if distance == 'manhattan': + return manhattan_distances + if distance == 'euclidean': + return euclidean_distances + raise ValueError( + f"unknown distance {distance!r}; valid aliases are 'manhattan' and 'euclidean'" + ) + if not hasattr(distance, '__call__'): + raise ValueError('distance must be a valid string alias or a callable function') + return distance + + def _is_pd(self, m): + """Check whether a symmetric matrix is positive definite.""" + return self._dpofa(m)[0] == 0 + + def _dpofa(self, m): + """Factor a symmetric positive definite matrix.""" + r = np.array(m, copy=True) + n = len(r) + for k in range(n): + s = 0.0 + if k >= 1: + for i in range(k): + t = r[i, k] + if i > 0: + t = t - np.sum(r[0:i, i] * r[0:i, k]) + t = t / r[i, i] + r[i, k] = t + s = s + t * t + s = r[k, k] - s + if s <= 0.0: + return k + 1, r + r[k, k] = np.sqrt(s) + return 0, r + + def _nearest_pd(self, A): + """Project a matrix onto the cone of positive-definite matrices.""" + B = (A + A.T) / 2 + _, s, V = np.linalg.svd(B) + H = V.T @ np.diag(s) @ V + A2 = (B + H) / 2 + A3 = (A2 + A2.T) / 2 + + if self._is_pd(A3): + return A3 + + spacing = np.spacing(np.linalg.norm(A)) + identity_matrix = np.eye(A.shape[0]) + k = 1 + while not self._is_pd(A3): + mineig = np.min(np.real(np.linalg.eigvals(A3))) + A3 += identity_matrix * (-mineig * k ** 2 + spacing) + k += 1 + + return A3 + + def _compute_ed_param_train(self, distance_func, train_distrib, n_cls_i): + """Pre-compute the training-side terms of the ED optimization problem.""" + n_classes = len(train_distrib) + K = np.zeros((n_classes, n_classes), dtype=float) + for i in range(n_classes): + K[i, i] = distance_func(train_distrib[i], train_distrib[i]).sum() + for j in range(i + 1, n_classes): + K[i, j] = distance_func(train_distrib[i], train_distrib[j]).sum() + K[j, i] = K[i, j] + + K = K / np.dot(n_cls_i, n_cls_i.T) + + B = np.zeros((n_classes - 1, n_classes - 1), dtype=float) + for i in range(n_classes - 1): + B[i, i] = -K[i, i] - K[-1, -1] + 2 * K[i, -1] + for j in range(n_classes - 1): + if j == i: + continue + B[i, j] = -K[i, j] - K[-1, -1] + K[i, -1] + K[j, -1] + + G = 2 * B + if not self._is_pd(G): + G = self._nearest_pd(G) + + C = -np.vstack([np.ones((1, n_classes - 1)), -np.eye(n_classes - 1)]).T + b = -np.array([1] + [0] * (n_classes - 1), dtype=float) + + return K, G, C, b + + def _compute_ed_param_test(self, distance_func, train_distrib, test_distrib, K, n_cls_i): + """Compute the test-dependent linear term of the ED objective.""" + n_classes = len(train_distrib) + Kt = np.zeros(n_classes, dtype=float) + for i in range(n_classes): + Kt[i] = distance_func(train_distrib[i], test_distrib).sum() + + Kt = Kt / (n_cls_i.squeeze() * float(test_distrib.shape[0])) + return 2 * (-Kt[:-1] + K[:-1, -1] + Kt[-1] - K[-1, -1]) + + def _solve_ed(self, G, a, C, b): + """Solve the energy-distance quadratic program.""" + quadprog = _get_quadprog() + sol = quadprog.solve_qp(G=G, a=a, C=C, b=b) + prevalences = sol[0] + prevalences = np.append(prevalences, 1 - prevalences.sum()) + return F.normalize_prevalence(prevalences, method='clip') + + def _fit_energy_model(self, train_distrib): + self.train_distrib_ = tuple(train_distrib) + self.train_n_cls_i_ = np.asarray( + [[distrib.shape[0]] for distrib in self.train_distrib_], + dtype=float, + ) + self.K_, self.G_, self.C_, self.b_ = self._compute_ed_param_train( + self.distance, + self.train_distrib_, + self.train_n_cls_i_, + ) + return self + + def _predict_energy(self, test_distrib): + self.a_ = self._compute_ed_param_test( + self.distance, + self.train_distrib_, + test_distrib, + self.K_, + self.train_n_cls_i_, + ) + return self._solve_ed(G=self.G_, a=self.a_, C=self.C_, b=self.b_) diff --git a/quapy/method/aggregative.py b/quapy/method/aggregative.py index 2aa160c..f678546 100644 --- a/quapy/method/aggregative.py +++ b/quapy/method/aggregative.py @@ -16,6 +16,7 @@ from quapy.functional import get_divergence from quapy.classification.svmperf import SVMperf from quapy.data import LabelledCollection from quapy.method.base import BaseQuantifier, BinaryQuantifier, OneVsAllGeneric +from quapy.method._energy import _EnergyDistanceCore from quapy.method._helper import ( _get_abstention_calibrators, _get_cvxpy, @@ -1712,6 +1713,106 @@ class AggregativeMedianEstimator(BinaryQuantifier): return np.median(prev_preds, axis=0) +class EDy(_EnergyDistanceCore, AggregativeSoftQuantifier): + """ + Energy Distance y (EDy), a posterior-space distribution-matching quantifier + based on energy distance. + + The method represents each class by the posterior-probability vectors + produced by a probabilistic classifier on validation data, and estimates the + test prevalence vector by matching the test posterior distribution against + the class-conditional validation distributions through an energy-distance + objective solved as a quadratic program. The method is therefore another + instance of the general mixture-matching view of quantification, but it + operates directly on posterior vectors rather than on histogram summaries. + + This implementation works for binary and multiclass single-label + quantification and relies on the optional ``quadprog`` dependency. It was + adapted to QuaPy's current aggregative API from the original implementation + available in `quantificationlib `_, + and now shares its numerical core with the classifier-free + :class:`quapy.method.non_aggregative.EDx` variant. + + The current implementation follows the energy-distance formulation discussed + in: + + * Alberto Castaño, Laura Morán-Fernández, Jaime Alonso, + Verónica Bolón-Canedo, Amparo Alonso-Betanzos, and Juan José del Coz. + *An analysis of quantification methods based on matching distributions*. + * Hideko Kawakubo, Marthinus Christoffel du Plessis, and Masashi Sugiyama + (2016). *Computationally efficient class-prior estimation under class + balance change using energy distance*. IEICE Transactions on Information + and Systems, 99(1):176-186. + + :param classifier: a scikit-learn ``BaseEstimator``, or ``None`` to use + ``qp.environ['DEFAULT_CLS']`` + :param fit_classifier: whether to train the learner (default ``True``). + Set to ``False`` if the learner has already been trained outside the + quantifier + :param val_split: specification of the data used for generating validation + posterior probabilities. This can be an integer (default ``5``) for + k-fold cross-validation, a float in ``(0, 1)`` for a held-out split, + or a tuple ``(X, y)`` with explicit validation data + :param distance: distance used to compare posterior vectors. Valid string + aliases are ``'manhattan'`` (default) and ``'euclidean'``; a custom + callable compatible with pairwise-distance signatures can also be used + :param n_jobs: number of parallel workers (default ``None``, meaning the + value is taken from the environment) + """ + + def __init__( + self, + classifier: BaseEstimator = None, + fit_classifier: bool = True, + val_split=5, + distance: Union[str, Callable] = 'manhattan', + n_jobs=None, + ): + super().__init__(classifier, fit_classifier, val_split) + self.distance = distance + self.n_jobs = qp._get_njobs(n_jobs) + self.train_n_cls_i_ = None + self.train_distrib_ = None + self.K_ = None + self.G_ = None + self.C_ = None + self.b_ = None + self.a_ = None + + def _check_init_parameters(self): + self._check_ed_init_parameters() + + def aggregation_fit(self, classif_predictions, labels): + """ + Estimate the class-conditional posterior distributions on validation + data and pre-compute the quadratic-program parameters that depend only + on the training side. + + In EDy, the validation posteriors are not discretized into histograms. + Instead, each class is represented by the cloud of posterior vectors + observed for that class, and these clouds are then compared through the + selected pairwise distance. + + :param classif_predictions: posterior probabilities returned by the + classifier on validation data + :param labels: true labels associated to each posterior vector + """ + posteriors = np.asarray(classif_predictions, dtype=float) + labels = np.asarray(labels) + train_distrib = [posteriors[labels == class_] for class_ in self.classes_] + return self._fit_energy_model(train_distrib) + + def aggregate(self, posteriors: np.ndarray): + """Estimate the prevalence vector for a test sample. + + :param posteriors: posterior probabilities returned by the classifier + for the instances in the test sample + :return: a prevalence vector of shape ``(n_classes,)`` + """ + posteriors = np.asarray(posteriors, dtype=float) + return self._predict_energy(posteriors) + + # --------------------------------------------------------------- # imports # --------------------------------------------------------------- @@ -1730,9 +1831,6 @@ KDEyML = _kdey.KDEyML KDEyHD = _kdey.KDEyHD KDEyCS = _kdey.KDEyCS -from . import _edy - -EDy = _edy.EDy # --------------------------------------------------------------- # aliases diff --git a/quapy/method/non_aggregative.py b/quapy/method/non_aggregative.py index e587f05..0250852 100644 --- a/quapy/method/non_aggregative.py +++ b/quapy/method/non_aggregative.py @@ -10,9 +10,11 @@ from quapy.method.confidence import WithConfidenceABC, ConfidenceRegionABC from quapy.functional import get_divergence from quapy.method.base import BaseQuantifier, BinaryQuantifier from quapy.method._helper import _labels_to_indices +from quapy.method._energy import _EnergyDistanceCore import quapy.functional as F from scipy.optimize import lsq_linear from scipy import sparse +import quapy as qp class MaximumLikelihoodPrevalenceEstimation(BaseQuantifier): @@ -160,6 +162,70 @@ class DMx(BaseQuantifier): return F.argmin_prevalence(loss, n_classes, method=self.search) +class EDx(_EnergyDistanceCore, BaseQuantifier): + """ + Energy Distance x (EDx), a covariate-space distribution-matching + quantifier based on energy distance. + + EDx is the classifier-free counterpart of :class:`quapy.method.aggregative.EDy`. + Instead of representing each class through posterior-probability vectors, it + represents each class by the cloud of raw feature vectors observed in the + training set and estimates the test prevalence vector by solving the same + energy-distance quadratic program directly in feature space. + + This implementation works for binary and multiclass single-label + quantification and relies on the optional ``quadprog`` dependency. The + current QuaPy adaptation shares its numerical core with EDy and keeps + credit to the original implementation available in + `quantificationlib `_. + + The formulation follows the same references as EDy, namely: + + * Alberto Castaño, Laura Morán-Fernández, Jaime Alonso, + Verónica Bolón-Canedo, Amparo Alonso-Betanzos, and Juan José del Coz. + *An analysis of quantification methods based on matching distributions*. + * Hideko Kawakubo, Marthinus Christoffel du Plessis, and Masashi Sugiyama + (2016). *Computationally efficient class-prior estimation under class + balance change using energy distance*. IEICE Transactions on Information + and Systems, 99(1):176-186. + + :param distance: distance used to compare feature vectors. Valid string + aliases are ``'manhattan'`` (default) and ``'euclidean'``; a custom + callable compatible with pairwise-distance signatures can also be used + :param n_jobs: number of parallel workers (default ``None``, meaning the + value is taken from the environment) + """ + + def __init__(self, distance: Union[str, Callable] = 'manhattan', n_jobs=None): + self.distance = distance + self.n_jobs = qp._get_njobs(n_jobs) + self.classes_ = None + self.n_features_in_ = None + self.train_distrib_ = None + self.train_n_cls_i_ = None + self.K_ = None + self.G_ = None + self.C_ = None + self.b_ = None + self.a_ = None + + def fit(self, X, y): + """Fit class-conditional feature-space distributions from training data.""" + self._check_ed_init_parameters() + labels = np.asarray(y) + self.classes_ = np.unique(labels) + self.n_features_in_ = X.shape[1] + train_distrib = [X[labels == class_] for class_ in self.classes_] + return self._fit_energy_model(train_distrib) + + def predict(self, X): + """Estimate class prevalences for a test sample of raw instances.""" + assert X.shape[1] == self.n_features_in_, ( + f'wrong shape; expected {self.n_features_in_}, found {X.shape[1]}' + ) + return self._predict_energy(X) + + class ReadMe(BaseQuantifier, WithConfidenceABC): """ ReadMe is a non-aggregative quantification system proposed by @@ -336,6 +402,8 @@ def _get_features_range(X): # aliases #--------------------------------------------------------------- + HDx = DMx.HDx DistributionMatchingX = DMx +EnergyDistanceX = EDx HellingerDistanceX = HDx \ No newline at end of file diff --git a/quapy/tests/test_methods.py b/quapy/tests/test_methods.py index bb9b060..e059a4d 100644 --- a/quapy/tests/test_methods.py +++ b/quapy/tests/test_methods.py @@ -7,7 +7,7 @@ import numpy as np from sklearn.linear_model import LogisticRegression from quapy.method import AGGREGATIVE_METHODS, BINARY_METHODS, NON_AGGREGATIVE_METHODS -from quapy.method.non_aggregative import DMx, HDx +from quapy.method.non_aggregative import DMx, EDx, HDx from quapy.method.aggregative import ACC, DMy, EDy, KDEyCS, RLLS from quapy.method.meta import Ensemble from quapy.functional import check_prevalence_vector @@ -23,6 +23,10 @@ OPTIONAL_AGGREGATIVE_METHODS = { 'EDy', } +OPTIONAL_NON_AGGREGATIVE_METHODS = { + 'EDx', +} + class TestMethods(unittest.TestCase): @@ -56,6 +60,8 @@ class TestMethods(unittest.TestCase): def test_non_aggregative(self): for dataset in TestMethods.datasets: for model in NON_AGGREGATIVE_METHODS: + if model.__name__ in OPTIONAL_NON_AGGREGATIVE_METHODS: + continue if not dataset.binary and model in BINARY_METHODS: continue @@ -153,6 +159,19 @@ class TestMethods(unittest.TestCase): self.assertTrue(check_prevalence_vector(estim_prevalences)) + def test_edx(self): + try: + import quadprog # noqa: F401 + except ImportError: + return + + dataset = TestMethods.tiny_dataset_multiclass + q = EDx() + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + + def test_dmy_noncanonical_labels(self): dataset = TestMethods.tiny_dataset_multiclass label_names = np.asarray(['class-a', 'class-c', 'class-z'])