merged from devel

.github/workflows/ci.yml

@@ -6,6 +6,8 @@ on:
     branches:
       - master
       - devel
+    tags:
+      - "[0-9]+.[0-9]+.[0-9]+"
 
 jobs:
 
@@ -28,7 +30,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip setuptools wheel
-        python -m pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.4"
+        python -m pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@main"
         python -m pip install -e .[bayes,tests]
     - name: Test with unittest
       run: python -m unittest
@@ -47,7 +49,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip setuptools wheel "jax[cpu]"
-        python -m pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.4"
+        python -m pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@main"
         python -m pip install -e .[neural,docs]
     - name: Build documentation
       run: sphinx-build -M html docs/source docs/build
@@ -66,3 +68,37 @@ jobs:
         branch: gh-pages
         directory: __gh-pages/
         github_token: ${{ secrets.GITHUB_TOKEN }}
+
+  release:
+    name: Build & Publish Release
+    runs-on: ubuntu-latest
+    if: startsWith(github.ref, 'refs/tags/')
+    steps:
+      - uses: actions/checkout@v4
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+      - name: Install build dependencies
+        run: |
+          python -m pip install --upgrade pip build twine
+      - name: Build package
+        run: python -m build
+      - name: Publish to PyPI
+        uses: pypa/gh-action-pypi-publish@release/v1
+        with:
+          user: __token__
+          password: ${{ secrets.PYPI_API_TOKEN }}
+      - name: Create GitHub Release
+        id: create_release
+        uses: actions/create-release@v1
+        with:
+          tag_name: ${{ github.ref_name }}
+          release_name: Release ${{ github.ref_name }}
+          body: |
+            Changes in this release: 
+            - see commit history for details
+          draft: false
+          prerelease: false
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

CHANGE_LOG.txt

@@ -1,3 +1,38 @@
+Change Log 0.2.0
+-----------------
+
+CLEAN TODO-FILE
+
+- Base code Refactor:
+    - Removing coupling between LabelledCollection and quantification methods; the fit interface changes:
+        def fit(data:LabelledCollection): -> def fit(X, y):
+    - Adding function "predict" (function "quantify" is still present as an alias, for the nostalgic)
+    - Aggregative methods's behavior in terms of fit_classifier and how to treat the val_split is now
+        indicated exclusively at construction time, and it is no longer possible to indicate it at fit time.
+        This is because, in v<=0.1.9, one could create a method (e.g., ACC) and then indicate:
+        my_acc.fit(tr_data, fit_classifier=False, val_split=val_data)
+        in which case the first argument is unused, and this was ambiguous with
+        my_acc.fit(the_data, fit_classifier=False)
+        in which case the_data is to be used for validation purposes. However, the val_split could be set as a fraction
+        indicating only part of the_data must be used for validation, and the rest wasted... it was certainly confusing.
+    - This change imposes a versioning constrain with qunfold, which now must be >= 0.1.6
+- EMQ has been modified, so that the representation function "classify" now only provides posterior
+    probabilities and, if required, these are recalibrated (e.g., by "bcts") during the aggregation function.
+    - A new parameter "on_calib_error" is passed to the constructor, which informs of the policy to follow
+        in case the abstention's calibration functions failed (which happens sometimes). Options include:
+            - 'raise': raises a RuntimeException (default)
+            - 'backup': reruns by silently avoiding calibration
+    - Parameter "recalib" has been renamed "calib"
+- Added aggregative bootstrap for deriving confidence regions (confidence intervals, ellipses in the simplex, or
+    ellipses in the CLR space). This method is efficient as it leverages the two-phases of the aggregative quantifiers.
+    This method applies resampling only to the aggregation phase, thus avoiding to train many quantifiers, or
+    classify multiple times the instances of a sample. See:
+    - quapy/method/confidence.py (new)
+    - the new example no. 16.confidence_regions.py
+- BayesianCC moved to confidence.py, where methods having to do with confidence intervals belong.
+- Improved documentation of qp.plot module.
+
+
 Change Log 0.1.9
 ----------------
 

README.md

@@ -13,8 +13,8 @@ for facilitating the analysis and interpretation of the experimental results.
 
 ### Last updates:
 
-* Version 0.1.9 is released! major changes can be consulted [here](CHANGE_LOG.txt).
+* Version 0.2.0 is released! major changes can be consulted [here](CHANGE_LOG.txt).
-* The developer API documentation is available [here](https://hlt-isti.github.io/QuaPy/index.html)
+* The developer API documentation is available [here](https://hlt-isti.github.io/QuaPy/build/html/modules.html)
 
 ### Installation
 
@@ -46,15 +46,15 @@ of the test set.
 ```python
 import quapy as qp
 
-dataset = qp.datasets.fetch_UCIBinaryDataset("yeast")
+training, test = qp.datasets.fetch_UCIBinaryDataset("yeast").train_test
-training, test = dataset.train_test
 
 # create an "Adjusted Classify & Count" quantifier
 model = qp.method.aggregative.ACC()
-model.fit(training)
+Xtr, ytr = training.Xy
+model.fit(Xtr, ytr)
 
-estim_prevalence = model.quantify(test.X)
+estim_prevalence = model.predict(test.X)
-true_prevalence  = test.prevalence()
+true_prevalence = test.prevalence()
 
 error = qp.error.mae(true_prevalence, estim_prevalence)
 print(f'Mean Absolute Error (MAE)={error:.3f}')
@@ -67,8 +67,7 @@ class prevalence of the training set. For this reason, any quantification model
 should be tested across many samples, even ones characterized by class prevalence 
 values different or very different from those found in the training set.
 QuaPy implements sampling procedures and evaluation protocols that automate this workflow.
-See the [documentation](https://hlt-isti.github.io/QuaPy/manuals/protocols.html) 
+See the [documentation](https://hlt-isti.github.io/QuaPy/build/html/) for detailed examples.
-and the [examples directory](https://github.com/HLT-ISTI/QuaPy/tree/master/examples) for detailed examples.
 
 ## Features
 
@@ -80,8 +79,8 @@ quantification methods based on structured output learning, HDy, QuaNet, quantif
     * 32 UCI Machine Learning datasets.
     * 11 Twitter quantification-by-sentiment datasets.
     * 3 product reviews quantification-by-sentiment datasets. 
-    * 4 tasks from LeQua 2022 competition
+    * 4 tasks from LeQua 2022 competition and 4 tasks from LeQua 2024 competition
-    * 4 tasks from LeQua 2024 competition (_new in v0.1.9!_)
+    * IFCB for Plancton quantification 
 * Native support for binary and single-label multiclass quantification scenarios.
 * Model selection functionality that minimizes quantification-oriented loss functions.
 * Visualization tools for analysing the experimental results.
@@ -102,22 +101,23 @@ In case you want to contribute improvements to quapy, please generate pull reque
   
 ## Documentation
 
-The developer API documentation is available [here](https://hlt-isti.github.io/QuaPy/). 
+The [developer API documentation](https://hlt-isti.github.io/QuaPy/build/html/modules.html) is available [here](https://hlt-isti.github.io/QuaPy/build/html/index.html). 
 
-Check out our [Manuals](https://hlt-isti.github.io/QuaPy/manuals.html), in which many examples
+Check out the [Manuals](https://hlt-isti.github.io/QuaPy/manuals.html), in which many code examples
 are provided:
 
 * [Datasets](https://hlt-isti.github.io/QuaPy/manuals/datasets.html)
 * [Evaluation](https://hlt-isti.github.io/QuaPy/manuals/evaluation.html)
-* [Explicit loss minimization](https://hlt-isti.github.io/QuaPy/manuals/explicit-loss-minimization.html)
-* [Methods](https://hlt-isti.github.io/QuaPy/manuals/methods.html)
-* [Model Selection](https://hlt-isti.github.io/QuaPy/manuals/datasets.html)
-* [Plotting](https://hlt-isti.github.io/QuaPy/manuals/plotting.html)
 * [Protocols](https://hlt-isti.github.io/QuaPy/manuals/protocols.html)
+* [Methods](https://hlt-isti.github.io/QuaPy/manuals/methods.html)
+* [SVMperf](https://hlt-isti.github.io/QuaPy/manuals/explicit-loss-minimization.html)
+* [Model Selection](https://hlt-isti.github.io/QuaPy/manuals/model-selection.html)
+* [Plotting](https://hlt-isti.github.io/QuaPy/manuals/plotting.html)
 
 ## Acknowledgments:
 
-This work has been funded by the QuaDaSh project (P2022TB5JF) "Finanziato dall’Unione europea- Next Generation EU, Missione 4 Componente 2 CUP B53D23026250001".
-
-<img src="docs/source/EUfooter.png" alt="EUcommission" width="1000"/>
 <img src="docs/source/SoBigData.png" alt="SoBigData++" width="250"/>
+
+This work has been supported by the QuaDaSh project 
+_"Finanziato dall’Unione europea---Next Generation EU, 
+Missione 4 Componente 2 CUP B53D23026250001"_.

TODO.txt

@@ -1,6 +1,60 @@
+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
+based on the number of actual positives and the required sample size; e.g., for sample_size=100 and actual
+positives [10, 100, 500] -> [0.1, 1.0, 1.0], meaning that class 0 can be sampled at most at 0.1 prevalence, while
+the others can be sampled up to 1. prevalence. Then, when a prevalence value is requested, e.g., [0.33, 0.33, 0.33],
+we may either clip each value and normalize (as you suggest for the extreme case, e.g., [0.1, 0.33, 0.33]/sum) or
+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.
+
+Para quitar el labelledcollection de los métodos:
+
+- El follón viene por la semántica confusa de fit en agregativos, que recibe 3 parámetros:
+    - data: LabelledCollection, que puede ser:
+        - el training set si hay que entrenar el clasificador
+        - None si no hay que entregar el clasificador
+        - el validation, que entra en conflicto con val_split, si no hay que entrenar clasificador
+    - fit_classifier: dice si hay que entrenar el clasificador o no, y estos cambia la semántica de los otros
+    - val_split: que puede ser:
+        - un número: el número de kfcv, lo cual implica fit_classifier=True y data=todo el training set
+        - una fración en [0,1]: que indica la parte que usamos para validation; implica fit_classifier=True y data=train+val
+        - un labelled collection: el conjunto de validación específico; no implica fit_classifier=True ni False
+- La forma de quitar la dependencia de los métodos con LabelledCollection debería ser así:
+    - En el constructor se dice si el clasificador que se recibe por parámetro hay que entrenarlo o ya está entrenado;
+        es decir, hay un fit_classifier=True o False.
+        - fit_classifier=True:
+            - data en fit es todo el training incluyendo el validation y todo
+            - val_split:
+                - int: número de folds en kfcv
+                - proporción en [0,1]
+        - fit_classifier=False:
+
+
+
+- [TODO] document confidence in manuals
+- [TODO] Test the return_type="index" in protocols and finish the "distributing_samples.py" example
+- [TODO] Add EDy (an implementation is available at quantificationlib)
 - [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

docs/source/manuals/datasets.md

@@ -340,10 +340,10 @@ and a set of test samples (for evaluation). QuaPy returns this data as a Labelle
 (training) and two generation protocols (for validation and test samples), as follows:
 
 ```python
-training, val_generator, test_generator = fetch_lequa2022(task=task)
+training, val_generator, test_generator = qp.datasets.fetch_lequa2022(task=task)
 ```
 
-See the `lequa2022_experiments.py` in the examples folder for further details on how to
+See the `5a.lequa2022_experiments.py` in the examples folder for further details on how to
 carry out experiments using these datasets.  
 
 The datasets are downloaded only once, and stored for fast reuse.
@@ -365,6 +365,53 @@ Esuli, A., Moreo, A., Sebastiani, F., & Sperduti, G. (2022).
 A Detailed Overview of LeQua@ CLEF 2022: Learning to Quantify.
 ```
 
+## LeQua 2024 Datasets
+
+QuaPy also provides the datasets used for the [LeQua 2024 competition](https://lequa2024.github.io/). 
+In brief, there are 4 tasks:
+* T1: binary quantification (by sentiment)
+* T2: multiclass quantification (28 classes, merchandise products)
+* T3: ordinal quantification (5-stars sentiment ratings)
+* T4: binary sentiment quantification under a combination of covariate shift and prior shift
+
+In all cases, the covariate space has 256 dimensions (extracted using the `ELECTRA-Small` model).
+
+Every task consists of a training set, a set of validation samples (for model selection)
+and a set of test samples (for evaluation). QuaPy returns this data as a LabelledCollection
+(training bags) and sampling generation protocols (for validation and test bags). 
+T3 also offers the possibility to obtain a series of training bags (in form of a 
+sampling generation protocol) instead of one single training bag. Use it as follows:
+
+```python
+training, val_generator, test_generator = qp.datasets.fetch_lequa2024(task=task)
+```
+
+See the `5b.lequa2024_experiments.py` in the examples folder for further details on how to
+carry out experiments using these datasets.  
+
+The datasets are downloaded only once, and stored for fast reuse.
+
+Some statistics are summarized below:
+
+| Dataset | classes | train size  | validation samples | test samples | docs by sample |   type   |
+|---------|:-------:|:-----------:|:------------------:|:------------:|:--------------:|:--------:| 
+| T1      |    2    |    5000     |        1000        |     5000     |      250       |  vector  | 
+| T2      |   28    |    20000    |        1000        |     5000     |      1000      |  vector  |
+| T3      |    5    | 100 samples |        1000        |     5000     |      200       |   vector   |
+| T4      |    2    |    5000     |        1000        |     5000     |      250       |   vector   |
+
+For further details on the datasets or the competition, we refer to 
+[the official site](https://lequa2024.github.io/data/) and
+[the overview paper](http://nmis.isti.cnr.it/sebastiani/Publications/LQ2024.pdf).
+
+```
+Esuli, A., Moreo, A., Sebastiani, F., & Sperduti, G. (2022).
+An Overview of LeQua 2024, the 2nd International Data Challenge on Learning to Quantify,
+Proceedings of the 4th International Workshop on Learning to Quantify (LQ 2024), 
+ECML-PKDD 2024, Vilnius, Lithuania.
+```
+
+
 ## IFCB Plankton dataset
 
 IFCB is a dataset of plankton species in water samples hosted in `Zenodo <https://zenodo.org/records/10036244>`_.
@@ -402,12 +449,20 @@ train, test_gen = qp.datasets.fetch_IFCB(for_model_selection=False, single_sampl
 # ... train and evaluation
 ```
 
+See also [Automatic plankton quantification using deep features
+P González, A Castaño, EE Peacock, J Díez, JJ Del Coz, HM Sosik
+Journal of Plankton Research 41 (4), 449-463](https://par.nsf.gov/servlets/purl/10172325).
+
 
 
 ## Adding Custom Datasets
 
+It is straightforward to import your own datasets into QuaPy. 
+I what follows, there are some code snippets for doing so; see also the example
+[3.custom_collection.py](https://github.com/HLT-ISTI/QuaPy/blob/master/examples/3.custom_collection.py).
+
 QuaPy provides data loaders for simple formats dealing with 
-text, following the format:
+text; for example, use `qp.data.reader.from_text` for the following the format:
 
 ```
 class-id \t first document's pre-processed text \n
@@ -415,13 +470,16 @@ class-id \t second document's pre-processed text \n
 ...
 ```
 
-and sparse representations of the form:
+or `qp.data.reader.from_sparse` for sparse representations of the form:
 
 ```
 {-1, 0, or +1} col(int):val(float) col(int):val(float) ... \n
 ...
 ```
 
+both functions return a tuple `X, y` containing a list of strings and the corresponding 
+labels, respectively.
+
 The code in charge in loading a LabelledCollection is:
 
 ```python
@@ -430,12 +488,13 @@ def load(cls, path:str, loader_func:callable):
     return LabelledCollection(*loader_func(path))
 ```
 
-indicating that any _loader_func_ (e.g., a user-defined one) which 
+indicating that any `loader_func` (e.g., `from_text`, `from_sparse`, `from_csv`, or a user-defined one) which 
 returns valid arguments for initializing a _LabelledCollection_ object will allow
-to load any collection. In particular, the _LabelledCollection_ receives as 
+to load any collection. More specifically, the _LabelledCollection_ receives as 
-arguments the instances (as an iterable) and the labels (as an iterable) and,
+arguments the _instances_ (iterable) and the _labels_ (iterable) and,
-additionally, the number of classes can be specified (it would otherwise be
+optionally, the number of classes (it would be
-inferred from the labels, but that requires at least one positive example for
+inferred from the labels if not indicated, but this requires at least one 
+positive example for
 all classes to be present in the collection).
 
 The same _loader_func_ can be passed to a Dataset, along with two 
@@ -448,20 +507,23 @@ import quapy as qp
 train_path = '../my_data/train.dat'
 test_path = '../my_data/test.dat'
 
-def my_custom_loader(path):
+def my_custom_loader(path, **custom_kwargs):
     with open(path, 'rb') as fin:
         ...
     return instances, labels
 
-data = qp.data.Dataset.load(train_path, test_path, my_custom_loader)
+data = qp.data.Dataset.load(train_path, test_path, my_custom_loader, **custom_kwargs)
 ```
 
 ### Data Processing
 
-QuaPy implements a number of preprocessing functions in the package _qp.data.preprocessing_, including:
+QuaPy implements a number of preprocessing functions in the package `qp.data.preprocessing`, including:
 
 * _text2tfidf_: tfidf vectorization 
 * _reduce_columns_: reducing the number of columns based on term frequency
 * _standardize_: transforms the column values into z-scores (i.e., subtract the mean and normalizes by the standard deviation, so
 that the column values have zero mean and unit variance).
-* _index_: transforms textual tokens into lists of numeric ids) 
+* _index_: transforms textual tokens into lists of numeric ids
+
+These functions are applied to `Dataset` objects, and offer the possibility to apply the transformation
+inline (thus modifying the original dataset), or to return a modified copy.

docs/source/manuals/evaluation.md

@@ -46,18 +46,18 @@ e.g.:
 
 ```python
 qp.environ['SAMPLE_SIZE'] = 100  # once for all
-true_prev = np.asarray([0.5, 0.3, 0.2])  # let's assume 3 classes
+true_prev = [0.5, 0.3, 0.2]  # let's assume 3 classes
-estim_prev = np.asarray([0.1, 0.3, 0.6])
+estim_prev = [0.1, 0.3, 0.6]
 error = qp.error.mrae(true_prev, estim_prev)
 print(f'mrae({true_prev}, {estim_prev}) = {error:.3f}')
 ```
 
 will print:
 ```
-mrae([0.500, 0.300, 0.200], [0.100, 0.300, 0.600]) = 0.914
+mrae([0.5, 0.3, 0.2], [0.1, 0.3, 0.6]) = 0.914
 ```
 
-Finally, it is possible to instantiate QuaPy's quantification
+It is also possible to instantiate QuaPy's quantification
 error functions from strings using, e.g.:
 
 ```python
@@ -85,7 +85,7 @@ print(f'MAE = {mae:.4f}')
 ```
 
 It is often desirable to evaluate our system using more than one
-single evaluatio measure. In this case, it is convenient to generate
+single evaluation measure. In this case, it is convenient to generate
 a _report_. A report in QuaPy is a dataframe accounting for all the
 true prevalence values with their corresponding prevalence values
 as estimated by the quantifier, along with the error each has given
@@ -104,7 +104,7 @@ report['estim-prev'] = report['estim-prev'].map(F.strprev)
 print(report)
 
 print('Averaged values:')
-print(report.mean())
+print(report.mean(numeric_only=True))
 ```
 
 This will produce an output like:
@@ -141,11 +141,14 @@ true_prevs, estim_prevs = qp.evaluation.prediction(quantifier, protocol=prot)
 
 All the evaluation functions implement specific optimizations for speeding-up 
 the evaluation of aggregative quantifiers (i.e., of instances of _AggregativeQuantifier_).
+
 The optimization comes down to generating classification predictions (either crisp or soft) 
 only once for the entire test set, and then applying the sampling procedure to the
 predictions, instead of generating samples of instances and then computing the 
 classification predictions every time. This is only possible when the protocol
-is an instance of _OnLabelledCollectionProtocol_. The optimization is only 
+is an instance of _OnLabelledCollectionProtocol_. 
+
+The optimization is only 
 carried out when the number of classification predictions thus generated would be
 smaller than the number of predictions required for the entire protocol; e.g., 
 if the original dataset contains 1M instances, but the protocol is such that it would
@@ -156,4 +159,4 @@ precompute all the predictions irrespectively of the number of instances and num
 Finally, this can be deactivated by setting _aggr_speedup=False_. Note that this optimization
 is not only applied for the final evaluation, but also for the internal evaluations carried
 out during _model selection_. Since these are typically many, the heuristic can help reduce the
-execution time a lot.
+execution time significatively.

docs/source/manuals/methods.md

@@ -1,7 +1,7 @@
 # Quantification Methods
 
 Quantification methods can be categorized as belonging to
-`aggregative` and `non-aggregative` groups. 
+`aggregative`, `non-aggregative`, and `meta-learning` groups. 
 Most methods included in QuaPy at the moment are of type `aggregative`
 (though we plan to add many more methods in the near future), i.e.,
 are methods characterized by the fact that
@@ -12,21 +12,17 @@ Any quantifier in QuaPy shoud extend the class `BaseQuantifier`,
 and implement some abstract methods:
 ```python
     @abstractmethod
-    def fit(self, data: LabelledCollection): ...
+    def fit(self, X, y): ...
 
     @abstractmethod
-    def quantify(self, instances): ...
+    def predict(self, X): ...
 ```
 The meaning of those functions should be familiar to those
 used to work with scikit-learn since the class structure of QuaPy
 is directly inspired by scikit-learn's _Estimators_. Functions
-`fit` and `quantify` are used to train the model and to provide
+`fit` and `predict` (for which there is an alias `quantify`) 
-class estimations (the reason why
+are used to train the model and to provide
-scikit-learn' structure has not been adopted _as is_ in QuaPy responds to 
+class estimations.
-the fact that scikit-learn's `predict` function is expected to return
-one output for each input element --e.g., a predicted label for each
-instance in a sample-- while in quantification the output for a sample
-is one single array of class prevalences).
 Quantifiers also extend from scikit-learn's `BaseEstimator`, in order
 to simplify the use of `set_params` and `get_params` used in 
 [model selection](./model-selection).
@@ -40,21 +36,26 @@ The methods that any `aggregative` quantifier must implement are:
 
 ```python
     @abstractmethod
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
 
     @abstractmethod
-    def aggregate(self, classif_predictions:np.ndarray): ...
+    def aggregate(self, classif_predictions): ...
 ```
 
-These two functions replace the `fit` and `quantify` methods, since those
+The argument `classif_predictions` is whatever the method `classify` returns. 
-come with default implementations. The `fit` function is provided and amounts to: 
+QuaPy comes with default implementations that cover most common cases, but you can
+override `classify` in case your method requires further or different information to work.
+
+These two functions replace the `fit` and `predict` methods, which
+come with default implementations. For instance, the `fit` function is 
+provided and amounts to: 
 
 ```python
-def fit(self, data: LabelledCollection, fit_classifier=True, val_split=None):
+    def fit(self, X, y):
-    self._check_init_parameters()
+        self._check_init_parameters()
-    classif_predictions = self.classifier_fit_predict(data, fit_classifier, predict_on=val_split)
+        classif_predictions, labels = self.classifier_fit_predict(X, y)
-    self.aggregation_fit(classif_predictions, data)
+        self.aggregation_fit(classif_predictions, labels)
-    return self
+        return self
 ```
 
 Note that this function fits the classifier, and generates the predictions. This is assumed
@@ -72,11 +73,11 @@ overriden (if needed) and allows the method to quickly raise any exception based
 found in the `__init__` arguments, thus avoiding to break after training the classifier and generating
 predictions.
 
-Similarly, the function `quantify` is provided, and amounts to:
+Similarly, the function `predict` (alias `quantify`) is provided, and amounts to:
 
 ```python
-def quantify(self, instances):
+def predict(self, X):
-    classif_predictions = self.classify(instances)
+    classif_predictions = self.classify(X)
     return self.aggregate(classif_predictions)
 ```
 
@@ -84,12 +85,14 @@ in which only the function `aggregate` is required to be overriden in most cases
 
 Aggregative quantifiers are expected to maintain a classifier (which is
 accessed through the `@property` `classifier`). This classifier is
-given as input to the quantifier, and can be already fit
+given as input to the quantifier, and will be trained by the quantifier's fit (default).
-on external data (in which case, the `fit_learner` argument should
+Alternatively, the classifier can be already fit on external data; in this case, the `fit_learner` 
-be set to False), or be fit by the quantifier's fit (default).
+argument in the `__init__` should be set to False (see [4.using_pretrained_classifier.py](https://github.com/HLT-ISTI/QuaPy/blob/master/examples/4.using_pretrained_classifier.py)
+for a full code example).
 
-The above patterns (in training: fit the classifier, then fit the aggregation; 
+The above patterns (in training: (i) fit the classifier, then (ii) fit the aggregation; 
-in test: classify, then aggregate) allows QuaPy to optimize many internal procedures.
+in test: (i) classify, then (ii) aggregate) allows QuaPy to optimize many internal procedures, 
+on the grounds that steps (i) are slower than steps (ii). 
 In particular, the model selection routing takes advantage of this two-step process
 and generates classifiers only for the valid combinations of hyperparameters of the 
 classifier, and then _clones_ these classifiers and explores the combinations
@@ -124,6 +127,7 @@ import quapy.functional as F
 from sklearn.svm import LinearSVC
 
 training, test = qp.datasets.fetch_twitter('hcr', pickle=True).train_test
+Xtr, ytr = training.Xy
 
 # instantiate a classifier learner, in this case a SVM
 svm = LinearSVC()
@@ -131,8 +135,8 @@ svm = LinearSVC()
 # instantiate a Classify & Count with the SVM
 # (an alias is available in qp.method.aggregative.ClassifyAndCount)
 model = qp.method.aggregative.CC(svm)
-model.fit(training)
+model.fit(Xtr, ytr)
-estim_prevalence = model.quantify(test.instances)
+estim_prevalence = model.predict(test.instances)
 ```
 
 The same code could be used to instantiate an ACC, by simply replacing
@@ -153,26 +157,14 @@ predictions. This parameters can also be set with an integer,
 indicating that the parameters should be estimated by means of
 _k_-fold cross-validation, for which the integer indicates the
 number _k_ of folds (the default value is 5). Finally, `val_split` can be set to a 
-specific held-out validation set (i.e., an instance of `LabelledCollection`).
+specific held-out validation set (i.e., an tuple `(X,y)`).
-
-The specification of `val_split` can be
-postponed to the invokation of the fit method (if `val_split` was also
-set in the constructor, the one specified at fit time would prevail), 
-e.g.:
-
-```python
-model = qp.method.aggregative.ACC(svm)
-# perform 5-fold cross validation for estimating ACC's parameters
-# (overrides the default val_split=0.4 in the constructor)
-model.fit(training, val_split=5)
-```
 
 The following code illustrates the case in which PCC is used:
 
 ```python
 model = qp.method.aggregative.PCC(svm)
-model.fit(training)
+model.fit(Xtr, ytr)
-estim_prevalence = model.quantify(test.instances)
+estim_prevalence = model.predict(Xte)
 print('classifier:', model.classifier)
 ```
 In this case, QuaPy will print:
@@ -185,11 +177,11 @@ is not a probabilistic classifier (i.e., it does not implement the
 `predict_proba` method) and so, the classifier will be converted to
 a probabilistic one through [calibration](https://scikit-learn.org/stable/modules/calibration.html).
 As a result, the classifier that is printed in the second line points
-to a `CalibratedClassifier` instance. Note that calibration can only
+to a `CalibratedClassifierCV` instance. Note that calibration can only
-be applied to hard classifiers when `fit_learner=True`; an exception 
+be applied to hard classifiers if `fit_learner=True`; an exception 
 will be raised otherwise.
 
-Lastly, everything we said aboud ACC and PCC
+Lastly, everything we said about ACC and PCC
 applies to PACC as well.
 
 _New in v0.1.9_: quantifiers ACC and PACC now have three additional arguments: `method`, `solver` and `norm`:
@@ -221,7 +213,7 @@ Options are:
   * `"condsoftmax"`  applies softmax normalization only if the prevalence vector lies outside of the probability simplex.
 
 
-#### BayesianCC (_New in v0.1.9_!)
+#### BayesianCC
 
 The `BayesianCC` is a variant of ACC introduced in 
 [Ziegler, A. and Czyż, P. "Bayesian quantification with black-box estimators", arXiv (2023)](https://arxiv.org/abs/2302.09159), 
@@ -259,29 +251,35 @@ An example of use can be found below:
 import quapy as qp
 from sklearn.linear_model import LogisticRegression
 
-dataset = qp.datasets.fetch_twitter('hcr', pickle=True)
+train, test = qp.datasets.fetch_twitter('hcr', pickle=True).train_test
 
 model = qp.method.aggregative.EMQ(LogisticRegression())
-model.fit(dataset.training)
+model.fit(*train.Xy)
-estim_prevalence = model.quantify(dataset.test.instances)
+estim_prevalence = model.predict(test.X)
 ```
 
-_New in v0.1.7_: EMQ now accepts two new parameters in the construction method, namely
+EMQ accepts additional parameters in the construction method:
-`exact_train_prev` which allows to use the true training prevalence as the departing
+* `exact_train_prev`: set to True for using the true training prevalence as the departing
-prevalence estimation (default behaviour), or instead an approximation of it as 
+prevalence estimation (default behaviour), or to False for using an approximation of it as
 suggested by [Alexandari et al. (2020)](http://proceedings.mlr.press/v119/alexandari20a.html) 
-(by setting `exact_train_prev=False`).
+* `calib`: allows to indicate a calibration method, among those
-The other parameter is `recalib` which allows to indicate a calibration method, among those
 proposed by [Alexandari et al. (2020)](http://proceedings.mlr.press/v119/alexandari20a.html),
-including the Bias-Corrected Temperature Scaling, Vector Scaling, etc.
+including the Bias-Corrected Temperature Scaling 
-See the API documentation for further details. 
+(`bcts`), Vector Scaling (`bcts`), No-Bias Temperature Scaling (`nbvs`), 
+or Temperature Scaling (`ts`); default is `None` (no calibration).
+* `on_calib_error`: indicates the policy to follow in case the calibrator fails at runtime.
+        Options include `raise` (default), in which case a RuntimeException is raised; and `backup`, in which
+        case the calibrator is silently skipped.
+
+You can use the class method `EMQ_BCTS` to effortlessly instantiate EMQ with the best performing
+heuristics found by [Alexandari et al. (2020)](http://proceedings.mlr.press/v119/alexandari20a.html). See the API documentation for further details. 
 
 
 ### Hellinger Distance y (HDy)
 
 Implementation of the method based on the Hellinger Distance y (HDy) proposed by
-[González-Castro, V., Alaiz-Rodrı́guez, R., and Alegre, E. (2013). Class distribution
+[González-Castro, V., Alaiz-Rodríguez, R., and Alegre, E. (2013). Class distribution
-estimation based on the Hellinger distance. Information Sciences, 218:146–164.](https://www.sciencedirect.com/science/article/pii/S0020025512004069)
+estimation based on the Hellinger distance. Information Sciences, 218:146-164.](https://www.sciencedirect.com/science/article/pii/S0020025512004069)
 
 It is implemented in `qp.method.aggregative.HDy` (also accessible
 through the allias `qp.method.aggregative.HellingerDistanceY`).
@@ -289,16 +287,16 @@ This method works with a probabilistic classifier (hard classifiers
 can be used as well and will be calibrated) and requires a validation
 set to estimate parameter for the mixture model. Just like 
 ACC and PACC, this quantifier receives a `val_split` argument
-in the constructor (or in the fit method, in which case the previous
+in the constructor that can either be a float indicating the proportion
-value is overridden) that can either be a float indicating the proportion
 of training data to be taken as the validation set (in a random
-stratified split), or a validation set (i.e., an instance of 
+stratified split), or the validation set itself (i.e., an tuple
-`LabelledCollection`) itself. 
+`(X,y)`). 
 
 HDy was proposed as a binary classifier and the implementation
 provided in QuaPy accepts only binary datasets. 
  
 The following code shows an example of use:
+
 ```python
 import quapy as qp
 from sklearn.linear_model import LogisticRegression
@@ -308,11 +306,11 @@ dataset = qp.datasets.fetch_reviews('hp', pickle=True)
 qp.data.preprocessing.text2tfidf(dataset, min_df=5, inplace=True)
 
 model = qp.method.aggregative.HDy(LogisticRegression())
-model.fit(dataset.training)
+model.fit(*dataset.training.Xy)
-estim_prevalence = model.quantify(dataset.test.instances)
+estim_prevalence = model.predict(dataset.test.X)
 ```
 
-_New in v0.1.7:_ QuaPy now provides an implementation of the generalized
+QuaPy also provides an implementation of the generalized
 "Distribution Matching" approaches for multiclass, inspired by the framework
 of [Firat (2016)](https://arxiv.org/abs/1606.00868). One can instantiate
 a variant of HDy for multiclass quantification as follows:
@@ -321,17 +319,22 @@ a variant of HDy for multiclass quantification as follows:
 mutliclassHDy = qp.method.aggregative.DMy(classifier=LogisticRegression(), divergence='HD', cdf=False)
 ``` 
 
-_New in v0.1.7:_ QuaPy now provides an implementation of the "DyS"
+QuaPy also provides an implementation of the "DyS"
 framework proposed by [Maletzke et al (2020)](https://ojs.aaai.org/index.php/AAAI/article/view/4376)
 and the "SMM" method proposed by [Hassan et al (2019)](https://ieeexplore.ieee.org/document/9260028)
 (thanks to _Pablo González_ for the contributions!)
 
 ### Threshold Optimization methods
 
-_New in v0.1.7:_ QuaPy now implements Forman's threshold optimization methods;
+QuaPy implements Forman's threshold optimization methods;
 see, e.g., [(Forman 2006)](https://dl.acm.org/doi/abs/10.1145/1150402.1150423) 
 and [(Forman 2008)](https://link.springer.com/article/10.1007/s10618-008-0097-y).
-These include: T50, MAX, X, Median Sweep (MS), and its variant MS2.
+These include: `T50`, `MAX`, `X`, Median Sweep (`MS`), and its variant `MS2`.
+
+These methods are binary-only and implement different heuristics for 
+improving the stability of the denominator of the ACC adjustment (`tpr-fpr`).
+The methods are called "threshold" since said heuristics have to do
+with different choices of the underlying classifier's threshold.
 
 ### Explicit Loss Minimization
 
@@ -411,19 +414,21 @@ qp.environ['SVMPERF_HOME'] = '../svm_perf_quantification'
 
 model = newOneVsAll(SVMQ(), n_jobs=-1)  # run them on parallel
 model.fit(dataset.training)
-estim_prevalence = model.quantify(dataset.test.instances)
+estim_prevalence = model.predict(dataset.test.instances)
 ```
 
-Check the examples on [explicit_loss_minimization](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/5.explicit_loss_minimization.py)
+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. 
 
 ### Kernel Density Estimation methods (KDEy)
 
-_New in v0.1.8_: QuaPy now provides implementations for the three variants
+QuaPy provides implementations for the three variants
 of KDE-based methods proposed in 
-_[Moreo, A., González, P. and del Coz, J.J., 2023. 
+_[Moreo, A., González, P. and del Coz, J.J.. 
 Kernel Density Estimation for Multiclass Quantification. 
-arXiv preprint arXiv:2401.00490.](https://arxiv.org/abs/2401.00490)_. 
+Machine Learning. Vol 114 (92), 2025](https://link.springer.com/article/10.1007/s10994-024-06726-5)_
+(a [preprint](https://arxiv.org/abs/2401.00490) is available online). 
 The variants differ in the divergence metric to be minimized:
 
 - KDEy-HD: minimizes the (squared) Hellinger Distance and solves the problem via a Monte Carlo approach
@@ -434,30 +439,42 @@ These methods are specifically devised for multiclass problems (although they ca
 binary problems too). 
 
 All KDE-based methods depend on the hyperparameter `bandwidth` of the kernel. Typical values
-that can be explored in model selection range in [0.01, 0.25]. The methods' performance
+that can be explored in model selection range in [0.01, 0.25]. Previous experiments reveal the methods' performance
-vary smoothing with smooth variations of this hyperparameter.
+varies smoothly at small variations of this hyperparameter.
 
 
 ## Composable Methods
 
-The [](quapy.method.composable) module allows the composition of quantification methods from loss functions and feature transformations. Any composed method solves a linear system of equations by minimizing the loss after transforming the data. Methods of this kind include ACC, PACC, HDx, HDy, and many other well-known methods, as well as an unlimited number of re-combinations of their building blocks.
+The `quapy.method.composable` module integrates [qunfold](https://github.com/mirkobunse/qunfold) allows the composition
+of quantification methods from loss functions and feature transformations (thanks to Mirko Bunse for the integration!). 
+
+Any composed method solves a linear system of equations by minimizing the loss after transforming the data. Methods of this kind include ACC, PACC, HDx, HDy, and many other well-known methods, as well as an unlimited number of re-combinations of their building blocks.
 
 ### Installation
 
 ```sh
 pip install --upgrade pip setuptools wheel
 pip install "jax[cpu]"
-pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.4"
+pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.5"
 ```
 
+**Note:** since version 0.2.0, QuaPy is only compatible with qunfold >=0.1.5. 
+
 ### Basics
 
 The composition of a method is implemented through the [](quapy.method.composable.ComposableQuantifier) class. Its documentation also features an example to get you started in composing your own methods.
 
 ```python
+from quapy.method.composable import (
+    ComposableQuantifier,
+    TikhonovRegularized,
+    LeastSquaresLoss,
+    ClassRepresentation,
+)
+
 ComposableQuantifier( # ordinal ACC, as proposed by Bunse et al., 2022
-  TikhonovRegularized(LeastSquaresLoss(), 0.01),
+    TikhonovRegularized(LeastSquaresLoss(), 0.01),
-  ClassTransformer(RandomForestClassifier(oob_score=True))
+    ClassRepresentation(RandomForestClassifier(oob_score=True))
 )
 ```
 
@@ -484,16 +501,16 @@ You can use the [](quapy.method.composable.CombinedLoss) to create arbitrary, we
 
 ### Feature transformations
 
-- [](quapy.method.composable.ClassTransformer)
+- [](quapy.method.composable.ClassRepresentation)
-- [](quapy.method.composable.DistanceTransformer)
+- [](quapy.method.composable.DistanceRepresentation)
-- [](quapy.method.composable.HistogramTransformer)
+- [](quapy.method.composable.HistogramRepresentation)
-- [](quapy.method.composable.EnergyKernelTransformer)
+- [](quapy.method.composable.EnergyKernelRepresentation)
-- [](quapy.method.composable.GaussianKernelTransformer)
+- [](quapy.method.composable.GaussianKernelRepresentation)
-- [](quapy.method.composable.LaplacianKernelTransformer)
+- [](quapy.method.composable.LaplacianKernelRepresentation)
-- [](quapy.method.composable.GaussianRFFKernelTransformer)
+- [](quapy.method.composable.GaussianRFFKernelRepresentation)
 
 ```{hint}
-The [](quapy.method.composable.ClassTransformer) requires the classifier to have a property `oob_score==True` and to produce a property `oob_decision_function` during fitting. In [scikit-learn](https://scikit-learn.org/), this requirement is fulfilled by any bagging classifier, such as random forests. Any other classifier needs to be cross-validated through the [](quapy.method.composable.CVClassifier).
+The [](quapy.method.composable.ClassRepresentation) requires the classifier to have a property `oob_score==True` and to produce a property `oob_decision_function` during fitting. In [scikit-learn](https://scikit-learn.org/), this requirement is fulfilled by any bagging classifier, such as random forests. Any other classifier needs to be cross-validated through the [](quapy.method.composable.CVClassifier).
 ```
 
 
@@ -528,10 +545,11 @@ from quapy.method.meta import Ensemble
 from sklearn.linear_model import LogisticRegression
 
 dataset = qp.datasets.fetch_UCIBinaryDataset('haberman')
+train, test = dataset.train_test
 
 model = Ensemble(quantifier=ACC(LogisticRegression()), size=30, policy='ave', n_jobs=-1)
-model.fit(dataset.training)
+model.fit(*train.Xy)
-estim_prevalence = model.quantify(dataset.test.instances)
+estim_prevalence = model.predict(test.X)
 ```
 
 Other aggregation policies implemented in QuaPy include:
@@ -578,7 +596,29 @@ learner = NeuralClassifierTrainer(cnn, device='cuda')
 
 # train QuaNet
 model = QuaNet(learner, device='cuda')
-model.fit(dataset.training)
+model.fit(*dataset.training.Xy)
-estim_prevalence = model.quantify(dataset.test.instances)
+estim_prevalence = model.predict(dataset.test.X)
 ```
 
+## Confidence Regions for Class Prevalence Estimation
+
+_(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 methods are currently implemented:
+
+* Aggregative Bootstrap: The Aggregative Bootstrap method extends any aggregative quantifier by generating confidence regions for class prevalence estimates through bootstrapping. Key features of this method include:
+
+    * Optimized Computation: The bootstrap is applied to pre-classified instances, significantly speeding up training and inference.
+During training, bootstrap repetitions are performed only after training the classifier once. These repetitions are used to train multiple aggregation functions.
+During inference, bootstrap is applied over pre-classified test instances.
+  * General Applicability: Aggregative Bootstrap can be applied to any aggregative quantifier.
+  For further information, check the [example](https://github.com/HLT-ISTI/QuaPy/tree/master/examples/16.confidence_regions.py) provided.
+
+* BayesianCC: is a Bayesian variant of the Adjusted Classify & Count (ACC) quantifier; see more details in the [example](https://github.com/HLT-ISTI/QuaPy/tree/master/examples/14.bayesian_quantification.py) provided.
+
+Confidence regions are constructed around a point estimate, which is typically computed as the mean value of a set of samples.
+The confidence region can be instantiated in three ways:
+* Confidence intervals: are standard confidence intervals generated for each class independently (_method="intervals"_).
+* Confidence ellipse in the simplex: an ellipse constructed around the mean point; the ellipse lies on the simplex and takes
+  into account possible inter-class dependencies in the data (_method="ellipse"_). 
+* Confidence ellipse in the Centered-Log Ratio (CLR) space: the underlying assumption of the ellipse is that the components are
+  normally distributed. However, we know elements from the simplex have an inner structure. A better approach is to first
+  transform the components into an unconstrained space (the CLR), and then construct the ellipse in such space (_method="ellipse-clr"_).

docs/source/manuals/model-selection.md

@@ -87,7 +87,7 @@ model = qp.model_selection.GridSearchQ(
     error='mae',  # the error to optimize is the MAE (a quantification-oriented loss)
     refit=True,  # retrain on the whole labelled set once done
     verbose=True  # show information as the process goes on
-).fit(training)
+).fit(*training.Xy)
 
 print(f'model selection ended: best hyper-parameters={model.best_params_}')
 model = model.best_model_
@@ -133,7 +133,7 @@ learner = GridSearchCV(
     LogisticRegression(),
     param_grid={'C': np.logspace(-4, 5, 10), 'class_weight': ['balanced', None]},
     cv=5)
-model = DistributionMatching(learner).fit(dataset.train)
+model = DistributionMatching(learner).fit(*dataset.train.Xy)
 ```
 
 However, this is conceptually flawed, since the model should be

docs/source/manuals/plotting.md

@@ -2,6 +2,9 @@
 
 The module _qp.plot_ implements some basic plotting functions
 that can help analyse the performance of a quantification method.
+See the provided 
+[code example](https://github.com/HLT-ISTI/QuaPy/blob/master/examples/13.plotting.py) 
+for a full example. 
 
 All plotting functions receive as inputs the outcomes of 
 some experiments and include, for each experiment, 
@@ -77,7 +80,7 @@ def gen_data():
     method_names, true_prevs, estim_prevs, tr_prevs = [], [], [], []
 
     for method_name, model in models():
-        model.fit(train)
+        model.fit(*train.Xy)
         true_prev, estim_prev = qp.evaluation.prediction(model, APP(test, repeats=100, random_state=0))
 
         method_names.append(method_name)
@@ -171,7 +174,7 @@ def gen_data():
         training_size = 5000
         # since the problem is binary, it suffices to specify the negative prevalence, since the positive is constrained
         train_sample = train.sampling(training_size, 1-training_prevalence)
-        model.fit(train_sample)
+        model.fit(*train_sample.Xy)
         true_prev, estim_prev = qp.evaluation.prediction(model, APP(test, repeats=100, random_state=0))
         method_name = 'CC$_{'+f'{int(100*training_prevalence)}' + '\%}$'
         method_data.append((method_name, true_prev, estim_prev, train_sample.prevalence()))

docs/source/manuals/protocols.md

@@ -1,7 +1,5 @@
 # Protocols
 
-_New in v0.1.7!_
-
 Quantification methods are expected to behave robustly in the presence of 
 shift. For this reason, quantification methods need to be confronted with
 samples exhibiting widely varying amounts of shift. 
@@ -106,15 +104,16 @@ train, test = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=5).train_test
 
 # model selection
 train, val = train.split_stratified(train_prop=0.75)
+Xtr, ytr = train.Xy
 quantifier = qp.model_selection.GridSearchQ(
     quantifier, 
     param_grid={'classifier__C': np.logspace(-2, 2, 5)}, 
     protocol=APP(val)  # <- this is the protocol we use for generating validation samples
-).fit(train)
+).fit(Xtr, ytr)
 
 # default values are n_prevalences=21, repeats=10, random_state=0; this is equialent to:
 # val_app = APP(val, n_prevalences=21, repeats=10, random_state=0)
-# quantifier = GridSearchQ(quantifier, param_grid, protocol=val_app).fit(train)
+# quantifier = GridSearchQ(quantifier, param_grid, protocol=val_app).fit(Xtr, ytr)
 
 # evaluation with APP
 mae = qp.evaluation.evaluate(quantifier, protocol=APP(test), error_metric='mae')

examples/0.basics.py

@@ -6,6 +6,7 @@ import numpy as np
 from sklearn.linear_model import LogisticRegression
 
 import quapy as qp
+from quapy.method.aggregative import PACC
 
 # let's fetch some dataset to run one experiment
 # datasets are available in the "qp.data.datasets" module (there is a shortcut in qp.datasets)
@@ -33,17 +34,15 @@ import quapy.functional as F  # <- this module has some functional utilities, li
 print(f'training prevalence = {F.strprev(train.prevalence())}')
 
 # let us train one quantifier, for example, PACC using a sklearn's Logistic Regressor as the underlying classifier
-# classifier = LogisticRegression()
+classifier = LogisticRegression()
-
+pacc = PACC(classifier)
-# pacc = qp.method.aggregative.PACC(classifier)
-pacc = qp.method.aggregative.PACC()
 
 print(f'training {pacc}')
-pacc.fit(train)
+pacc.fit(X, y)
 
 # let's now test our quantifier on the test data (of course, we should not use the test labels y at this point, only X)
 X_test = test.X
-estim_prevalence = pacc.quantify(X_test)
+estim_prevalence = pacc.predict(X_test)
 
 print(f'estimated test prevalence = {F.strprev(estim_prevalence)}')
 print(f'true test prevalence = {F.strprev(test.prevalence())}')

examples/1.model_selection.py

@@ -12,15 +12,24 @@ In this example, we show how to perform model selection on a DistributionMatchin
 model = DMy()
 
 qp.environ['SAMPLE_SIZE'] = 100
+qp.environ['N_JOBS'] = -1
 
 print(f'running model selection with N_JOBS={qp.environ["N_JOBS"]}; '
-      f'to increase the number of jobs use:\n> N_JOBS=-1 python3 1.model_selection.py\n'
+      f'to increase/decrease the number of jobs use:\n'
+      f'> N_JOBS=-1 python3 1.model_selection.py\n'
       f'alternatively, you can set this variable within the script as:\n'
       f'import quapy as qp\n'
       f'qp.environ["N_JOBS"]=-1')
 
 training, test = qp.datasets.fetch_UCIMulticlassDataset('letter').train_test
 
+# evaluation in terms of MAE with default hyperparameters
+Xtr, ytr = training.Xy
+model.fit(Xtr, ytr)
+mae_score = qp.evaluation.evaluate(model, protocol=UPP(test), error_metric='mae')
+print(f'MAE (non optimized)={mae_score:.5f}')
+
+
 with qp.util.temp_seed(0):
 
     # The model will be returned by the fit method of GridSearchQ.
@@ -50,6 +59,7 @@ with qp.util.temp_seed(0):
 
     tinit = time()
 
+    Xtr, ytr = training.Xy
     model = qp.model_selection.GridSearchQ(
         model=model,
         param_grid=param_grid,
@@ -58,7 +68,7 @@ with qp.util.temp_seed(0):
         refit=False,   # retrain on the whole labelled set once done
         # raise_errors=False,
         verbose=True  # show information as the process goes on
-    ).fit(training)
+    ).fit(Xtr, ytr)
 
 tend = time()
 

examples/10.one_vs_all.py

@@ -9,6 +9,11 @@ import numpy as np
 """
 In this example, we will create a quantifier for tweet sentiment analysis considering three classes: negative, neutral,
 and positive. We will use a one-vs-all approach using a binary quantifier for demonstration purposes.
+
+Caveat: the one-vs-all approach is deemed inadequate under prior probability shift conditions. The reasons
+are discussed in:
+Donyavi, Z., Serapio, A., & Batista, G. (2023). MC-SQ: A highly accurate ensemble for multi-class quantifi-
+cation. In: Proceedings of the 2023 SIAM International Conference on Data Mining (SDM), SIAM, pp. 622–630
 """
 
 qp.environ['SAMPLE_SIZE'] = 100
@@ -40,11 +45,11 @@ param_grid = {
 }
 print('starting model selection')
 model_selection = GridSearchQ(quantifier, param_grid, protocol=UPP(val), verbose=True, refit=False)
-quantifier = model_selection.fit(train_modsel).best_model()
+quantifier = model_selection.fit(*train_modsel.Xy).best_model()
 
 print('training on the whole training set')
 train, test = qp.datasets.fetch_twitter('hcr', for_model_selection=False, pickle=True).train_test
-quantifier.fit(train)
+quantifier.fit(*train.Xy)
 
 # evaluation
 mae = qp.evaluation.evaluate(quantifier, protocol=UPP(test), error_metric='mae')

examples/11.comparing_HDy_HDx.py

@@ -23,8 +23,9 @@ qp.environ['SAMPLE_SIZE']=100
 
 df = pd.DataFrame(columns=['method', 'dataset', 'MAE', 'MRAE', 'tr-time', 'te-time'])
 
+datasets = qp.datasets.UCI_BINARY_DATASETS
 
-for dataset_name in tqdm(qp.datasets.UCI_BINARY_DATASETS, total=len(qp.datasets.UCI_BINARY_DATASETS)):
+for dataset_name in tqdm(datasets, total=len(datasets), desc='datasets processed'):
     if dataset_name in ['acute.a', 'acute.b', 'balance.2', 'iris.1']:
         # these datasets tend to produce either too good or too bad results...
         continue
@@ -32,23 +33,25 @@ for dataset_name in tqdm(qp.datasets.UCI_BINARY_DATASETS, total=len(qp.datasets.
     collection = qp.datasets.fetch_UCIBinaryLabelledCollection(dataset_name, verbose=False)
     train, test = collection.split_stratified()
 
+    Xtr, ytr = train.Xy
+
     # HDy............................................
     tinit = time()
-    hdy = HDy(LogisticRegression()).fit(train)
+    hdy = HDy(LogisticRegression()).fit(Xtr, ytr)
     t_hdy_train = time()-tinit
 
     tinit = time()
-    hdy_report = qp.evaluation.evaluation_report(hdy, APP(test), error_metrics=['mae', 'mrae']).mean()
+    hdy_report = qp.evaluation.evaluation_report(hdy, APP(test), error_metrics=['mae', 'mrae']).mean(numeric_only=True)
     t_hdy_test = time() - tinit
     df.loc[len(df)] = ['HDy', dataset_name, hdy_report['mae'], hdy_report['mrae'], t_hdy_train, t_hdy_test]
 
     # HDx............................................
     tinit = time()
-    hdx = DMx.HDx(n_jobs=-1).fit(train)
+    hdx = DMx.HDx(n_jobs=-1).fit(Xtr, ytr)
     t_hdx_train = time() - tinit
 
     tinit = time()
-    hdx_report = qp.evaluation.evaluation_report(hdx, APP(test), error_metrics=['mae', 'mrae']).mean()
+    hdx_report = qp.evaluation.evaluation_report(hdx, APP(test), error_metrics=['mae', 'mrae']).mean(numeric_only=True)
     t_hdx_test = time() - tinit
     df.loc[len(df)] = ['HDx', dataset_name, hdx_report['mae'], hdx_report['mrae'], t_hdx_train, t_hdx_test]
 

examples/12.custom_protocol.py

@@ -3,14 +3,13 @@ from sklearn.linear_model import LogisticRegression
 
 import quapy as qp
 from quapy.method.aggregative import PACC
-from quapy.data import LabelledCollection
 from quapy.protocol import AbstractStochasticSeededProtocol
 import quapy.functional as F
 
 """
 In this example, we create a custom protocol.
-The protocol generates samples of a Gaussian mixture model with random mixture parameter (the sample prevalence).
+The protocol generates synthetic samples of a Gaussian mixture model with random mixture parameter 
-Datapoints are univariate and we consider 2 classes only.
+(the sample prevalence). Datapoints are univariate and we consider 2 classes only for simplicity.
 """
 class GaussianMixProtocol(AbstractStochasticSeededProtocol):
     # We need to extend AbstractStochasticSeededProtocol if we want the samples to be replicable
@@ -81,10 +80,9 @@ with qp.util.temp_seed(0):
     Xpos = np.random.normal(loc=mu_2, scale=std_2, size=100)
     X = np.concatenate([Xneg, Xpos]).reshape(-1,1)
     y = [0]*100 + [1]*100
-    training = LabelledCollection(X, y)
 
     pacc = PACC(LogisticRegression())
-    pacc.fit(training)
+    pacc.fit(X, y)
 
 
 mae = qp.evaluation.evaluate(pacc, protocol=gm, error_metric='mae', verbose=True)

examples/13.plotting.py

@@ -0,0 +1,73 @@
+import quapy as qp
+import numpy as np
+
+from protocol import APP
+from quapy.method.aggregative import CC, ACC, PCC, PACC
+from sklearn.svm import LinearSVC
+
+qp.environ['SAMPLE_SIZE'] = 500
+
+
+'''
+In this example, we show how to create some plots for the analysis of experimental results.
+The main functions are included in qp.plot but, before, we will generate some basic experimental data
+'''
+
+def gen_data():
+    # this function generates some experimental data to plot
+
+    def base_classifier():
+        return LinearSVC(class_weight='balanced')
+
+    def datasets():
+        # the plots can handle experiments in different datasets
+        yield qp.datasets.fetch_reviews('kindle', tfidf=True, min_df=5).train_test
+        # by uncommenting thins line, the experiments will be carried out in more than one dataset
+        # yield qp.datasets.fetch_reviews('hp', tfidf=True, min_df=5).train_test
+
+    def models():
+        yield 'CC', CC(base_classifier())
+        yield 'ACC', ACC(base_classifier())
+        yield 'PCC', PCC(base_classifier())
+        yield 'PACC', PACC(base_classifier())
+
+    # these are the main parameters we need to fill for generating the plots;
+    # note that each these list must have the same number of elements, since the ith entry of each list regards
+    # an independent experiment
+    method_names, true_prevs, estim_prevs, tr_prevs = [], [], [], []
+
+    for train, test in datasets():
+        for method_name, model in models():
+            model.fit(*train.Xy)
+            true_prev, estim_prev = qp.evaluation.prediction(model, APP(test, repeats=100, random_state=0))
+
+            # gather all the data for this experiment
+            method_names.append(method_name)
+            true_prevs.append(true_prev)
+            estim_prevs.append(estim_prev)
+            tr_prevs.append(train.prevalence())
+
+    return method_names, true_prevs, estim_prevs, tr_prevs
+
+# generate some experimental data
+method_names, true_prevs, estim_prevs, tr_prevs = gen_data()
+# if you want to play around with the different plots and parameters, you might prefer to generate the data only once,
+# so you better replace the above line of code with this one, that pickles the experimental results for faster reuse
+# method_names, true_prevs, estim_prevs, tr_prevs = qp.util.pickled_resource('./plots/data.pickle', gen_data)
+
+# if there is only one training prevalence, we can display it
+only_train_prev = tr_prevs[0] if len(np.unique(tr_prevs, axis=0))==1 else None
+
+# diagonal plot (useful for analyzing the performance of quantifiers on binary data)
+qp.plot.binary_diagonal(method_names, true_prevs, estim_prevs,
+                        train_prev=only_train_prev, savepath='./plots/bin_diag.png')
+
+# bias plot (box plots displaying the bias of each method)
+qp.plot.binary_bias_global(method_names, true_prevs, estim_prevs, savepath='./plots/bin_bias.png')
+
+# error by drift allows to plot the quantification error as a function of the amount of prior probability shift, and
+# is preferable than diagonal plots for multiclass datasets
+qp.plot.error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
+                       error_name='ae', n_bins=10, savepath='./plots/err_drift.png')
+
+# each functions return (fig, ax) objects from matplotlib; use them to customize the plots to your liking

examples/14.bayesian_quantification.py

@@ -13,7 +13,7 @@ $ pip install quapy[bayesian]
 Running the script via:
 
 ```
-$ python examples/13.bayesian_quantification.py
+$ python examples/14.bayesian_quantification.py
 ```
 
 will produce a plot `bayesian_quantification.pdf`.
@@ -29,7 +29,8 @@ import quapy as qp
 
 from sklearn.ensemble import RandomForestClassifier
 
-from quapy.method.aggregative import BayesianCC, ACC, PACC
+from quapy.method.aggregative import ACC, PACC
+from method.confidence import BayesianCC
 from quapy.data import LabelledCollection, Dataset
 
 
@@ -121,18 +122,18 @@ def get_random_forest() -> RandomForestClassifier:
 def _get_estimate(estimator_class, training: LabelledCollection, test: np.ndarray) -> None:
     """Auxiliary method for running ACC and PACC."""
     estimator = estimator_class(get_random_forest())
-    estimator.fit(training)
+    estimator.fit(*training.Xy)
-    return estimator.quantify(test)
+    return estimator.predict(test)
 
 
 def train_and_plot_bayesian_quantification(ax: plt.Axes, training: LabelledCollection, test: LabelledCollection) -> None:
     """Fits Bayesian quantification and plots posterior mean as well as individual samples"""
     print('training model Bayesian CC...', end='')
     quantifier = BayesianCC(classifier=get_random_forest())
-    quantifier.fit(training)
+    quantifier.fit(*training.Xy)
 
     # Obtain mean prediction
-    mean_prediction = quantifier.quantify(test.X)
+    mean_prediction = quantifier.predict(test.X)
     mae = qp.error.mae(test.prevalence(), mean_prediction)
     x_ax = np.arange(training.n_classes)
     ax.plot(x_ax, mean_prediction, c="salmon", linewidth=2, linestyle=":", label="Bayesian")

examples/15.composable_methods.py

@@ -1,6 +1,6 @@
 """
 This example illustrates the composition of quantification methods from
-arbitrary loss functions and feature transformations. It will extend the basic
+arbitrary loss functions and feature representations. It will extend the basic
 example on the usage of quapy with this composition.
 
 This example requires the installation of qunfold, the back-end of QuaPy's
@@ -8,7 +8,7 @@ composition module:
 
     pip install --upgrade pip setuptools wheel
     pip install "jax[cpu]"
-    pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.4"
+    pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.5"
 """
 
 import numpy as np
@@ -22,22 +22,23 @@ data = qp.data.preprocessing.text2tfidf(
     min_df = 5,
 )
 training, testing = data.train_test
+Xtr, ytr = training.Xy
 
 # We start by recovering PACC from its building blocks, a LeastSquaresLoss and
-# a probabilistic ClassTransformer. A 5-fold cross-validation is implemented
+# a probabilistic ClassRepresentation. A 5-fold cross-validation is implemented
 # through a CVClassifier.
 
 from quapy.method.composable import (
     ComposableQuantifier,
     LeastSquaresLoss,
-    ClassTransformer,
+    ClassRepresentation,
     CVClassifier,
 )
 from sklearn.linear_model import LogisticRegression
 
 pacc = ComposableQuantifier(
     LeastSquaresLoss(),
-    ClassTransformer(
+    ClassRepresentation(
         CVClassifier(LogisticRegression(random_state=0), 5),
         is_probabilistic = True
     ),
@@ -46,7 +47,7 @@ pacc = ComposableQuantifier(
 # Let's evaluate this quantifier.
 
 print(f"Evaluating PACC: {pacc}")
-pacc.fit(training)
+pacc.fit(Xtr, ytr)
 app = qp.protocol.APP(testing, sample_size=100, n_prevalences=21, repeats=1)
 absolute_errors = qp.evaluation.evaluate(
     model = pacc,
@@ -63,14 +64,14 @@ from quapy.method.composable import HellingerSurrogateLoss
 
 model = ComposableQuantifier(
     HellingerSurrogateLoss(), # the loss is different from before
-    ClassTransformer( # we use the same transformer
+    ClassRepresentation( # we use the same representation
         CVClassifier(LogisticRegression(random_state=0), 5),
         is_probabilistic = True
     ),
 )
 
 print(f"Evaluating {model}")
-model.fit(training)
+model.fit(Xtr, ytr)
 absolute_errors = qp.evaluation.evaluate(
     model = model,
     protocol = app, # use the same protocol for evaluation
@@ -79,7 +80,7 @@ absolute_errors = qp.evaluation.evaluate(
 print(f"MAE = {np.mean(absolute_errors):.4f}+-{np.std(absolute_errors):.4f}")
 
 # In general, any composed method solves a linear system of equations by
-# minimizing the loss after transforming the data. Methods of this kind include
+# minimizing the loss after representing the data. Methods of this kind include
 # ACC, PACC, HDx, HDy, and many other well-known methods, as well as an
 # unlimited number of re-combinations of their building blocks.
 
@@ -93,18 +94,18 @@ from quapy.method.composable import CombinedLoss
 
 model = ComposableQuantifier(
     CombinedLoss(HellingerSurrogateLoss(), LeastSquaresLoss()),
-    ClassTransformer(
+    ClassRepresentation(
         CVClassifier(LogisticRegression(random_state=0), 5),
         is_probabilistic = True
     ),
 )
 
-from qunfold.quapy import QuaPyWrapper
+from quapy.method.composable import QUnfoldWrapper
-from qunfold import GenericMethod
+from qunfold import LinearMethod
 
-model = QuaPyWrapper(GenericMethod(
+model = QUnfoldWrapper(LinearMethod(
     CombinedLoss(HellingerSurrogateLoss(), LeastSquaresLoss()),
-    ClassTransformer(
+    ClassRepresentation(
         CVClassifier(LogisticRegression(random_state=0), 5),
         is_probabilistic = True
     ),
@@ -115,7 +116,7 @@ model = QuaPyWrapper(GenericMethod(
 
 param_grid = {
     "loss__weights": [ (w, 1-w) for w in [.1, .5, .9] ],
-    "transformer__classifier__estimator__C": [1e-1, 1e1],
+    "representation__classifier__estimator__C": [1e-1, 1e1],
 }
 
 grid_search = qp.model_selection.GridSearchQ(
@@ -125,7 +126,7 @@ grid_search = qp.model_selection.GridSearchQ(
     error = "mae",
     refit = False,
     verbose = True,
-).fit(training)
+).fit(Xtr, ytr)
 print(
     f"Best hyper-parameters = {grid_search.best_params_}",
     f"Best MAE = {grid_search.best_score_}",

examples/16.KDEy_bandwidth.py

@@ -0,0 +1,83 @@
+import quapy as qp
+import numpy as np
+from quapy.protocol import UPP
+from quapy.method.aggregative import KDEyML
+import quapy.functional as F
+from time import time
+
+"""
+Let see one example:  
+"""
+
+# load some data
+qp.environ['SAMPLE_SIZE'] = 100
+data = qp.datasets.fetch_UCIMulticlassDataset('molecular')
+training, test = data.train_test
+training, validation = training.split_stratified(train_prop=0.7, random_state=0)
+protocol = UPP(validation)
+
+hyper_C = np.logspace(-3, 3, 7)
+
+model = KDEyML()
+
+with qp.util.temp_seed(0):
+
+    param_grid = {
+        'classifier__C': hyper_C,
+        'bandwidth': np.linspace(0.01, 0.20, 20) # [0.01, 0.02, 0.03, ..., 0.20]
+    }
+
+    model = qp.model_selection.GridSearchQ(
+        model=model,
+        param_grid=param_grid,
+        protocol=protocol,
+        error='mae',  # the error to optimize is the MAE (a quantification-oriented loss)
+        refit=False,  # retrain on the whole labelled set once done
+        n_jobs=-1,
+        verbose=True  # show information as the process goes on
+    ).fit(training)
+
+best_params = model.best_params_
+took = model.fit_time_
+model = model.best_model_
+print(f'model selection ended: best hyper-parameters={best_params}')
+
+# evaluation in terms of MAE
+# we use the same evaluation protocol (APP) on the test set
+mae_score = qp.evaluation.evaluate(model, protocol=UPP(test), error_metric='mae')
+
+print(f'MAE={mae_score:.5f}')
+print(f'model selection took {took:.1f}s')
+
+
+model = KDEyML(bandwidth='auto')
+
+with qp.util.temp_seed(0):
+
+    param_grid = {
+        'classifier__C': hyper_C,
+    }
+
+    model = qp.model_selection.GridSearchQ(
+        model=model,
+        param_grid=param_grid,
+        protocol=protocol,
+        error='mae',  # the error to optimize is the MAE (a quantification-oriented loss)
+        refit=False,  # retrain on the whole labelled set once done
+        n_jobs=-1,
+        verbose=True  # show information as the process goes on
+    ).fit(training)
+
+best_params = model.best_params_
+took = model.fit_time_
+model = model.best_model_
+bandwidth = model.bandwidth_val
+print(f'model selection ended: best hyper-parameters={best_params} ({bandwidth=})')
+
+# evaluation in terms of MAE
+# we use the same evaluation protocol (APP) on the test set
+mae_score = qp.evaluation.evaluate(model, protocol=UPP(test), error_metric='mae')
+
+print(f'MAE={mae_score:.5f}')
+print(f'model selection took {took:.1f}s')
+

examples/16.confidence_regions.py

@@ -0,0 +1,81 @@
+from quapy.method.confidence import AggregativeBootstrap
+from quapy.method.aggregative import PACC
+import quapy.functional as F
+import quapy as qp
+
+"""
+Just like any other type of estimator, quantifier predictions are affected by error. It is therefore useful to provide,
+along with the point estimate (the class prevalence values) a measure of uncertainty. These, typically come in the 
+form of credible regions around the point estimate. 
+
+QuaPy implements a method for deriving confidence regions around point estimates of class prevalence based on bootstrap.
+
+Bootstrap method comes down to resampling the population several times, thus generating a series of point estimates.
+QuaPy provides a variant of bootstrap for aggregative quantifiers, that only applies resampling to the pre-classified
+instances.
+
+Let see one example:  
+"""
+
+# load some data
+data = qp.datasets.fetch_UCIMulticlassDataset('molecular')
+train, test = data.train_test
+Xtr, ytr = train.Xy
+
+# by simply wrapping an aggregative quantifier within the AggregativeBootstrap class, we can obtain confidence
+# intervals around the point estimate, in this case, at 95% of confidence
+pacc = AggregativeBootstrap(PACC(), n_test_samples=500, confidence_level=0.95)
+
+with qp.util.temp_seed(0):
+    # we train the quantifier the usual way
+    pacc.fit(Xtr, ytr)
+
+    # let us simulate some shift in the test data
+    random_prevalence = F.uniform_prevalence_sampling(n_classes=test.n_classes)
+    shifted_test = test.sampling(200, *random_prevalence)
+    true_prev = shifted_test.prevalence()
+
+    # by calling "quantify_conf", we obtain the point estimate and the confidence intervals around it
+    pred_prev, conf_intervals = pacc.quantify_conf(shifted_test.X)
+
+    # conf_intervals is an instance of ConfidenceRegionABC, which provides some useful utilities like:
+    # - coverage: a function which computes the fraction of true values that belong to the confidence region
+    # - simplex_proportion: estimates the proportion of the simplex covered by the confidence region (amplitude)
+    # ideally, we are interested in obtaining confidence regions with high level of coverage and small amplitude
+
+    # the point estimate is computed as the mean of all bootstrap predictions; let us see the prediction error
+    error = qp.error.ae(true_prev, pred_prev)
+
+    # some useful outputs
+    print(f'train prevalence: {F.strprev(train.prevalence())}')
+    print(f'test prevalence:  {F.strprev(true_prev)}')
+    print(f'point-estimate:   {F.strprev(pred_prev)}')
+    print(f'absolute error:   {error:.3f}')
+    print(f'Is the true value in the confidence region?: {conf_intervals.coverage(true_prev)==1}')
+    print(f'Proportion of simplex covered at confidence level {pacc.confidence_level*100:.1f}%: {conf_intervals.simplex_portion()*100:.2f}%')
+
+"""
+Final remarks: 
+There are various ways for performing bootstrap:
+- the population-based approach (default): performs resampling of the test instances
+    e.g., use  AggregativeBootstrap(PACC(), n_train_samples=1, n_test_samples=100, confidence_level=0.95)
+- the model-based approach: performs resampling of the training instances, thus training several quantifiers
+    e.g., use  AggregativeBootstrap(PACC(), n_train_samples=100, n_test_samples=1, confidence_level=0.95)
+    this implementation avoids retraining the classifier, and performs resampling only to train different aggregation functions 
+- the combined approach: a combination of the above
+    e.g., use  AggregativeBootstrap(PACC(), n_train_samples=100, n_test_samples=100, confidence_level=0.95)
+    this example will generate 100 x 100 predictions
+    
+There are different ways for constructing confidence regions implemented in QuaPy:
+- confidence intervals: the simplest way, and one that typically works well in practice
+    use: AggregativeBootstrap(PACC(), confidence_level=0.95, method='intervals')
+- confidence ellipse in the simplex: creates an ellipse, which lies on the probability simplex, around the point estimate
+    use: AggregativeBootstrap(PACC(), confidence_level=0.95, method='ellipse')
+- confidence ellipse in the Centered-Log Ratio (CLR) space: creates an ellipse in the CLR space (this should be 
+    convenient for taking into account the inner structure of the probability simplex)
+    use: AggregativeBootstrap(PACC(), confidence_level=0.95, method='ellipse-clr')
+    
+Other methods that return confidence regions in QuaPy include the BayesianCC method.
+"""
+
+

examples/17.explicit_loss_minimization.py

@@ -50,7 +50,7 @@ train_modsel, val = qp.datasets.fetch_twitter('hcr', for_model_selection=True, p
 model selection: 
 We explore the classifier's loss and the classifier's C hyperparameters.
 Since our model is actually an instance of OneVsAllAggregative, we need to add the prefix "binary_quantifier", and
-since our binary quantifier is an instance of CC, we need to add the prefix "classifier".
+since our binary quantifier is an instance of CC (an aggregative quantifier), we need to add the prefix "classifier".
 """
 param_grid = {
     'binary_quantifier__classifier__loss': ['q', 'kld', 'mae'],  # classifier-dependent hyperparameter
@@ -58,11 +58,11 @@ param_grid = {
 }
 print('starting model selection')
 model_selection = GridSearchQ(quantifier, param_grid, protocol=UPP(val), verbose=True, refit=False)
-quantifier = model_selection.fit(train_modsel).best_model()
+quantifier = model_selection.fit(*train_modsel.Xy).best_model()
 
 print('training on the whole training set')
 train, test = qp.datasets.fetch_twitter('hcr', for_model_selection=False, pickle=True).train_test
-quantifier.fit(train)
+quantifier.fit(*train.Xy)
 
 # evaluation
 mae = qp.evaluation.evaluate(quantifier, protocol=UPP(test), error_metric='mae')

examples/2.custom_quantifier.py

@@ -4,6 +4,7 @@ from quapy.method.base import BinaryQuantifier, BaseQuantifier
 from quapy.model_selection import GridSearchQ
 from quapy.method.aggregative import AggregativeSoftQuantifier
 from quapy.protocol import APP
+import quapy.functional as F
 import numpy as np
 from sklearn.linear_model import LogisticRegression
 from time import time
@@ -30,19 +31,19 @@ class MyQuantifier(BaseQuantifier):
         self.alpha = alpha
         self.classifier = classifier
 
-    # in general, we would need to implement the method fit(self, data: LabelledCollection, fit_classifier=True,
+    # in general, we would need to implement the method fit(self, X, y); this would amount to:
-    # val_split=None); this would amount to:
+    def fit(self, X, y):
-    def fit(self, data: LabelledCollection):
+        n_classes = F.num_classes_from_labels(y)
-        assert data.n_classes==2, \
+        assert n_classes==2, \
             'this quantifier is only valid for binary problems [abort]'
-        self.classifier.fit(*data.Xy)
+        self.classifier.fit(X, y)
         return self
 
     # in general, we would need to implement the method quantify(self, instances); this would amount to:
-    def quantify(self, instances):
+    def predict(self, X):
         assert hasattr(self.classifier, 'predict_proba'), \
             'the underlying classifier is not probabilistic! [abort]'
-        posterior_probabilities = self.classifier.predict_proba(instances)
+        posterior_probabilities = self.classifier.predict_proba(X)
         positive_probabilities = posterior_probabilities[:, 1]
         crisp_decisions = positive_probabilities > self.alpha
         pos_prev = crisp_decisions.mean()
@@ -57,9 +58,11 @@ class MyQuantifier(BaseQuantifier):
 # of the method, now adhering to the AggregativeSoftQuantifier:
 
 class MyAggregativeSoftQuantifier(AggregativeSoftQuantifier, BinaryQuantifier):
+
     def __init__(self, classifier, alpha=0.5):
-        # aggregative quantifiers have an internal attribute called self.classifier
+        # aggregative quantifiers have an internal attribute called self.classifier, but this is defined
-        self.classifier = classifier
+        # within the super's init
+        super().__init__(classifier, fit_classifier=True, val_split=None)
         self.alpha = alpha
 
     # since this method is of type aggregative, we can simply implement the method aggregation_fit, which
@@ -68,7 +71,7 @@ class MyAggregativeSoftQuantifier(AggregativeSoftQuantifier, BinaryQuantifier):
     # k-fold cross validation strategy). What remains ahead is to learn an aggregation function. In our case
     # this amounts to doing... nothing, since our method was pretty basic. BinaryQuantifier also add some
     # basic functionality for checking binary consistency.
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
         pass
 
     # since this method is of type aggregative, we can simply implement the method aggregate (i.e., we should
@@ -94,7 +97,7 @@ if __name__ == '__main__':
     train, test = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=5).train_test
     train, val = train.split_stratified(train_prop=0.75)  # let's create a validation set for optimizing hyperparams
 
-    def test_implementation(quantifier):
+    def try_implementation(quantifier):
         class_name = quantifier.__class__.__name__
         print(f'\ntesting implementation {class_name}...')
         # model selection
@@ -104,7 +107,7 @@ if __name__ == '__main__':
             'alpha': np.linspace(0, 1, 11),         # quantifier-dependent hyperparameter
             'classifier__C': np.logspace(-2, 2, 5)  # classifier-dependent hyperparameter
         }
-        gridsearch = GridSearchQ(quantifier, param_grid, protocol=APP(val), n_jobs=-1, verbose=False).fit(train)
+        gridsearch = GridSearchQ(quantifier, param_grid, protocol=APP(val), n_jobs=-1, verbose=True).fit(*train.Xy)
         t_modsel = time() - tinit
         print(f'\tmodel selection took {t_modsel:.2f}s', flush=True)
 
@@ -112,7 +115,7 @@ if __name__ == '__main__':
         optimized_model = gridsearch.best_model_
         mae = qp.evaluation.evaluate(
             optimized_model,
-            protocol=APP(test, repeats=5000, sanity_check=None),  # disable the check, we want to generate many tests!
+            protocol=APP(test, repeats=500, sanity_check=None),  # disable the check, we want to generate many tests!
             error_metric='mae',
             verbose=True)
 
@@ -121,11 +124,11 @@ if __name__ == '__main__':
 
     # define an instance of our custom quantifier and test it!
     quantifier = MyQuantifier(LogisticRegression(), alpha=0.5)
-    test_implementation(quantifier)
+    try_implementation(quantifier)
 
     # define an instance of our custom quantifier, with the second implementation, and test it!
     quantifier = MyAggregativeSoftQuantifier(LogisticRegression(), alpha=0.5)
-    test_implementation(quantifier)
+    try_implementation(quantifier)
 
     # the output should look like this:
     """
@@ -141,7 +144,7 @@ if __name__ == '__main__':
         evaluation took 4.66s [MAE = 0.0630]
     """
     # Note that the first implementation is much slower, both in terms of grid-search optimization and in terms of
-    # evaluation. The reason why is that QuaPy is highly optimized for aggregative quantifiers (by far, the most
+    # evaluation. The reason why, is that QuaPy is highly optimized for aggregative quantifiers (by far, the most
     # popular type of quantification methods), thus significantly speeding up model selection and test routines.
     # Furthermore, it is simpler to extend an aggregation type since QuaPy implements boilerplate functions for you.
 

examples/3.custom_collection.py

@@ -0,0 +1,103 @@
+import quapy as qp
+from quapy.method.aggregative import PACC
+from quapy.data import LabelledCollection, Dataset
+from quapy.protocol import ArtificialPrevalenceProtocol
+import quapy.functional as F
+import os
+from os.path import join
+
+# While quapy comes with ready-to-use datasets for experimental purposes, you may prefer to run experiments using
+# your own data. Most of the quapy's functionality relies on an internal class called LabelledCollection, for fast
+# indexing and sampling, and so this example provides guidance on how to convert your datasets into a LabelledCollection
+# so all the functionality becomes available. This includes procedures for tuning the hyperparameters of your methods,
+# evaluating the performance using high level sampling protocols, etc.
+
+# Let us assume that we have a binary sentiment dataset of opinions in natural language. We will use the "IMDb"
+# dataset of reviews, which can be downloaded as follows
+URL_TRAIN = f'https://zenodo.org/record/4117827/files/imdb_train.txt'
+URL_TEST = f'https://zenodo.org/record/4117827/files/imdb_test.txt'
+os.makedirs('./reviews', exist_ok=True)
+train_path = join('reviews', 'hp_train.txt')
+test_path = join('reviews', 'hp_test.txt')
+qp.util.download_file_if_not_exists(URL_TRAIN, train_path)
+qp.util.download_file_if_not_exists(URL_TEST, test_path)
+
+# these files contain 2 columns separated by a \t:
+# the first one is a binary value (0=negative, 1=positive), and the second is the text
+# Everything we need is to implement a function returning the instances and the labels as follows
+def my_data_loader(path):
+    with open(path, 'rt') as fin:
+        labels, texts = zip(*[line.split('\t') for line in fin.readlines()])
+        labels = list(map(int, labels))  # convert string numbers to int
+        return texts, labels
+
+# check that our function is working properly...
+train_texts, train_labels = my_data_loader(train_path)
+for i, (text, label) in enumerate(zip(train_texts, train_labels)):
+    print(f'#{i}: {label=}\t{text=}')
+    if i>=5:
+        print('...')
+        break
+
+# We can now instantiate a LabelledCollection simply as
+train_lc = LabelledCollection(instances=train_texts, labels=train_labels)
+print('my training collection:', train_lc)
+
+# We can instantiate directly a LabelledCollection using the data loader function,
+# without having to load the data ourselves:
+train_lc = LabelledCollection.load(train_path, loader_func=my_data_loader)
+print('my training collection:', train_lc)
+
+# We can do the same for the test set, or we can instead directly instantiate a Dataset object (this is by and large
+# simply a tuple with training and test LabelledCollections) as follows:
+my_data = Dataset.load(train_path, test_path, loader_func=my_data_loader)
+print('my dataset:', my_data)
+
+# However, since this is a textual dataset, we must vectorize it prior to training any quantification algorithm.
+# We can do this in several ways in quapy. For example, manually...
+# from sklearn.feature_extraction.text import TfidfVectorizer
+# tfidf = TfidfVectorizer(min_df=5)
+# Xtr = tfidf.fit_transform(my_data.training.instances)
+# Xte = tfidf.transform(my_data.test.instances)
+# ... or using some preprocessing functionality of quapy (recommended):
+my_data_tfidf = qp.data.preprocessing.text2tfidf(my_data, min_df=5)
+
+training, test = my_data_tfidf.train_test
+
+# Once you have loaded your training and test data, you have access to a series of quapy's utilities, e.g.:
+print(f'the training prevalence is {F.strprev(training.prevalence())}')
+print(f'the test prevalence is {F.strprev(test.prevalence())}')
+print(f'let us generate a small balanced training sample:')
+desired_size = 200
+desired_prevalence = [0.5, 0.5]
+small_training_balanced = training.sampling(desired_size, *desired_prevalence, shuffle=True, random_state=0)
+print(small_training_balanced)
+print(f'or generating train/val splits such as: {training.split_stratified(train_prop=0.7)}')
+
+# training
+print('let us train a simple quantifier:...')
+Xtr, ytr = training.Xy
+quantifier = PACC()
+quantifier.fit(Xtr, ytr)  # or: quantifier.fit(*training.Xy)
+
+# test
+print("and use quapy' evaluation functions")
+evaluation_protocol = ArtificialPrevalenceProtocol(
+    data=test,
+    sample_size=200,
+    random_state=0
+)
+
+report = qp.evaluation.evaluation_report(quantifier, protocol=evaluation_protocol, error_metrics=['ae'])
+print(report)
+print(f'mean absolute error across {len(report)} experiments: {report.mean(numeric_only=True)}')
+
+
+
+
+
+
+
+
+
+

examples/4.using_pretrained_classifier.py

@@ -0,0 +1,75 @@
+"""
+Aggregative quantifiers use an underlying classifier. Often, one has one pre-trained classifier available, and
+needs to use this classifier at the basis of a quantification system. In such cases, the classifier should not
+be retrained, but only used to issue classifier predictions for the quantifier.
+In this example, we show how to instantiate a quantifier with a pre-trained classifier.
+"""
+from typing import List, Dict
+
+import quapy as qp
+from quapy.method.aggregative import PACC
+from sklearn.base import BaseEstimator, ClassifierMixin
+from transformers import pipeline
+import numpy as np
+import quapy.functional as F
+
+
+# A scikit-learn's style wrapper for a huggingface-based pre-trained transformer for binary sentiment classification
+class HFTextClassifier(BaseEstimator, ClassifierMixin):
+    def __init__(self, model_name='distilbert-base-uncased-finetuned-sst-2-english'):
+        self.pipe = pipeline("sentiment-analysis", model=model_name)
+        self.classes_ = np.asarray([0,1])
+
+    def fit(self, X, y=None):
+        return self
+
+    def _binary_decisions(self, transformer_output: List[Dict]):
+        return np.array([(1 if p['label']=='POSITIVE' else 0) for p in transformer_output], dtype=int)
+
+    def predict(self, X):
+        X = list(map(str, X))
+        preds = self.pipe(X, truncation=True)
+        return self._binary_decisions(preds)
+
+    def predict_proba(self, X):
+        X = list(map(str, X))
+        n_examples = len(X)
+        preds = self.pipe(X, truncation=True)
+        decisions = self._binary_decisions(preds)
+        scores = np.array([p['score'] for p in preds], dtype=float)
+        probas = np.zeros(shape=(len(X), 2), dtype=float)
+        probas[np.arange(n_examples),decisions] = scores
+        probas[np.arange(n_examples),~decisions] = 1-scores
+        return probas
+
+# load a sentiment dataset
+dataset = qp.datasets.fetch_reviews('imdb', tfidf=False)  # raw text
+train, test = dataset.training, dataset.test
+
+# instantiate a pre-trained classifier
+clf = HFTextClassifier()
+
+# Let us fit a quantifier based on our pre-trained classifier.
+# Note that, since the classifier is already fit, we will use the entire training set for
+# learning the aggregation function of the quantifier.
+# To do so, we only need to indicate "fit_classifier"=False, as follows:
+quantifier = PACC(clf, fit_classifier=False)   # Probabilistic Classify & Count using a pre-trained model
+
+print('training PACC...')
+quantifier.fit(*train.Xy)
+
+# let us simulate some shifted test data...
+new_prevalence = [0.75, 0.25]
+shifted_test = test.sampling(500, *new_prevalence, random_state=0)
+
+# and do some evaluation
+print('predicting with PACC...')
+estim_prevalence = quantifier.predict(shifted_test.X)
+
+print('Result:\n'+('='*20))
+print(f'training prevalence: {F.strprev(train.prevalence())}')
+print(f'(shifted) test prevalence: {F.strprev(shifted_test.prevalence())}')
+print(f'estimated prevalence: {F.strprev(estim_prevalence)}')
+
+absolute_error = qp.error.ae(new_prevalence, estim_prevalence)
+print(f'absolute error={absolute_error:.4f}')

examples/5a.lequa2022_experiments.py

@@ -15,7 +15,7 @@ https://lequa2022.github.io/index (the site of the competition)
 https://ceur-ws.org/Vol-3180/paper-146.pdf (the overview paper)
 """
 
-# there are 4 tasks (T1A, T1B, T2A, T2B)
+# there are 4 tasks (T1A, T1B, T2A, T2B), let us symply consider T1A (binary quantification, vector form)
 task = 'T1A'
 
 # set the sample size in the environment. The sample size is task-dendendent and can be consulted by doing:
@@ -28,18 +28,19 @@ qp.environ['N_JOBS'] = -1
 # of SamplesFromDir, a protocol that simply iterates over pre-generated samples (those provided for the competition)
 # stored in a directory.
 training, val_generator, test_generator = fetch_lequa2022(task=task)
+Xtr, ytr = training.Xy
 
 # define the quantifier
-quantifier = EMQ(classifier=LogisticRegression())
+quantifier = EMQ(classifier=LogisticRegression(), val_split=5)
 
 # model selection
 param_grid = {
     'classifier__C': np.logspace(-3, 3, 7),          # classifier-dependent: inverse of regularization strength
     'classifier__class_weight': ['balanced', None],  # classifier-dependent: weights of each class
-    'recalib': ['bcts', 'platt', None]               # quantifier-dependent: recalibration method (new in v0.1.7)
+    'calib': ['bcts', None]                          # quantifier-dependent: recalibration method (new in v0.1.7)
 }
 model_selection = GridSearchQ(quantifier, param_grid, protocol=val_generator, error='mrae', refit=False, verbose=True)
-quantifier = model_selection.fit(training)
+quantifier = model_selection.fit(Xtr, ytr)
 
 # evaluation
 report = evaluation_report(quantifier, protocol=test_generator, error_metrics=['mae', 'mrae', 'mkld'], verbose=True)
@@ -50,4 +51,4 @@ report['estim-prev'] = report['estim-prev'].map(F.strprev)
 print(report)
 
 print('Averaged values:')
-print(report.mean())
+print(report.mean(numeric_only=True))

examples/5b.lequa2024_experiments.py

@@ -1,6 +1,6 @@
+import quapy as qp
 import numpy as np
 from sklearn.linear_model import LogisticRegression
-import quapy as qp
 import quapy.functional as F
 from quapy.data.datasets import LEQUA2024_SAMPLE_SIZE, fetch_lequa2024
 from quapy.evaluation import evaluation_report
@@ -14,6 +14,7 @@ LeQua competition itself, check:
 https://lequa2024.github.io/index (the site of the competition)
 """
 
+
 # there are 4 tasks: T1 (binary), T2 (multiclass), T3 (ordinal), T4 (binary - covariate & prior shift)
 task = 'T2'
 
@@ -27,6 +28,7 @@ qp.environ['N_JOBS'] = -1
 # of SamplesFromDir, a protocol that simply iterates over pre-generated samples (those provided for the competition)
 # stored in a directory.
 training, val_generator, test_generator = fetch_lequa2024(task=task)
+Xtr, ytr = training.Xy
 
 # define the quantifier
 quantifier = KDEyML(classifier=LogisticRegression())
@@ -37,8 +39,9 @@ param_grid = {
     'classifier__class_weight': ['balanced', None],         # classifier-dependent: weights of each class
     'bandwidth': np.linspace(0.01, 0.2, 20)  # quantifier-dependent: bandwidth of the kernel
 }
+
 model_selection = GridSearchQ(quantifier, param_grid, protocol=val_generator, error='mrae', refit=False, verbose=True)
-quantifier = model_selection.fit(training)
+quantifier = model_selection.fit(Xtr, ytr)
 
 # evaluation
 report = evaluation_report(quantifier, protocol=test_generator, error_metrics=['mae', 'mrae'], verbose=True)

examples/6.quanet_example.py

@@ -20,14 +20,13 @@ train, test = dataset.train_test
 # train the text classifier:
 cnn_module = CNNnet(dataset.vocabulary_size, dataset.training.n_classes)
 cnn_classifier = NeuralClassifierTrainer(cnn_module, device='cuda')
-cnn_classifier.fit(*dataset.training.Xy)
 
 # train QuaNet (alternatively, we can set fit_classifier=True and let QuaNet train the classifier)
 quantifier = QuaNet(cnn_classifier, device='cuda')
-quantifier.fit(train, fit_classifier=False)
+quantifier.fit(*train.Xy)
 
 # prediction and evaluation
-estim_prevalence = quantifier.quantify(test.instances)
+estim_prevalence = quantifier.predict(test.instances)
 mae = qp.error.mae(test.prevalence(), estim_prevalence)
 
 print(f'true prevalence: {F.strprev(test.prevalence())}')

examples/7.uci_binary_experiments.py

@@ -1,4 +1,7 @@
 from copy import deepcopy
+from pathlib import Path
+
+import pandas as pd
 
 import quapy as qp
 from sklearn.calibration import CalibratedClassifierCV
@@ -15,6 +18,18 @@ import itertools
 import argparse
 import torch
 import shutil
+from glob import glob
+
+
+"""
+This example shows how to generate experiments for the UCI ML repository binary datasets following the protocol
+proposed in "Pérez-Gállego , P., Quevedo , J. R., and del Coz, J. J. Using ensembles for problems with characteriz-
+able changes in data distribution: A case study on quantification. Information Fusion 34 (2017), 87–100."
+
+This example covers most important steps in the experimentation pipeline, namely, the training and optimization
+of the hyperparameters of different quantifiers, and the evaluation of these quantifiers based on standard 
+prevalence sampling protocols aimed at simulating different levels of prior probability shift.
+"""
 
 
 N_JOBS = -1
@@ -28,10 +43,6 @@ def newLR():
     return LogisticRegression(max_iter=1000, solver='lbfgs', n_jobs=-1)
 
 
-def calibratedLR():
-    return CalibratedClassifierCV(newLR())
-
-
 __C_range = np.logspace(-3, 3, 7)
 lr_params = {
     'classifier__C': __C_range,
@@ -50,7 +61,7 @@ def quantification_models():
     yield 'MAX', MAX(newLR()), lr_params
     yield 'MS', MS(newLR()), lr_params
     yield 'MS2', MS2(newLR()), lr_params
-    yield 'sldc', EMQ(newLR(), recalib='platt'), lr_params
+    yield 'sldc', EMQ(newLR()), lr_params
     yield 'svmmae', newSVMAE(), svmperf_params
     yield 'hdy', HDy(newLR()), lr_params
 
@@ -74,6 +85,13 @@ def result_path(path, dataset_name, model_name, run, optim_loss):
     return os.path.join(path, f'{dataset_name}-{model_name}-run{run}-{optim_loss}.pkl')
 
 
+def parse_result_path(path):
+    *dataset, method, run, metric = Path(path).name.split('-')
+    dataset = '-'.join(dataset)
+    run = int(run.replace('run',''))
+    return dataset, method, run, metric
+
+
 def is_already_computed(dataset_name, model_name, run, optim_loss):
     return os.path.exists(result_path(args.results, dataset_name, model_name, run, optim_loss))
 
@@ -98,8 +116,8 @@ def run(experiment):
         print(f'running dataset={dataset_name} model={model_name} loss={optim_loss} run={run+1}/5')
         # model selection (hyperparameter optimization for a quantification-oriented loss)
         train, test = data.train_test
-        train, val = train.split_stratified()
         if hyperparams is not None:
+            train, val = train.split_stratified()
             model_selection = qp.model_selection.GridSearchQ(
                 deepcopy(model),
                 param_grid=hyperparams,
@@ -107,13 +125,13 @@ def run(experiment):
                 error=optim_loss,
                 refit=True,
                 timeout=60*60,
-                verbose=True
+                verbose=False
             )
-            model_selection.fit(train)
+            model_selection.fit(*train.Xy)
             model = model_selection.best_model()
             best_params = model_selection.best_params_
         else:
-            model.fit(data.training)
+            model.fit(*train.Xy)
             best_params = {}
 
         # model evaluation
@@ -121,19 +139,37 @@ def run(experiment):
             model,
             protocol=APP(test, n_prevalences=21, repeats=100)
         )
-        test_true_prevalence = data.test.prevalence()
+        test_true_prevalence = test.prevalence()
 
         evaluate_experiment(true_prevalences, estim_prevalences)
         save_results(dataset_name, model_name, run, optim_loss,
                      true_prevalences, estim_prevalences,
-                     data.training.prevalence(), test_true_prevalence,
+                     train.prevalence(), test_true_prevalence,
                      best_params)
 
 
+def show_results(result_folder):
+    result_data = []
+    for file in glob(os.path.join(result_folder,'*.pkl')):
+        true_prevalences, estim_prevalences, *_ = pickle.load(open(file, 'rb'))
+        dataset, method, run, metric = parse_result_path(file)
+        mae = qp.error.mae(true_prevalences, estim_prevalences)
+        result_data.append({
+            'dataset': dataset,
+            'method': method,
+            'run': run,
+            metric: mae
+        })
+    df = pd.DataFrame(result_data)
+    pd.set_option("display.max_columns", None)
+    pd.set_option("display.expand_frame_repr", False)
+    print(df.pivot_table(index='dataset', columns='method', values=metric))
+
+
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='Run experiments for Tweeter Sentiment Quantification')
-    parser.add_argument('results', metavar='RESULT_PATH', type=str,
+    parser.add_argument('--results', metavar='RESULT_PATH', type=str,
-                        help='path to the directory where to store the results')
+                        help='path to the directory where to store the results', default='./results/uci_binary')
     parser.add_argument('--svmperfpath', metavar='SVMPERF_PATH', type=str, default='../svm_perf_quantification',
                         help='path to the directory with svmperf')
     parser.add_argument('--checkpointdir', metavar='PATH', type=str, default='./checkpoint',
@@ -155,3 +191,5 @@ if __name__ == '__main__':
     qp.util.parallel(run, itertools.product(optim_losses, datasets, models), n_jobs=CUDA_N_JOBS)
 
     shutil.rmtree(args.checkpointdir, ignore_errors=True)
+
+    show_results(args.results)

examples/8.uci_multiclass_experiments.py

@@ -1,4 +1,3 @@
-import pickle
 import os
 from time import time
 from collections import defaultdict
@@ -7,11 +6,16 @@ import numpy as np
 from sklearn.linear_model import LogisticRegression
 
 import quapy as qp
-from quapy.method.aggregative import PACC, EMQ
+from quapy.method.aggregative import PACC, EMQ, KDEyML
 from quapy.model_selection import GridSearchQ
 from quapy.protocol import UPP
 from pathlib import Path
 
+"""
+This example is the analogous counterpart of example 7 but involving multiclass quantification problems
+using datasets from the UCI ML repository.
+"""
+
 
 SEED = 1
 
@@ -31,7 +35,7 @@ def wrap_hyper(classifier_hyper_grid:dict):
 METHODS = [
     ('PACC', PACC(newLR()), wrap_hyper(logreg_grid)),
     ('EMQ',  EMQ(newLR()), wrap_hyper(logreg_grid)),
-    # ('KDEy-ML',  KDEyML(newLR()), {**wrap_hyper(logreg_grid), **{'bandwidth': np.linspace(0.01, 0.2, 20)}}),
+    ('KDEy-ML',  KDEyML(newLR()), {**wrap_hyper(logreg_grid), **{'bandwidth': np.linspace(0.01, 0.2, 20)}}),
 ]
 
 
@@ -43,6 +47,7 @@ def show_results(result_path):
     pv = df.pivot_table(index='Dataset', columns="Method", values=["MAE", "MRAE", "t_train"], margins=True)
     print(pv)
 
+
 def load_timings(result_path):
     import pandas as pd
     timings = defaultdict(lambda: {})
@@ -59,7 +64,7 @@ if __name__ == '__main__':
     qp.environ['N_JOBS'] = -1
     n_bags_val = 250
     n_bags_test = 1000
-    result_dir = f'results/ucimulti'
+    result_dir = f'results/uci_multiclass'
 
     os.makedirs(result_dir, exist_ok=True)
 
@@ -100,7 +105,7 @@ if __name__ == '__main__':
                         
                         t_init = time()
                         try:
-                            modsel.fit(train)
+                            modsel.fit(*train.Xy)
 
                             print(f'best params {modsel.best_params_}')
                             print(f'best score {modsel.best_score_}')
@@ -108,7 +113,8 @@ if __name__ == '__main__':
                             quantifier = modsel.best_model()
                         except:
                             print('something went wrong... trying to fit the default model')
-                            quantifier.fit(train)
+                            quantifier.fit(*train.Xy)
+
                         timings[method_name][dataset] = time() - t_init
                         
 

examples/9.ifcb_experiments.py

@@ -6,6 +6,18 @@ from sklearn.linear_model import LogisticRegression
 from quapy.model_selection import GridSearchQ
 from quapy.evaluation import evaluation_report
 
+"""
+This example shows a complete experiment using the IFCB Plankton dataset;
+see https://hlt-isti.github.io/QuaPy/manuals/datasets.html#ifcb-plankton-dataset
+
+Note that this dataset can be downloaded in two modes: for model selection or for evaluation.
+
+See also:
+Automatic plankton quantification using deep features
+P González, A Castaño, EE Peacock, J Díez, JJ Del Coz, HM Sosik
+Journal of Plankton Research 41 (4), 449-463
+"""
+
 
 print('Quantifying the IFCB dataset with PACC\n')
 
@@ -30,7 +42,7 @@ mod_sel = GridSearchQ(
     n_jobs=-1,
     verbose=True,
     raise_errors=True
-).fit(train)
+).fit(*train.Xy)
 
 print(f'model selection chose hyperparameters: {mod_sel.best_params_}')
 quantifier = mod_sel.best_model_
@@ -42,7 +54,7 @@ print(f'\ttraining size={len(train)}, features={train.X.shape[1]}, classes={trai
 print(f'\ttest samples={test_gen.total()}')
 
 print('training on the whole dataset before test')
-quantifier.fit(train)
+quantifier.fit(*train.Xy)
 
 print('testing...')
 report = evaluation_report(quantifier, protocol=test_gen, error_metrics=['mae'], verbose=True)

examples/distributing_samples.py

@@ -0,0 +1,38 @@
+"""
+Imagine we want to generate many samples out of a collection, that we want to distribute for others to run their
+own experiments in the very same test samples. One naive solution would come down to applying a given protocol to
+our collection (say the artificial prevalence protocol on the 'academic-success' UCI dataset), store all those samples
+on disk and make them available online. Distributing many such samples is undesirable.
+In this example, we generate the indexes that allow anyone to regenerate the samples out of the original collection.
+"""
+
+import quapy as qp
+from quapy.method.aggregative import PACC
+from quapy.protocol import UPP
+
+data = qp.datasets.fetch_UCIMulticlassDataset('academic-success')
+train, test = data.train_test
+
+# let us train a quantifier to check whether we can actually replicate the results
+quantifier = PACC()
+quantifier.fit(train)
+
+# let us simulate our experimental results
+protocol = UPP(test, sample_size=100, repeats=100, random_state=0)
+our_mae = qp.evaluation.evaluate(quantifier, protocol=protocol, error_metric='mae')
+
+print(f'We have obtained a MAE={our_mae:.3f}')
+
+# let us distribute the indexes; we specify that we want the indexes, not the samples
+protocol = UPP(test, sample_size=100, repeats=100, random_state=0, return_type='index')
+indexes = protocol.samples_parameters()
+
+# Imagine we distribute the indexes; now we show how to replicate our experiments.
+from quapy.protocol import ProtocolFromIndex
+data = qp.datasets.fetch_UCIMulticlassDataset('academic-success')
+train, test = data.train_test
+protocol = ProtocolFromIndex(data=test, indexes=indexes)
+their_mae = qp.evaluation.evaluate(quantifier, protocol=protocol, error_metric='mae')
+
+print(f'Another lab obtains a MAE={our_mae:.3f}')
+

examples/ensembles.py

@@ -0,0 +1,56 @@
+from sklearn.exceptions import ConvergenceWarning
+from sklearn.linear_model import LogisticRegression
+from sklearn.naive_bayes import MultinomialNB
+from sklearn.neighbors import KNeighborsClassifier
+from statsmodels.sandbox.distributions.genpareto import quant
+
+import quapy as qp
+from quapy.protocol import UPP
+from quapy.method.aggregative import PACC, DMy, EMQ, KDEyML
+from quapy.method.meta import SCMQ, MCMQ, MCSQ
+import warnings
+warnings.filterwarnings("ignore", category=DeprecationWarning)
+warnings.filterwarnings("ignore", category=ConvergenceWarning)
+
+qp.environ["SAMPLE_SIZE"]=100
+
+def train_and_test_model(quantifier, train, test):
+    quantifier.fit(train)
+    report = qp.evaluation.evaluation_report(quantifier, UPP(test), error_metrics=['mae', 'mrae'])
+    print(quantifier.__class__.__name__)
+    print(report.mean(numeric_only=True))
+
+
+quantifiers = [
+    PACC(),
+    DMy(),
+    EMQ(),
+    KDEyML()
+]
+
+classifier = LogisticRegression()
+
+dataset_name = qp.datasets.UCI_MULTICLASS_DATASETS[0]
+data = qp.datasets.fetch_UCIMulticlassDataset(dataset_name)
+train, test = data.train_test
+
+scmq = SCMQ(classifier, quantifiers)
+
+train_and_test_model(scmq, train, test)
+
+# for quantifier in quantifiers:
+#     train_and_test_model(quantifier, train, test)
+
+classifiers = [
+    LogisticRegression(),
+    KNeighborsClassifier(),
+    # MultinomialNB()
+]
+
+mcmq = MCMQ(classifiers, quantifiers)
+
+train_and_test_model(mcmq, train, test)
+
+mcsq = MCSQ(classifiers, PACC())
+
+train_and_test_model(mcsq, train, test)

prepare_svmperf.sh

@@ -11,13 +11,5 @@ rm $FILE
 patch -s -p0 < svm-perf-quantification-ext.patch
 mv svm_perf svm_perf_quantification
 cd svm_perf_quantification
-make
+make CFLAGS="-O3 -Wall -Wno-unused-result -fcommon"
-
-
-
-
-
-
-
-
 

quapy/__init__.py

@@ -1,5 +1,4 @@
 """QuaPy module for quantification"""
-from sklearn.linear_model import LogisticRegression
 
 from quapy.data import datasets
 from . import error
@@ -14,7 +13,13 @@ from . import model_selection
 from . import classification
 import os
 
-__version__ = '0.1.9'
+__version__ = '0.2.0'
+
+
+def _default_cls():
+    from sklearn.linear_model import LogisticRegression
+    return LogisticRegression()
+
 
 environ = {
     'SAMPLE_SIZE': None,
@@ -24,7 +29,7 @@ environ = {
     'PAD_INDEX': 1,
     'SVMPERF_HOME': './svm_perf_quantification',
     'N_JOBS': int(os.getenv('N_JOBS', 1)),
-    'DEFAULT_CLS': LogisticRegression(max_iter=3000)
+    'DEFAULT_CLS': _default_cls()
 }
 
 
@@ -68,3 +73,5 @@ def _get_classifier(classifier):
     if classifier is None:
         raise ValueError('neither classifier nor qp.environ["DEFAULT_CLS"] have been specified')
     return classifier
+
+

quapy/classification/svmperf.py

@@ -33,27 +33,16 @@ class SVMperf(BaseEstimator, ClassifierMixin):
     valid_losses = {'01':0, 'f1':1, 'kld':12, 'nkld':13, 'q':22, 'qacc':23, 'qf1':24, 'qgm':25, 'mae':26, 'mrae':27}
 
     def __init__(self, svmperf_base, C=0.01, verbose=False, loss='01', host_folder=None):
-        assert exists(svmperf_base), f'path {svmperf_base} does not seem to point to a valid path'
+        assert exists(svmperf_base), \
+            (f'path {svmperf_base} does not seem to point to a valid path;'
+             f'did you install svm-perf? '
+             f'see instructions in https://hlt-isti.github.io/QuaPy/manuals/explicit-loss-minimization.html')
         self.svmperf_base = svmperf_base
         self.C = C
         self.verbose = verbose
         self.loss = loss
         self.host_folder = host_folder
 
-    # def set_params(self, **parameters):
-    #     """
-    #     Set the hyper-parameters for svm-perf. Currently, only the `C` and `loss` parameters are supported
-    #
-    #     :param parameters: a `**kwargs` dictionary `{'C': <float>}`
-    #     """
-    #     assert sorted(list(parameters.keys())) == ['C', 'loss'], \
-    #         'currently, only the C and loss parameters are supported'
-    #     self.C = parameters.get('C', self.C)
-    #     self.loss = parameters.get('loss', self.loss)
-    #
-    # def get_params(self, deep=True):
-    #     return {'C': self.C, 'loss': self.loss}
-
     def fit(self, X, y):
         """
         Trains the SVM for the multivariate performance loss

quapy/data/base.py

@@ -9,6 +9,7 @@ from sklearn.model_selection import train_test_split, RepeatedStratifiedKFold
 from numpy.random import RandomState
 from quapy.functional import strprev
 from quapy.util import temp_seed
+import quapy.functional as F
 
 
 class LabelledCollection:
@@ -34,8 +35,7 @@ class LabelledCollection:
         self.labels = np.asarray(labels)
         n_docs = len(self)
         if classes is None:
-            self.classes_ = np.unique(self.labels)
+            self.classes_ = F.classes_from_labels(self.labels)
-            self.classes_.sort()
         else:
             self.classes_ = np.unique(np.asarray(classes))
             self.classes_.sort()
@@ -95,6 +95,15 @@ class LabelledCollection:
         """
         return len(self.classes_)
 
+    @property
+    def n_instances(self):
+        """
+        The number of instances
+
+        :return: integer
+        """
+        return len(self.labels)
+
     @property
     def binary(self):
         """
@@ -232,11 +241,11 @@ class LabelledCollection:
         :return: two instances of :class:`LabelledCollection`, the first one with `train_prop` elements, and the
             second one with `1-train_prop` elements
         """
-        tr_docs, te_docs, tr_labels, te_labels = train_test_split(
+        tr_X, te_X, tr_y, te_y = train_test_split(
             self.instances, self.labels, train_size=train_prop, stratify=self.labels, random_state=random_state
         )
-        training = LabelledCollection(tr_docs, tr_labels, classes=self.classes_)
+        training = LabelledCollection(tr_X, tr_y, classes=self.classes_)
-        test = LabelledCollection(te_docs, te_labels, classes=self.classes_)
+        test = LabelledCollection(te_X, te_y, classes=self.classes_)
         return training, test
 
     def split_random(self, train_prop=0.6, random_state=None):
@@ -318,6 +327,15 @@ class LabelledCollection:
         classes = np.unique(labels).sort()
         return LabelledCollection(instances, labels, classes=classes)
 
+    @property
+    def classes(self):
+        """
+        Gets an array-like with the classes used in this collection
+
+        :return: array-like
+        """
+        return self.classes_
+
     @property
     def Xy(self):
         """
@@ -414,6 +432,11 @@ class LabelledCollection:
             test = self.sampling_from_index(test_index)
             yield train, test
 
+    def __repr__(self):
+        repr=f'<{self.n_instances} instances (dtype={type(self.instances[0])}), '
+        repr+=f'n_classes={self.n_classes} {self.classes_}, prevalence={F.strprev(self.prevalence())}>'
+        return repr
+
 
 class Dataset:
     """
@@ -568,3 +591,6 @@ class Dataset:
             random_state = random_state
         )
         return self
+
+    def __repr__(self):
+        return f'training={self.training}; test={self.test}'

quapy/data/datasets.py

@@ -114,7 +114,8 @@ def fetch_reviews(dataset_name, tfidf=False, min_df=None, data_home=None, pickle
     """
     Loads a Reviews dataset as a Dataset instance, as used in
     `Esuli, A., Moreo, A., and Sebastiani, F. "A recurrent neural network for sentiment quantification."
-    Proceedings of the 27th ACM International Conference on Information and Knowledge Management. 2018. <https://dl.acm.org/doi/abs/10.1145/3269206.3269287>`_.
+    Proceedings of the 27th ACM International Conference on Information and Knowledge Management. 2018.
+    <https://dl.acm.org/doi/abs/10.1145/3269206.3269287>`_.
     The list of valid dataset names can be accessed in `quapy.data.datasets.REVIEWS_SENTIMENT_DATASETS`
 
     :param dataset_name: the name of the dataset: valid ones are 'hp', 'kindle', 'imdb'
@@ -499,7 +500,7 @@ def fetch_UCIBinaryLabelledCollection(dataset_name, data_home=None, standardize=
             y = df["NSP"].astype(int).values
         elif group == "semeion":
             with download_tmp_file("semeion", "semeion.data") as tmp:
-                df = pd.read_csv(tmp, header=None, delim_whitespace=True)
+                df = pd.read_csv(tmp, header=None, sep='\s+')
             X = df.iloc[:, 0:256].astype(float).values
             y = df[263].values  # 263 stands for digit 8 (labels are one-hot vectors from col 256-266)
         else:
@@ -548,25 +549,20 @@ def fetch_UCIBinaryLabelledCollection(dataset_name, data_home=None, standardize=
         """
         if name == "acute.a":
             X, y = data["X"], data["y"][:, 0]
-            # X, y = Xy[:, :-2], Xy[:, -2]
         elif name == "acute.b":
             X, y = data["X"], data["y"][:, 1]
-            # X, y = Xy[:, :-2], Xy[:, -1]
         elif name == "wine-q-red":
             X, y, color = data["X"], data["y"], data["color"]
-            # X, y, color = Xy[:, :-2], Xy[:, -2], Xy[:, -1]
             red_idx = color == "red"
             X, y = X[red_idx, :], y[red_idx]
             y = (y > 5).astype(int)
         elif name == "wine-q-white":
             X, y, color = data["X"], data["y"], data["color"]
-            # X, y, color = Xy[:, :-2], Xy[:, -2], Xy[:, -1]
             white_idx = color == "white"
             X, y = X[white_idx, :], y[white_idx]
             y = (y > 5).astype(int)
         else:
             X, y = data["X"], data["y"]
-            # X, y = Xy[:, :-1], Xy[:, -1]
 
         y = binarize(y, pos_class=pos_class[name])
 
@@ -797,7 +793,7 @@ def _array_replace(arr, repl={"yes": 1, "no": 0}):
 
 def fetch_lequa2022(task, data_home=None):
     """
-    Loads the official datasets provided for the `LeQua <https://lequa2022.github.io/index>`_ competition.
+    Loads the official datasets provided for the `LeQua 2022 <https://lequa2022.github.io/index>`_ competition.
     In brief, there are 4 tasks (T1A, T1B, T2A, T2B) having to do with text quantification
     problems. Tasks T1A and T1B provide documents in vector form, while T2A and T2B provide raw documents instead.
     Tasks T1A and T2A are binary sentiment quantification problems, while T2A and T2B are multiclass quantification
@@ -817,7 +813,7 @@ def fetch_lequa2022(task, data_home=None):
         ~/quay_data/ directory)
     :return: a tuple `(train, val_gen, test_gen)` where `train` is an instance of
         :class:`quapy.data.base.LabelledCollection`, `val_gen` and `test_gen` are instances of
-        :class:`quapy.data._lequa2022.SamplesFromDir`, a subclass of :class:`quapy.protocol.AbstractProtocol`,
+        :class:`quapy.data._lequa.SamplesFromDir`, a subclass of :class:`quapy.protocol.AbstractProtocol`,
         that return a series of samples stored in a directory which are labelled by prevalence.
     """
 
@@ -839,7 +835,9 @@ def fetch_lequa2022(task, data_home=None):
         tmp_path = join(lequa_dir, task + '_tmp.zip')
         download_file_if_not_exists(url, tmp_path)
         with zipfile.ZipFile(tmp_path) as file:
+            print(f'Unzipping {tmp_path}...', end='')
             file.extractall(unzipped_path)
+            print(f'[done]')
         os.remove(tmp_path)
 
     if not os.path.exists(join(lequa_dir, task)):
@@ -867,6 +865,35 @@ def fetch_lequa2022(task, data_home=None):
 
 
 def fetch_lequa2024(task, data_home=None, merge_T3=False):
+    """
+    Loads the official datasets provided for the `LeQua 2024 <https://lequa2024.github.io/index>`_ competition.
+    LeQua 2024 defines four tasks (T1, T2, T3, T4) related to the problem of quantification;
+    all tasks are affected by some type of dataset shift. Tasks T1 and T2 are akin to tasks T1A and T1B of LeQua 2022,
+    while T3 and T4 are new tasks introduced in LeQua 2024.
+
+    - Task T1 evaluates binary quantifiers under prior probability shift (akin to T1A of LeQua 2022).
+    - Task T2 evaluates single-label multi-class quantifiers (for n > 2 classes) under prior probability shift (akin to T1B of LeQua 2022).
+    - Task T3 evaluates ordinal quantifiers, where the classes are totally ordered.
+    - Task T4 also evaluates binary quantifiers, but under some mix of covariate shift and prior probability shift.
+
+    For a broader discussion, we refer to the `online official documentation <https://lequa2024.github.io/tasks/>`_
+
+    The datasets are downloaded only once, and stored locally for future reuse.
+
+    See `4b.lequa2024_experiments.py` provided in the example folder, which can serve as a guide on how to use these
+    datasets.
+
+    :param task: a string representing the task name; valid ones are T1, T2, T3, and T4
+    :param data_home: specify the quapy home directory where collections will be dumped (leave empty to use the default
+        ~/quapy_data/ directory)
+    :param merge_T3: bool, if False (default), returns a generator of training collections, corresponding to natural
+        groups of reviews; if True, returns one single :class:`quapy.data.base.LabelledCollection` representing the
+        entire training set, as a concatenation of all the training collections
+    :return: a tuple `(train, val_gen, test_gen)` where `train` is an instance of
+        :class:`quapy.data.base.LabelledCollection`, `val_gen` and `test_gen` are instances of
+        :class:`quapy.data._lequa.SamplesFromDir`, a subclass of :class:`quapy.protocol.AbstractProtocol`,
+        that return a series of samples stored in a directory which are labelled by prevalence.
+    """
 
     from quapy.data._lequa import load_vector_documents_2024, SamplesFromDir, LabelledCollectionsFromDir
 
@@ -909,11 +936,7 @@ def fetch_lequa2024(task, data_home=None, merge_T3=False):
     test_true_prev_path = join(lequa_dir, task, 'public', 'test_prevalences.txt')
     test_gen = SamplesFromDir(test_samples_path, test_true_prev_path, load_fn=load_fn)
 
-    if task != 'T3':
+    if task == 'T3':
-        tr_path = join(lequa_dir, task, 'public', 'training_data.txt')
-        train = LabelledCollection.load(tr_path, loader_func=load_fn)
-        return train, val_gen, test_gen
-    else:
         training_samples_path = join(lequa_dir, task, 'public', 'training_samples')
         training_true_prev_path = join(lequa_dir, task, 'public', 'training_prevalences.txt')
         train_gen = LabelledCollectionsFromDir(training_samples_path, training_true_prev_path, load_fn=load_fn)
@@ -922,7 +945,10 @@ def fetch_lequa2024(task, data_home=None, merge_T3=False):
             return train, val_gen, test_gen
         else:
             return train_gen, val_gen, test_gen
-
+    else:
+        tr_path = join(lequa_dir, task, 'public', 'training_data.txt')
+        train = LabelledCollection.load(tr_path, loader_func=load_fn)
+        return train, val_gen, test_gen
 
 
 def fetch_IFCB(single_sample_train=True, for_model_selection=False, data_home=None):

quapy/error.py

@@ -45,89 +45,95 @@ def acce(y_true, y_pred):
     return 1. - (y_true == y_pred).mean()
 
 
-def mae(prevs, prevs_hat):
+def mae(prevs_true, prevs_hat):
     """Computes the mean absolute error (see :meth:`quapy.error.ae`) across the sample pairs.
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
         prevalence values
     :return: mean absolute error
     """
-    return ae(prevs, prevs_hat).mean()
+    return ae(prevs_true, prevs_hat).mean()
 
 
-def ae(prevs, prevs_hat):
+def ae(prevs_true, prevs_hat):
     """Computes the absolute error between the two prevalence vectors.
      Absolute error between two prevalence vectors :math:`p` and :math:`\\hat{p}`  is computed as
      :math:`AE(p,\\hat{p})=\\frac{1}{|\\mathcal{Y}|}\\sum_{y\\in \\mathcal{Y}}|\\hat{p}(y)-p(y)|`,
      where :math:`\\mathcal{Y}` are the classes of interest.
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :return: absolute error
     """
-    assert prevs.shape == prevs_hat.shape, f'wrong shape {prevs.shape} vs. {prevs_hat.shape}'
+    prevs_true = np.asarray(prevs_true)
-    return abs(prevs_hat - prevs).mean(axis=-1)
+    prevs_hat = np.asarray(prevs_hat)
+    assert prevs_true.shape == prevs_hat.shape, f'wrong shape {prevs_true.shape} vs. {prevs_hat.shape}'
+    return abs(prevs_hat - prevs_true).mean(axis=-1)
 
 
-def nae(prevs, prevs_hat):
+def nae(prevs_true, prevs_hat):
     """Computes the normalized absolute error between the two prevalence vectors.
      Normalized absolute error between two prevalence vectors :math:`p` and :math:`\\hat{p}`  is computed as
      :math:`NAE(p,\\hat{p})=\\frac{AE(p,\\hat{p})}{z_{AE}}`,
      where :math:`z_{AE}=\\frac{2(1-\\min_{y\\in \\mathcal{Y}} p(y))}{|\\mathcal{Y}|}`, and :math:`\\mathcal{Y}`
      are the classes of interest.
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :return: normalized absolute error
     """
-    assert prevs.shape == prevs_hat.shape, f'wrong shape {prevs.shape} vs. {prevs_hat.shape}'
+    prevs_true = np.asarray(prevs_true)
-    return abs(prevs_hat - prevs).sum(axis=-1)/(2*(1-prevs.min(axis=-1)))
+    prevs_hat = np.asarray(prevs_hat)
+    assert prevs_true.shape == prevs_hat.shape, f'wrong shape {prevs_true.shape} vs. {prevs_hat.shape}'
+    return abs(prevs_hat - prevs_true).sum(axis=-1)/(2 * (1 - prevs_true.min(axis=-1)))
 
 
-def mnae(prevs, prevs_hat):
+def mnae(prevs_true, prevs_hat):
     """Computes the mean normalized absolute error (see :meth:`quapy.error.nae`) across the sample pairs.
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
         prevalence values
     :return: mean normalized absolute error
     """
-    return nae(prevs, prevs_hat).mean()
+    return nae(prevs_true, prevs_hat).mean()
 
 
-def mse(prevs, prevs_hat):
+def mse(prevs_true, prevs_hat):
     """Computes the mean squared error (see :meth:`quapy.error.se`) across the sample pairs.
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the
         true prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the
         predicted prevalence values
     :return: mean squared error
     """
-    return se(prevs, prevs_hat).mean()
+    return se(prevs_true, prevs_hat).mean()
 
 
-def se(prevs, prevs_hat):
+def se(prevs_true, prevs_hat):
     """Computes the squared error between the two prevalence vectors.
      Squared error between two prevalence vectors :math:`p` and :math:`\\hat{p}`  is computed as
      :math:`SE(p,\\hat{p})=\\frac{1}{|\\mathcal{Y}|}\\sum_{y\\in \\mathcal{Y}}(\\hat{p}(y)-p(y))^2`,
      where
      :math:`\\mathcal{Y}` are the classes of interest.
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :return: absolute error
     """
-    return ((prevs_hat - prevs) ** 2).mean(axis=-1)
+    prevs_true = np.asarray(prevs_true)
+    prevs_hat = np.asarray(prevs_hat)
+    return ((prevs_hat - prevs_true) ** 2).mean(axis=-1)
 
 
-def mkld(prevs, prevs_hat, eps=None):
+def mkld(prevs_true, prevs_hat, eps=None):
     """Computes the mean Kullback-Leibler divergence (see :meth:`quapy.error.kld`) across the
     sample pairs. The distributions are smoothed using the `eps` factor
     (see :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the true
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true
         prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
         prevalence values
@@ -137,10 +143,10 @@ def mkld(prevs, prevs_hat, eps=None):
         (which has thus to be set beforehand).
     :return: mean Kullback-Leibler distribution
     """
-    return kld(prevs, prevs_hat, eps).mean()
+    return kld(prevs_true, prevs_hat, eps).mean()
 
 
-def kld(prevs, prevs_hat, eps=None):
+def kld(prevs_true, prevs_hat, eps=None):
     """Computes the Kullback-Leibler divergence between the two prevalence distributions.
      Kullback-Leibler divergence between two prevalence distributions :math:`p` and :math:`\\hat{p}`
      is computed as
@@ -149,7 +155,7 @@ def kld(prevs, prevs_hat, eps=None):
      where :math:`\\mathcal{Y}` are the classes of interest.
      The distributions are smoothed using the `eps` factor (see :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :param eps: smoothing factor. KLD is not defined in cases in which the distributions contain
         zeros; `eps` is typically set to be :math:`\\frac{1}{2T}`, with :math:`T` the sample size.
@@ -158,17 +164,17 @@ def kld(prevs, prevs_hat, eps=None):
     :return: Kullback-Leibler divergence between the two distributions
     """
     eps = __check_eps(eps)
-    smooth_prevs = smooth(prevs, eps)
+    smooth_prevs = smooth(prevs_true, eps)
     smooth_prevs_hat = smooth(prevs_hat, eps)
     return (smooth_prevs*np.log(smooth_prevs/smooth_prevs_hat)).sum(axis=-1)
 
 
-def mnkld(prevs, prevs_hat, eps=None):
+def mnkld(prevs_true, prevs_hat, eps=None):
     """Computes the mean Normalized Kullback-Leibler divergence (see :meth:`quapy.error.nkld`)
     across the sample pairs. The distributions are smoothed using the `eps` factor
     (see :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
         prevalence values
     :param eps: smoothing factor. NKLD is not defined in cases in which the distributions contain
@@ -177,10 +183,10 @@ def mnkld(prevs, prevs_hat, eps=None):
         (which has thus to be set beforehand).
     :return: mean Normalized Kullback-Leibler distribution
     """
-    return nkld(prevs, prevs_hat, eps).mean()
+    return nkld(prevs_true, prevs_hat, eps).mean()
 
 
-def nkld(prevs, prevs_hat, eps=None):
+def nkld(prevs_true, prevs_hat, eps=None):
     """Computes the Normalized Kullback-Leibler divergence between the two prevalence distributions.
      Normalized Kullback-Leibler divergence between two prevalence distributions :math:`p` and
      :math:`\\hat{p}` is computed as
@@ -189,7 +195,7 @@ def nkld(prevs, prevs_hat, eps=None):
      :math:`\\mathcal{Y}` are the classes of interest.
      The distributions are smoothed using the `eps` factor (see :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :param eps: smoothing factor. NKLD is not defined in cases in which the distributions
         contain zeros; `eps` is typically set to be :math:`\\frac{1}{2T}`, with :math:`T` the sample
@@ -197,16 +203,16 @@ def nkld(prevs, prevs_hat, eps=None):
         `SAMPLE_SIZE` (which has thus to be set beforehand).
     :return: Normalized Kullback-Leibler divergence between the two distributions
     """
-    ekld = np.exp(kld(prevs, prevs_hat, eps))
+    ekld = np.exp(kld(prevs_true, prevs_hat, eps))
     return 2. * ekld / (1 + ekld) - 1.
 
 
-def mrae(prevs, prevs_hat, eps=None):
+def mrae(prevs_true, prevs_hat, eps=None):
     """Computes the mean relative absolute error (see :meth:`quapy.error.rae`) across
     the sample pairs. The distributions are smoothed using the `eps` factor (see
     :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the true
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true
         prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
         prevalence values
@@ -216,10 +222,10 @@ def mrae(prevs, prevs_hat, eps=None):
         the environment variable `SAMPLE_SIZE` (which has thus to be set beforehand).
     :return: mean relative absolute error
     """
-    return rae(prevs, prevs_hat, eps).mean()
+    return rae(prevs_true, prevs_hat, eps).mean()
 
 
-def rae(prevs, prevs_hat, eps=None):
+def rae(prevs_true, prevs_hat, eps=None):
     """Computes the absolute relative error between the two prevalence vectors.
      Relative absolute error between two prevalence vectors :math:`p` and :math:`\\hat{p}`
      is computed as
@@ -228,7 +234,7 @@ def rae(prevs, prevs_hat, eps=None):
      where :math:`\\mathcal{Y}` are the classes of interest.
      The distributions are smoothed using the `eps` factor (see :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :param eps: smoothing factor. `rae` is not defined in cases in which the true distribution
         contains zeros; `eps` is typically set to be :math:`\\frac{1}{2T}`, with :math:`T` the
@@ -237,12 +243,12 @@ def rae(prevs, prevs_hat, eps=None):
     :return: relative absolute error
     """
     eps = __check_eps(eps)
-    prevs = smooth(prevs, eps)
+    prevs_true = smooth(prevs_true, eps)
     prevs_hat = smooth(prevs_hat, eps)
-    return (abs(prevs - prevs_hat) / prevs).mean(axis=-1)
+    return (abs(prevs_true - prevs_hat) / prevs_true).mean(axis=-1)
 
 
-def nrae(prevs, prevs_hat, eps=None):
+def nrae(prevs_true, prevs_hat, eps=None):
     """Computes the normalized absolute relative error between the two prevalence vectors.
      Relative absolute error between two prevalence vectors :math:`p` and :math:`\\hat{p}`
      is computed as
@@ -252,7 +258,7 @@ def nrae(prevs, prevs_hat, eps=None):
      and :math:`\\mathcal{Y}` are the classes of interest.
      The distributions are smoothed using the `eps` factor (see :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
     :param eps: smoothing factor. `nrae` is not defined in cases in which the true distribution
         contains zeros; `eps` is typically set to be :math:`\\frac{1}{2T}`, with :math:`T` the
@@ -261,18 +267,18 @@ def nrae(prevs, prevs_hat, eps=None):
     :return: normalized relative absolute error
     """
     eps = __check_eps(eps)
-    prevs = smooth(prevs, eps)
+    prevs_true = smooth(prevs_true, eps)
     prevs_hat = smooth(prevs_hat, eps)
-    min_p = prevs.min(axis=-1)
+    min_p = prevs_true.min(axis=-1)
-    return (abs(prevs - prevs_hat) / prevs).sum(axis=-1)/(prevs.shape[-1]-1+(1-min_p)/min_p)
+    return (abs(prevs_true - prevs_hat) / prevs_true).sum(axis=-1)/(prevs_true.shape[-1] - 1 + (1 - min_p) / min_p)
 
 
-def mnrae(prevs, prevs_hat, eps=None):
+def mnrae(prevs_true, prevs_hat, eps=None):
     """Computes the mean normalized relative absolute error (see :meth:`quapy.error.nrae`) across
     the sample pairs. The distributions are smoothed using the `eps` factor (see
     :meth:`quapy.error.smooth`).
 
-    :param prevs: array-like of shape `(n_samples, n_classes,)` with the true
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true
         prevalence values
     :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
         prevalence values
@@ -282,32 +288,61 @@ def mnrae(prevs, prevs_hat, eps=None):
         the environment variable `SAMPLE_SIZE` (which has thus to be set beforehand).
     :return: mean normalized relative absolute error
     """
-    return nrae(prevs, prevs_hat, eps).mean()
+    return nrae(prevs_true, prevs_hat, eps).mean()
 
 
-def nmd(prevs, prevs_hat):
+def nmd(prevs_true, prevs_hat):
     """
     Computes the Normalized Match Distance; which is the Normalized Distance multiplied by the factor
     `1/(n-1)` to guarantee the measure ranges between 0 (best prediction) and 1 (worst prediction).
 
-    :param prevs: array-like of shape `(n_classes,)` or `(n_instances, n_classes)`  with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` or `(n_instances, n_classes)`  with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` or `(n_instances, n_classes)` with the predicted prevalence values
     :return: float in [0,1]
     """
-    n = prevs.shape[-1]
+    prevs_true = np.asarray(prevs_true)
-    return (1./(n-1))*np.mean(match_distance(prevs, prevs_hat))
+    prevs_hat = np.asarray(prevs_hat)
+    n = prevs_true.shape[-1]
+    return (1./(n-1))*np.mean(match_distance(prevs_true, prevs_hat))
 
 
-def md(prevs, prevs_hat, ERROR_TOL=1E-3):
+def bias_binary(prevs_true, prevs_hat):
+    """
+    Computes the (positive) bias in a binary problem. The bias is simply the difference between the
+    predicted positive value and the true positive value, so that a positive such value indicates the
+    prediction has positive bias (i.e., it tends to overestimate) the true value, and negative otherwise.
+    :math:`bias(p,\\hat{p})=\\hat{p}_1-p_1`,
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
+    :param prevs_hat: array-like of shape `(n_samples, n_classes,)` with the predicted
+        prevalence values
+    :return: binary bias
+    """
+    prevs_true = np.asarray(prevs_true)
+    prevs_hat = np.asarray(prevs_hat)
+    assert prevs_true.shape[-1] == 2 and prevs_true.shape[-1] == 2, f'bias_binary can only be applied to binary problems'
+    return prevs_hat[...,1]-prevs_true[...,1]
+
+
+def mean_bias_binary(prevs_true, prevs_hat):
+    """
+    Computes the mean of the (positive) bias in a binary problem.
+    :param prevs_true: array-like of shape `(n_classes,)` with the true prevalence values
+    :param prevs_hat: array-like of shape `(n_classes,)` with the predicted prevalence values
+    :return: mean binary bias
+    """
+    return np.mean(bias_binary(prevs_true, prevs_hat))
+
+
+def md(prevs_true, prevs_hat, ERROR_TOL=1E-3):
     """
     Computes the Match Distance, under the assumption that the cost in mistaking class i with class i+1 is 1 in
     all cases.
 
-    :param prevs: array-like of shape `(n_classes,)` or `(n_instances, n_classes)`  with the true prevalence values
+    :param prevs_true: array-like of shape `(n_classes,)` or `(n_instances, n_classes)`  with the true prevalence values
     :param prevs_hat: array-like of shape `(n_classes,)` or `(n_instances, n_classes)` with the predicted prevalence values
     :return: float
     """
-    P = np.cumsum(prevs, axis=-1)
+    P = np.cumsum(prevs_true, axis=-1)
     P_hat = np.cumsum(prevs_hat, axis=-1)
     assert np.all(np.isclose(P_hat[..., -1], 1.0, rtol=ERROR_TOL)), \
         'arg error in match_distance: the array does not represent a valid distribution'
@@ -324,6 +359,7 @@ def smooth(prevs, eps):
     :param eps: smoothing factor
     :return: array-like of shape `(n_classes,)` with the smoothed distribution
     """
+    prevs = np.asarray(prevs)
     n_classes = prevs.shape[-1]
     return (prevs + eps) / (eps * n_classes + 1)
 

quapy/evaluation.py

@@ -63,7 +63,7 @@ def prediction(
         protocol_with_predictions = protocol.on_preclassified_instances(pre_classified)
         return __prediction_helper(model.aggregate, protocol_with_predictions, verbose)
     else:
-        return __prediction_helper(model.quantify, protocol, verbose)
+        return __prediction_helper(model.predict, protocol, verbose)
 
 
 def __prediction_helper(quantification_fn, protocol: AbstractProtocol, verbose=False):

quapy/functional.py

@@ -7,6 +7,29 @@ import scipy
 import numpy as np
 
 
+# ------------------------------------------------------------------------------------------
+# General utils
+# ------------------------------------------------------------------------------------------
+
+def classes_from_labels(labels):
+    """
+    Obtains a np.ndarray with the (sorted) classes
+    :param labels: array-like with the instances' labels
+    :return: a sorted np.ndarray with the class labels
+    """
+    classes = np.unique(labels)
+    classes.sort()
+    return classes
+
+
+def num_classes_from_labels(labels):
+    """
+    Obtains the number of classes from an array-like of instance's labels
+    :param labels: array-like with the instances' labels
+    :return: int, the number of classes
+    """
+    return len(classes_from_labels(labels))
+
 # ------------------------------------------------------------------------------------------
 # Counter utils
 # ------------------------------------------------------------------------------------------
@@ -416,7 +439,7 @@ def argmin_prevalence(loss: Callable,
         raise NotImplementedError()
 
 
-def optim_minimize(loss: Callable, n_classes: int):
+def optim_minimize(loss: Callable, n_classes: int, 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
@@ -424,18 +447,24 @@ def optim_minimize(loss: Callable, n_classes: int):
 
     :param loss: (callable) the function to minimize
     :param n_classes: (int) the number of classes, i.e., the dimensionality of the prevalence vector
-    :return: (ndarray) the best prevalence vector found
+    :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 = np.full(fill_value=1 / n_classes, shape=(n_classes,))
+    uniform_distribution = 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)
-    return r.x
+    
+    if return_loss:
+        return r.x, r.fun
+    else:
+        return r.x
 
 
 def linear_search(loss: Callable, n_classes: int):

quapy/method/__init__.py

@@ -1,3 +1,8 @@
+import warnings
+from sklearn.exceptions import ConvergenceWarning
+warnings.simplefilter("ignore", ConvergenceWarning)
+
+from . import confidence
 from . import base
 from . import aggregative
 from . import non_aggregative
@@ -22,7 +27,8 @@ AGGREGATIVE_METHODS = {
     aggregative.KDEyML,
     aggregative.KDEyCS,
     aggregative.KDEyHD,
-    aggregative.BayesianCC
+    # aggregative.OneVsAllAggregative,
+    confidence.BayesianCC,
 }
 
 BINARY_METHODS = {
@@ -45,7 +51,7 @@ MULTICLASS_METHODS = {
     aggregative.KDEyML,
     aggregative.KDEyCS,
     aggregative.KDEyHD,
-    aggregative.BayesianCC
+    confidence.BayesianCC
 }
 
 NON_AGGREGATIVE_METHODS = {
@@ -62,3 +68,5 @@ QUANTIFICATION_METHODS = AGGREGATIVE_METHODS | NON_AGGREGATIVE_METHODS | META_ME
 
 
 
+
+

quapy/method/_kdey.py

@@ -1,10 +1,8 @@
-from typing import Union
 import numpy as np
 from sklearn.base import BaseEstimator
 from sklearn.neighbors import KernelDensity
 
 import quapy as qp
-from quapy.data import LabelledCollection
 from quapy.method.aggregative import AggregativeSoftQuantifier
 import quapy.functional as F
 
@@ -99,26 +97,27 @@ class KDEyML(AggregativeSoftQuantifier, KDEBase):
 
     which corresponds to the maximum likelihood estimate.
 
-    :param classifier: a sklearn's Estimator that generates a binary classifier.
+    :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 collection defining the specific set of data to use for validation.
+        for `k`); or as a tuple (X,y) defining the specific set of data to use for validation.
-        Alternatively, this set can be specified at fit time by indicating the exact set of data
-        on which the predictions are to be generated.
     :param bandwidth: float, the bandwidth of the Kernel
     :param random_state: a seed to be set before fitting any base quantifier (default None)
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=5, bandwidth=0.1, random_state=None):
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5, bandwidth=0.1,
-        self.classifier = qp._get_classifier(classifier)
+                 random_state=None):
-        self.val_split = val_split
+        super().__init__(classifier, fit_classifier, val_split)
         self.bandwidth = KDEBase._check_bandwidth(bandwidth)
         self.random_state=random_state
 
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
-        self.mix_densities = self.get_mixture_components(*classif_predictions.Xy, data.classes_, self.bandwidth)
+        self.mix_densities = self.get_mixture_components(classif_predictions, labels, self.classes_, self.bandwidth)
         return self
 
     def aggregate(self, posteriors: np.ndarray):
@@ -173,35 +172,35 @@ class KDEyHD(AggregativeSoftQuantifier, KDEBase):
     where the datapoints (trials) :math:`x_1,\\ldots,x_t\\sim_{\\mathrm{iid}} r` with :math:`r`  the
     uniform distribution.
 
-    :param classifier: a sklearn's Estimator that generates a binary classifier.
+    :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 collection defining the specific set of data to use for validation.
+        for `k`); or as a tuple (X,y) defining the specific set of data to use for validation.
-        Alternatively, this set can be specified at fit time by indicating the exact set of data
-        on which the predictions are to be generated.
     :param bandwidth: float, the bandwidth of the Kernel
     :param random_state: a seed to be set before fitting any base quantifier (default None)
     :param montecarlo_trials: number of Monte Carlo trials (default 10000)
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=5, divergence: str='HD',
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5, divergence: str='HD',
                  bandwidth=0.1, random_state=None, montecarlo_trials=10000):
 
-        self.classifier = qp._get_classifier(classifier)
+        super().__init__(classifier, fit_classifier, val_split)
-        self.val_split = val_split
         self.divergence = divergence
         self.bandwidth = KDEBase._check_bandwidth(bandwidth)
         self.random_state=random_state
         self.montecarlo_trials = montecarlo_trials
 
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
-        self.mix_densities = self.get_mixture_components(*classif_predictions.Xy, data.classes_, self.bandwidth)
+        self.mix_densities = self.get_mixture_components(classif_predictions, labels, self.classes_, self.bandwidth)
 
         N = self.montecarlo_trials
         rs = self.random_state
-        n = data.n_classes
+        n = len(self.classes_)
         self.reference_samples = np.vstack([kde_i.sample(N//n, random_state=rs) for kde_i in self.mix_densities])
         self.reference_classwise_densities = np.asarray([self.pdf(kde_j, self.reference_samples) for kde_j in self.mix_densities])
         self.reference_density = np.mean(self.reference_classwise_densities, axis=0)  # equiv. to (uniform @ self.reference_classwise_densities)
@@ -265,20 +264,20 @@ class KDEyCS(AggregativeSoftQuantifier):
 
     The authors showed that this distribution matching admits a closed-form solution
 
-    :param classifier: a sklearn's Estimator that generates a binary classifier.
+    :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 collection defining the specific set of data to use for validation.
+        for `k`); or as a tuple (X,y) defining the specific set of data to use for validation.
-        Alternatively, this set can be specified at fit time by indicating the exact set of data
-        on which the predictions are to be generated.
     :param bandwidth: float, the bandwidth of the Kernel
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=5, bandwidth=0.1):
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5, bandwidth=0.1):
-        self.classifier = qp._get_classifier(classifier)
+        super().__init__(classifier, fit_classifier, val_split)
-        self.val_split = val_split
         self.bandwidth = KDEBase._check_bandwidth(bandwidth)
 
     def gram_matrix_mix_sum(self, X, Y=None):
@@ -293,17 +292,17 @@ class KDEyCS(AggregativeSoftQuantifier):
         gram = norm_factor * rbf_kernel(X, Y, gamma=gamma)
         return gram.sum()
 
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
 
-        P, y = classif_predictions.Xy
+        P, y = classif_predictions, labels
-        n = data.n_classes
+        n = len(self.classes_)
 
         assert all(sorted(np.unique(y)) == np.arange(n)), \
             'label name gaps not allowed in current implementation'
 
         # counts_inv keeps track of the relative weight of each datapoint within its class
         # (i.e., the weight in its KDE model)
-        counts_inv = 1 / (data.counts())
+        counts_inv = 1 / (F.counts_from_labels(y, classes=self.classes_))
 
         # tr_tr_sums corresponds to symbol \overline{B} in the paper
         tr_tr_sums = np.zeros(shape=(n,n), dtype=float)

quapy/method/_neural.py

@@ -21,13 +21,13 @@ class QuaNetTrainer(BaseQuantifier):
     Example:
 
     >>> import quapy as qp
-    >>> from quapy.method_name.meta import QuaNet
+    >>> from quapy.method.meta import QuaNet
     >>> from quapy.classification.neural import NeuralClassifierTrainer, CNNnet
     >>>
     >>> # use samples of 100 elements
     >>> qp.environ['SAMPLE_SIZE'] = 100
     >>>
-    >>> # load the kindle dataset as text, and convert words to numerical indexes
+    >>> # load the Kindle dataset as text, and convert words to numerical indexes
     >>> dataset = qp.datasets.fetch_reviews('kindle', pickle=True)
     >>> qp.train.preprocessing.index(dataset, min_df=5, inplace=True)
     >>>
@@ -37,12 +37,14 @@ class QuaNetTrainer(BaseQuantifier):
     >>>
     >>> # train QuaNet (QuaNet is an alias to QuaNetTrainer)
     >>> model = QuaNet(classifier, qp.environ['SAMPLE_SIZE'], device='cuda')
-    >>> model.fit(dataset.training)
+    >>> model.fit(*dataset.training.Xy)
-    >>> estim_prevalence = model.quantify(dataset.test.instances)
+    >>> estim_prevalence = model.predict(dataset.test.instances)
 
     :param classifier: an object implementing `fit` (i.e., that can be trained on labelled data),
         `predict_proba` (i.e., that can generate posterior probabilities of unlabelled examples) and
         `transform` (i.e., that can generate embedded representations of the unlabelled instances).
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
+        learner has been trained outside the quantifier.
     :param sample_size: integer, the sample size; default is None, meaning that the sample size should be
         taken from qp.environ["SAMPLE_SIZE"]
     :param n_epochs: integer, maximum number of training epochs
@@ -64,6 +66,7 @@ class QuaNetTrainer(BaseQuantifier):
 
     def __init__(self,
                  classifier,
+                 fit_classifier=True,
                  sample_size=None,
                  n_epochs=100,
                  tr_iter_per_poch=500,
@@ -86,6 +89,7 @@ class QuaNetTrainer(BaseQuantifier):
             f'the classifier {classifier.__class__.__name__} does not seem to be able to produce posterior probabilities ' \
                 f'since it does not implement the method "predict_proba"'
         self.classifier = classifier
+        self.fit_classifier = fit_classifier
         self.sample_size = qp._get_sample_size(sample_size)
         self.n_epochs = n_epochs
         self.tr_iter = tr_iter_per_poch
@@ -111,20 +115,21 @@ class QuaNetTrainer(BaseQuantifier):
         self.__check_params_colision(self.quanet_params, self.classifier.get_params())
         self._classes_ = None
 
-    def fit(self, data: LabelledCollection, fit_classifier=True):
+    def fit(self, X, y):
         """
         Trains QuaNet.
 
-        :param data: the training data on which to train QuaNet. If `fit_classifier=True`, the data will be split in
+        :param X: the training instances on which to train QuaNet. If `fit_classifier=True`, the data will be split in
             40/40/20 for training the classifier, training QuaNet, and validating QuaNet, respectively. If
             `fit_classifier=False`, the data will be split in 66/34 for training QuaNet and validating it, respectively.
-        :param fit_classifier: if True, trains the classifier on a split containing 40% of the data
+        :param y: the labels of X
         :return: self
         """
+        data = LabelledCollection(X, y)
         self._classes_ = data.classes_
         os.makedirs(self.checkpointdir, exist_ok=True)
 
-        if fit_classifier:
+        if self.fit_classifier:
             classifier_data, unused_data = data.split_stratified(0.4)
             train_data, valid_data = unused_data.split_stratified(0.66)  # 0.66 split of 60% makes 40% and 20%
             self.classifier.fit(*classifier_data.Xy)
@@ -144,13 +149,13 @@ class QuaNetTrainer(BaseQuantifier):
         train_data_embed = LabelledCollection(self.classifier.transform(train_data.instances), train_data.labels, self._classes_)
 
         self.quantifiers = {
-            'cc': CC(self.classifier).fit(None, fit_classifier=False),
+            'cc': CC(self.classifier, fit_classifier=False).fit(*valid_data.Xy),
-            'acc': ACC(self.classifier).fit(None, fit_classifier=False, val_split=valid_data),
+            'acc': ACC(self.classifier, fit_classifier=False).fit(*valid_data.Xy),
-            'pcc': PCC(self.classifier).fit(None, fit_classifier=False),
+            'pcc': PCC(self.classifier, fit_classifier=False).fit(*valid_data.Xy),
-            'pacc': PACC(self.classifier).fit(None, fit_classifier=False, val_split=valid_data),
+            'pacc': PACC(self.classifier, fit_classifier=False).fit(*valid_data.Xy),
         }
         if classifier_data is not None:
-            self.quantifiers['emq'] = EMQ(self.classifier).fit(classifier_data, fit_classifier=False)
+            self.quantifiers['emq'] = EMQ(self.classifier, fit_classifier=False).fit(*valid_data.Xy)
 
         self.status = {
             'tr-loss': -1,
@@ -201,9 +206,9 @@ class QuaNetTrainer(BaseQuantifier):
 
         return prevs_estim
 
-    def quantify(self, instances):
+    def predict(self, X):
-        posteriors = self.classifier.predict_proba(instances)
+        posteriors = self.classifier.predict_proba(X)
-        embeddings = self.classifier.transform(instances)
+        embeddings = self.classifier.transform(X)
         quant_estims = self._get_aggregative_estims(posteriors)
         self.quanet.eval()
         with torch.no_grad():

quapy/method/_threshold_optim.py

@@ -18,18 +18,23 @@ class ThresholdOptimization(BinaryAggregativeQuantifier):
     that would allow for more true positives and many more false positives, on the grounds this
     would deliver larger denominators.
 
-    :param classifier: a sklearn's Estimator that generates a classifier
+    :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be
-    :param val_split: indicates the proportion of data to be used as a stratified held-out validation set in which the
+        the one indicated in `qp.environ['DEFAULT_CLS']`
-        misclassification rates are to be estimated.
+
-        This parameter can be indicated as a real value (between 0 and 1), representing a proportion of
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
-        validation data, or as an integer, indicating that the misclassification rates should be estimated via
+        learner has been trained outside the quantifier.
-        `k`-fold cross validation (this integer stands for the number of folds `k`, defaults 5), or as a
+
-        :class:`quapy.data.base.LabelledCollection` (the split itself).
+    :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 n_jobs: number of parallel workers
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=None, n_jobs=None):
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=None, n_jobs=None):
-        self.classifier = qp._get_classifier(classifier)
+        super().__init__(classifier, fit_classifier, val_split)
-        self.val_split = val_split
         self.n_jobs = qp._get_njobs(n_jobs)
 
     @abstractmethod
@@ -115,8 +120,8 @@ class ThresholdOptimization(BinaryAggregativeQuantifier):
             return 0
         return FP / (FP + TN)
 
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
-        decision_scores, y = classif_predictions.Xy
+        decision_scores, y = classif_predictions, labels
         # the standard behavior is to keep the best threshold only
         self.tpr, self.fpr, self.threshold = self._eval_candidate_thresholds(decision_scores, y)[0]
         return self
@@ -134,17 +139,22 @@ class T50(ThresholdOptimization):
     for the threshold that makes `tpr` closest to 0.5.
     The goal is to bring improved stability to the denominator of the adjustment.
 
-    :param classifier: a sklearn's Estimator that generates a classifier
+    :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be
-    :param val_split: indicates the proportion of data to be used as a stratified held-out validation set in which the
+        the one indicated in `qp.environ['DEFAULT_CLS']`
-        misclassification rates are to be estimated.
+
-        This parameter can be indicated as a real value (between 0 and 1), representing a proportion of
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
-        validation data, or as an integer, indicating that the misclassification rates should be estimated via
+        learner has been trained outside the quantifier.
-        `k`-fold cross validation (this integer stands for the number of folds `k`, defaults 5), or as a
+
-        :class:`quapy.data.base.LabelledCollection` (the split itself).
+    :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.
+
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=5):
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5):
-        super().__init__(classifier, val_split)
+        super().__init__(classifier, fit_classifier, val_split)
 
     def condition(self, tpr, fpr) -> float:
         return abs(tpr - 0.5)
@@ -158,17 +168,20 @@ class MAX(ThresholdOptimization):
     for the threshold that maximizes `tpr-fpr`.
     The goal is to bring improved stability to the denominator of the adjustment.
 
-    :param classifier: a sklearn's Estimator that generates a classifier
+    :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be
-    :param val_split: indicates the proportion of data to be used as a stratified held-out validation set in which the
+        the one indicated in `qp.environ['DEFAULT_CLS']`
-        misclassification rates are to be estimated.
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
-        This parameter can be indicated as a real value (between 0 and 1), representing a proportion of
+        learner has been trained outside the quantifier.
-        validation data, or as an integer, indicating that the misclassification rates should be estimated via
+    :param val_split: specifies the data used for generating classifier predictions. This specification
-        `k`-fold cross validation (this integer stands for the number of folds `k`, defaults 5), or as a
+        can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to
-        :class:`quapy.data.base.LabelledCollection` (the split itself).
+        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.
+
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=5):
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5):
-        super().__init__(classifier, val_split)
+        super().__init__(classifier, fit_classifier, val_split)
 
     def condition(self, tpr, fpr) -> float:
         # MAX strives to maximize (tpr - fpr), which is equivalent to minimize (fpr - tpr)
@@ -183,17 +196,20 @@ class X(ThresholdOptimization):
     for the threshold that yields `tpr=1-fpr`.
     The goal is to bring improved stability to the denominator of the adjustment.
 
-    :param classifier: a sklearn's Estimator that generates a classifier
+    :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be
-    :param val_split: indicates the proportion of data to be used as a stratified held-out validation set in which the
+        the one indicated in `qp.environ['DEFAULT_CLS']`
-        misclassification rates are to be estimated.
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
-        This parameter can be indicated as a real value (between 0 and 1), representing a proportion of
+        learner has been trained outside the quantifier.
-        validation data, or as an integer, indicating that the misclassification rates should be estimated via
+    :param val_split: specifies the data used for generating classifier predictions. This specification
-        `k`-fold cross validation (this integer stands for the number of folds `k`, defaults 5), or as a
+        can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to
-        :class:`quapy.data.base.LabelledCollection` (the split itself).
+        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.
+
     """
 
-    def __init__(self, classifier: BaseEstimator=None, val_split=5):
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5):
-        super().__init__(classifier, val_split)
+        super().__init__(classifier, fit_classifier, val_split)
 
     def condition(self, tpr, fpr) -> float:
         return abs(1 - (tpr + fpr))
@@ -207,22 +223,25 @@ class MS(ThresholdOptimization):
     class prevalence estimates for all decision thresholds and returns the median of them all.
     The goal is to bring improved stability to the denominator of the adjustment.
 
-    :param classifier: a sklearn's Estimator that generates a classifier
+    :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be
-    :param val_split: indicates the proportion of data to be used as a stratified held-out validation set in which the
+        the one indicated in `qp.environ['DEFAULT_CLS']`
-        misclassification rates are to be estimated.
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
-        This parameter can be indicated as a real value (between 0 and 1), representing a proportion of
+        learner has been trained outside the quantifier.
-        validation data, or as an integer, indicating that the misclassification rates should be estimated via
+    :param val_split: specifies the data used for generating classifier predictions. This specification
-        `k`-fold cross validation (this integer stands for the number of folds `k`, defaults 5), or as a
+        can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to
-        :class:`quapy.data.base.LabelledCollection` (the split itself).
+        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.
     """
-    def __init__(self, classifier: BaseEstimator=None, val_split=5):
+
-        super().__init__(classifier, val_split)
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5):
+        super().__init__(classifier, fit_classifier, val_split)
 
     def condition(self, tpr, fpr) -> float:
         return 1
 
-    def aggregation_fit(self, classif_predictions: LabelledCollection, data: LabelledCollection):
+    def aggregation_fit(self, classif_predictions, labels):
-        decision_scores, y = classif_predictions.Xy
+        decision_scores, y = classif_predictions, labels
         # keeps all candidates
         tprs_fprs_thresholds = self._eval_candidate_thresholds(decision_scores, y)
         self.tprs = tprs_fprs_thresholds[:, 0]
@@ -246,16 +265,19 @@ class MS2(MS):
     which `tpr-fpr>0.25`
     The goal is to bring improved stability to the denominator of the adjustment.
 
-    :param classifier: a sklearn's Estimator that generates a classifier
+    :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be
-    :param val_split: indicates the proportion of data to be used as a stratified held-out validation set in which the
+        the one indicated in `qp.environ['DEFAULT_CLS']`
-        misclassification rates are to be estimated.
+    :param fit_classifier: whether to train the learner (default is True). Set to False if the
-        This parameter can be indicated as a real value (between 0 and 1), representing a proportion of
+        learner has been trained outside the quantifier.
-        validation data, or as an integer, indicating that the misclassification rates should be estimated via
+    :param val_split: specifies the data used for generating classifier predictions. This specification
-        `k`-fold cross validation (this integer stands for the number of folds `k`, defaults 5), or as a
+        can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to
-        :class:`quapy.data.base.LabelledCollection` (the split itself).
+        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.
     """
-    def __init__(self, classifier: BaseEstimator=None, val_split=5):
+
-        super().__init__(classifier, val_split)
+    def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5):
+        super().__init__(classifier, fit_classifier, val_split)
 
     def discard(self, tpr, fpr) -> bool:
         return (tpr-fpr) <= 0.25

quapy/method/base.py

@@ -14,30 +14,40 @@ import numpy as np
 class BaseQuantifier(BaseEstimator):
     """
     Abstract Quantifier. A quantifier is defined as an object of a class that implements the method :meth:`fit` on
-    :class:`quapy.data.base.LabelledCollection`, the method :meth:`quantify`, and the :meth:`set_params` and
+    a pair X, y, the method :meth:`predict`, and the :meth:`set_params` and
     :meth:`get_params` for model selection (see :meth:`quapy.model_selection.GridSearchQ`)
     """
 
     @abstractmethod
-    def fit(self, data: LabelledCollection):
+    def fit(self, X, y):
         """
-        Trains a quantifier.
+        Generates a quantifier.
 
-        :param data: a :class:`quapy.data.base.LabelledCollection` consisting of the training data
+        :param X: array-like, the training instances
+        :param y: array-like, the labels
         :return: self
         """
         ...
 
     @abstractmethod
-    def quantify(self, instances):
+    def predict(self, X):
         """
         Generate class prevalence estimates for the sample's instances
 
-        :param instances: array-like
+        :param X: array-like, the test instances
         :return: `np.ndarray` of shape `(n_classes,)` with class prevalence estimates.
         """
         ...
 
+    def quantify(self, X):
+        """
+        Alias to :meth:`predict`, for old compatibility
+
+        :param X: array-like
+        :return: `np.ndarray` of shape `(n_classes,)` with class prevalence estimates.
+        """
+        return self.predict(X)
+
 
 class BinaryQuantifier(BaseQuantifier):
     """
@@ -45,8 +55,9 @@ class BinaryQuantifier(BaseQuantifier):
     (typically, to be interpreted as one class and its complement).
     """
 
-    def _check_binary(self, data: LabelledCollection, quantifier_name):
+    def _check_binary(self, y, quantifier_name):
-        assert data.binary, f'{quantifier_name} works only on problems of binary classification. ' \
+        n_classes = len(set(y))
+        assert n_classes==2, f'{quantifier_name} works only on problems of binary classification. ' \
                             f'Use the class OneVsAll to enable {quantifier_name} work on single-label data.'
 
 
@@ -66,7 +77,7 @@ def newOneVsAll(binary_quantifier: BaseQuantifier, n_jobs=None):
 class OneVsAllGeneric(OneVsAll, BaseQuantifier):
     """
     Allows any binary quantifier to perform quantification on single-label datasets. The method maintains one binary
-    quantifier for each class, and then l1-normalizes the outputs so that the class prevelence values sum up to 1.
+    quantifier for each class, and then l1-normalizes the outputs so that the class prevalence values sum up to 1.
     """
 
     def __init__(self, binary_quantifier: BaseQuantifier, n_jobs=None):
@@ -78,32 +89,32 @@ class OneVsAllGeneric(OneVsAll, BaseQuantifier):
         self.binary_quantifier = binary_quantifier
         self.n_jobs = qp._get_njobs(n_jobs)
 
-    def fit(self, data: LabelledCollection, fit_classifier=True):
+    def fit(self, X, y):
-        assert not data.binary, f'{self.__class__.__name__} expect non-binary data'
+        self.classes = sorted(np.unique(y))
-        assert fit_classifier == True, 'fit_classifier must be True'
+        assert len(self.classes)!=2, f'{self.__class__.__name__} expect non-binary data'
 
-        self.dict_binary_quantifiers = {c: deepcopy(self.binary_quantifier) for c in data.classes_}
+        self.dict_binary_quantifiers = {c: deepcopy(self.binary_quantifier) for c in self.classes}
-        self._parallel(self._delayed_binary_fit, data)
+        self._parallel(self._delayed_binary_fit, X, y)
         return self
 
     def _parallel(self, func, *args, **kwargs):
         return np.asarray(
             Parallel(n_jobs=self.n_jobs, backend='threading')(
-                delayed(func)(c, *args, **kwargs) for c in self.classes_
+                delayed(func)(c, *args, **kwargs) for c in self.classes
             )
         )
 
-    def quantify(self, instances):
+    def predict(self, X):
-        prevalences = self._parallel(self._delayed_binary_predict, instances)
+        prevalences = self._parallel(self._delayed_binary_predict, X)
         return qp.functional.normalize_prevalence(prevalences)
 
-    @property
+    # @property
-    def classes_(self):
+    # def classes_(self):
-        return sorted(self.dict_binary_quantifiers.keys())
+    #     return sorted(self.dict_binary_quantifiers.keys())
 
     def _delayed_binary_predict(self, c, X):
-        return self.dict_binary_quantifiers[c].quantify(X)[1]
+        return self.dict_binary_quantifiers[c].predict(X)[1]
 
-    def _delayed_binary_fit(self, c, data):
+    def _delayed_binary_fit(self, c, X, y):
-        bindata = LabelledCollection(data.instances, data.labels == c, classes=[False, True])
+        bindata = LabelledCollection(X, y == c, classes=[False, True])
-        self.dict_binary_quantifiers[c].fit(bindata)
+        self.dict_binary_quantifiers[c].fit(*bindata.Xy)

quapy/method/composable.py

@@ -1,19 +1,28 @@
 """This module allows the composition of quantification methods from loss functions and feature transformations. This functionality is realized through an integration of the qunfold package: https://github.com/mirkobunse/qunfold."""
 
-_import_error_message = """qunfold, the back-end of quapy.method.composable, is not properly installed.
+from dataclasses import dataclass
+from packaging.version import Version
+
+from .base import BaseQuantifier
+
+# what to display when an ImportError is thrown
+_IMPORT_ERROR_MESSAGE = """qunfold, the back-end of quapy.method.composable, is not properly installed.
 
 To fix this error, call:
 
     pip install --upgrade pip setuptools wheel
     pip install "jax[cpu]"
-    pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.4"
+    pip install "qunfold @ git+https://github.com/mirkobunse/qunfold@v0.1.5"
 """
 
+# try to import members of qunfold as members of this module
 try:
     import qunfold
-    from qunfold.quapy import QuaPyWrapper
+    from qunfold.base import BaseMixin
+    from qunfold.methods import AbstractMethod
     from qunfold.sklearn import CVClassifier
     from qunfold import (
+        LinearMethod, # methods
         LeastSquaresLoss, # losses
         BlobelLoss,
         EnergyLoss,
@@ -21,46 +30,95 @@ try:
         CombinedLoss,
         TikhonovRegularization,
         TikhonovRegularized,
-        ClassTransformer, # transformers
+        ClassRepresentation, # representations
-        HistogramTransformer,
+        HistogramRepresentation,
-        DistanceTransformer,
+        DistanceRepresentation,
-        KernelTransformer,
+        KernelRepresentation,
-        EnergyKernelTransformer,
+        EnergyKernelRepresentation,
-        LaplacianKernelTransformer,
+        LaplacianKernelRepresentation,
-        GaussianKernelTransformer,
+        GaussianKernelRepresentation,
-        GaussianRFFKernelTransformer,
+        GaussianRFFKernelRepresentation,
     )
-
-    __all__ = [ # control public members, e.g., for auto-documentation in sphinx; omit QuaPyWrapper
-        "ComposableQuantifier",
-        "CVClassifier",
-        "LeastSquaresLoss",
-        "BlobelLoss",
-        "EnergyLoss",
-        "HellingerSurrogateLoss",
-        "CombinedLoss",
-        "TikhonovRegularization",
-        "TikhonovRegularized",
-        "ClassTransformer",
-        "HistogramTransformer",
-        "DistanceTransformer",
-        "KernelTransformer",
-        "EnergyKernelTransformer",
-        "LaplacianKernelTransformer",
-        "GaussianKernelTransformer",
-        "GaussianRFFKernelTransformer",
-    ]
 except ImportError as e:
-    raise ImportError(_import_error_message) from e
+    raise ImportError(_IMPORT_ERROR_MESSAGE) from e
 
-def ComposableQuantifier(loss, transformer, **kwargs):
+__all__ = [ # control public members, e.g., for auto-documentation in sphinx
+    "QUnfoldWrapper",
+    "ComposableQuantifier",
+    "CVClassifier",
+    "LeastSquaresLoss",
+    "BlobelLoss",
+    "EnergyLoss",
+    "HellingerSurrogateLoss",
+    "CombinedLoss",
+    "TikhonovRegularization",
+    "TikhonovRegularized",
+    "ClassRepresentation",
+    "HistogramRepresentation",
+    "DistanceRepresentation",
+    "KernelRepresentation",
+    "EnergyKernelRepresentation",
+    "LaplacianKernelRepresentation",
+    "GaussianKernelRepresentation",
+    "GaussianRFFKernelRepresentation",
+]
+
+
+def check_compatible_qunfold_version():
+    try:
+        version_str = qunfold.__version__
+    except AttributeError:
+        # versions of qunfold <= 0.1.4 did not declare __version__ in the __init__.py but only in the setup.py
+        version_str = "0.1.4"
+
+    installed_ver = Version(version_str)
+    required_ver = Version("0.1.5")
+    compatible = installed_ver.base_version == required_ver.base_version or installed_ver>=required_ver
+    return compatible
+
+
+@dataclass
+class QUnfoldWrapper(BaseQuantifier,BaseMixin):
+    """A thin wrapper for using qunfold methods in QuaPy.
+
+    Args:
+      _method: An instance of `qunfold.methods.AbstractMethod` to wrap.
+
+    Examples:
+      Here, we wrap an instance of ACC to perform a grid search with QuaPy.
+
+        >>> from qunfold import ACC
+        >>> qunfold_method = QUnfoldWrapper(ACC(RandomForestClassifier(obb_score=True)))
+        >>> quapy.model_selection.GridSearchQ(
+        >>>     model = qunfold_method,
+        >>>     param_grid = { # try both splitting criteria
+        >>>         "representation__classifier__estimator__criterion": ["gini", "entropy"],
+        >>>     },
+        >>>     # ...
+        >>> )
+    """
+    _method: AbstractMethod
+    def fit(self, X, y): # data is a qp.LabelledCollection
+        self._method.fit(X, y)
+        return self
+    def predict(self, X):
+        return self._method.predict(X)
+    def set_params(self, **params):
+        self._method.set_params(**params)
+        return self
+    def get_params(self, deep=True):
+        return self._method.get_params(deep)
+    def __str__(self):
+        return self._method.__str__()
+
+def ComposableQuantifier(loss, representation, **kwargs):
     """A generic quantification / unfolding method that solves a linear system of equations.
 
     This class represents any quantifier that can be described in terms of a loss function, a feature transformation, and a regularization term. In this implementation, the loss is minimized through unconstrained second-order minimization. Valid probability estimates are ensured through a soft-max trick by Bunse (2022).
 
     Args:
         loss: An instance of a loss class from `quapy.methods.composable`.
-        transformer: An instance of a transformer class from `quapy.methods.composable`.
+        representation: An instance of a representation class from `quapy.methods.composable`.
         solver (optional): The `method` argument in `scipy.optimize.minimize`. Defaults to `"trust-ncg"`.
         solver_options (optional): The `options` argument in `scipy.optimize.minimize`. Defaults to `{"gtol": 1e-8, "maxiter": 1000}`.
         seed (optional): A random number generator seed from which a numpy RandomState is created. Defaults to `None`.
@@ -72,12 +130,12 @@ def ComposableQuantifier(loss, transformer, **kwargs):
             >>>     ComposableQuantifier,
             >>>     TikhonovRegularized,
             >>>     LeastSquaresLoss,
-            >>>     ClassTransformer,
+            >>>     ClassRepresentation,
             >>> )
             >>> from sklearn.ensemble import RandomForestClassifier
             >>> o_acc = ComposableQuantifier(
             >>>     TikhonovRegularized(LeastSquaresLoss(), 0.01),
-            >>>     ClassTransformer(RandomForestClassifier(oob_score=True))
+            >>>     ClassRepresentation(RandomForestClassifier(oob_score=True))
             >>> )
 
         Here, we perform hyper-parameter optimization with the ordinal ACC.
@@ -85,7 +143,7 @@ def ComposableQuantifier(loss, transformer, **kwargs):
             >>> quapy.model_selection.GridSearchQ(
             >>>     model = o_acc,
             >>>     param_grid = { # try both splitting criteria
-            >>>         "transformer__classifier__estimator__criterion": ["gini", "entropy"],
+            >>>         "representation__classifier__estimator__criterion": ["gini", "entropy"],
             >>>     },
             >>>     # ...
             >>> )
@@ -96,7 +154,7 @@ def ComposableQuantifier(loss, transformer, **kwargs):
             >>> from sklearn.linear_model import LogisticRegression
             >>> acc_lr = ComposableQuantifier(
             >>>     LeastSquaresLoss(),
-            >>>     ClassTransformer(CVClassifier(LogisticRegression(), 10))
+            >>>     ClassRepresentation(CVClassifier(LogisticRegression(), 10))
             >>> )
         """
-    return QuaPyWrapper(qunfold.GenericMethod(loss, transformer, **kwargs))
+    return QUnfoldWrapper(LinearMethod(loss, representation, **kwargs))

quapy/method/confidence.py

@@ -0,0 +1,551 @@
+import numpy as np
+from sklearn.base import BaseEstimator
+from sklearn.metrics import confusion_matrix
+
+import quapy as qp
+import quapy.functional as F
+from quapy.method import _bayesian
+from quapy.method.aggregative import AggregativeCrispQuantifier
+from quapy.data import LabelledCollection
+from quapy.method.aggregative import AggregativeQuantifier
+from scipy.stats import chi2
+from sklearn.utils import resample
+from abc import ABC, abstractmethod
+from scipy.special import softmax, factorial
+import copy
+from functools import lru_cache
+
+"""
+This module provides implementation of different types of confidence regions, and the implementation of Bootstrap
+for AggregativeQuantifiers.
+"""
+
+class ConfidenceRegionABC(ABC):
+    """
+    Abstract class of confidence regions
+    """
+
+    @abstractmethod
+    def point_estimate(self) -> np.ndarray:
+        """
+        Returns the point estimate corresponding to a set of bootstrap estimates.
+
+        :return: np.ndarray
+        """
+        ...
+
+    def ndim(self) -> int:
+        """
+        Number of dimensions of the region. This number corresponds to the total number of classes. The dimensionality
+        of the simplex is therefore ndim-1
+
+        :return: int
+        """
+        return len(self.point_estimate())
+
+    @abstractmethod
+    def coverage(self, true_value) -> float:
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        ...
+
+    @lru_cache
+    def simplex_portion(self):
+        """
+        Computes the fraction of the simplex which is covered by the region. This is not the volume of the region
+        itself (which could lie outside the boundaries of the simplex), but the actual fraction of the simplex
+        contained in the region. A default implementation, based on Monte Carlo approximation, is provided.
+
+        :return: float, the fraction of the simplex covered by the region
+        """
+        return self.montecarlo_proportion()
+
+    @lru_cache
+    def montecarlo_proportion(self, n_trials=10_000):
+        """
+        Estimates, via a Monte Carlo approach, the fraction of the simplex covered by the region. This is carried
+        out by returning the fraction of the `n_trials` points, uniformly drawn at random from the simplex, that
+        are included in the region. The value is only computed once when multiple calls are made.
+
+        :return: float in [0,1]
+        """
+        with qp.util.temp_seed(0):
+            uniform_simplex = F.uniform_simplex_sampling(n_classes=self.ndim(), size=n_trials)
+        proportion = np.clip(self.coverage(uniform_simplex), 0., 1.)
+        return proportion
+
+
+class WithConfidenceABC(ABC):
+    """
+    Abstract class for confidence regions.
+    """
+    METHODS = ['intervals', 'ellipse', 'ellipse-clr']
+
+    @abstractmethod
+    def quantify_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
+        """
+        Adds the method `quantify_conf` to the interface. This method returns not only the point-estimate, but
+        also the confidence region around it.
+
+        :param instances: a np.ndarray of shape (n_instances, n_features,)
+        :confidence_level: float in (0, 1)
+        :return: a tuple (`point_estimate`, `conf_region`), where `point_estimate` is a np.ndarray of shape
+            (n_classes,) and  `conf_region` is an object from :class:`ConfidenceRegionABC`
+        """
+        ...
+
+    @classmethod
+    def construct_region(cls, prev_estims, confidence_level=0.95, method='intervals'):
+        """
+        Construct a confidence region given many prevalence estimations.
+
+        :param prev_estims: np.ndarray of shape (n_estims, n_classes)
+        :param confidence_level: float, the confidence level for the region (default 0.95)
+        :param method: str, indicates the method for constructing regions. Set to `intervals` for
+            constructing confidence intervals (default), or to `ellipse` for constructing an
+            ellipse in the probability simplex, or to `ellipse-clr` for constructing an ellipse
+            in the Centered-Log Ratio (CLR) unconstrained space.
+        """
+        region = None
+        if method == 'intervals':
+            region = ConfidenceIntervals(prev_estims, confidence_level=confidence_level)
+        elif method == 'ellipse':
+            region = ConfidenceEllipseSimplex(prev_estims, confidence_level=confidence_level)
+        elif method == 'ellipse-clr':
+            region = ConfidenceEllipseCLR(prev_estims, confidence_level=confidence_level)
+
+        if region is None:
+            raise NotImplementedError(f'unknown method {method}')
+
+        return region
+
+def simplex_volume(n):
+    """
+    Computes the volume of the n-dimensional simplex. For n classes, the corresponding volume
+    is :meth:`simplex_volume(n-1)` since the simplex has one degree of freedom less.
+
+    :param n: int, the dimensionality of the simplex
+    :return: float, the volume of the n-dimensional simplex
+    """
+    return 1 / factorial(n)
+
+
+def within_ellipse_prop(values, mean, prec_matrix, chi2_critical):
+    """
+    Checks the proportion of values that belong to the ellipse with center `mean` and precision matrix `prec_matrix`
+    at a distance `chi2_critical`.
+
+    :param values: a np.ndarray of shape (n_dim,) or (n_values, n_dim,)
+    :param mean: a np.ndarray of shape (n_dim,) with the center of the ellipse
+    :param prec_matrix: a np.ndarray with the precision matrix (inverse of the
+        covariance matrix) of the ellipse. If this inverse cannot be computed
+        then None must be passed
+    :param chi2_critical: float, the chi2 critical value
+
+    :return: float in [0,1], the fraction of values that are contained in the ellipse
+        defined by the mean (center), the precision matrix (shape), and the chi2_critical value (distance).
+        If `values` is only one value, then either 0. (not contained) or 1. (contained) is returned.
+    """
+    if prec_matrix is None:
+        return 0.
+
+    diff = values - mean  # Mahalanobis distance
+
+    d_M_squared = diff @ prec_matrix @ diff.T  # d_M^2
+    if d_M_squared.ndim == 2:
+        d_M_squared = np.diag(d_M_squared)
+
+    within_elipse = (d_M_squared <= chi2_critical)
+
+    if isinstance(within_elipse, np.ndarray):
+        within_elipse = np.mean(within_elipse)
+
+    return within_elipse * 1.0
+
+
+class ConfidenceEllipseSimplex(ConfidenceRegionABC):
+    """
+    Instantiates a Confidence Ellipse in the probability simplex.
+
+    :param X: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    """
+
+    def __init__(self, X, confidence_level=0.95):
+
+        assert 0. < confidence_level < 1., f'{confidence_level=} must be in range(0,1)'
+
+        X = np.asarray(X)
+
+        self.mean_ = X.mean(axis=0)
+        self.cov_ = np.cov(X, rowvar=False, ddof=1)
+
+        try:
+            self.precision_matrix_ = np.linalg.inv(self.cov_)
+        except:
+            self.precision_matrix_ = None
+
+        self.dim = X.shape[-1]
+        self.ddof = self.dim - 1
+
+        # critical chi-square value
+        self.confidence_level = confidence_level
+        self.chi2_critical_ = chi2.ppf(confidence_level, df=self.ddof)
+
+    def point_estimate(self):
+        """
+        Returns the point estimate, the center of the ellipse.
+
+        :return: np.ndarray of shape (n_classes,)
+        """
+        return self.mean_
+
+    def coverage(self, true_value):
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        return within_ellipse_prop(true_value, self.mean_, self.precision_matrix_, self.chi2_critical_)
+
+
+class ConfidenceEllipseCLR(ConfidenceRegionABC):
+    """
+    Instantiates a Confidence Ellipse in the Centered-Log Ratio (CLR) space.
+
+    :param X: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    """
+
+    def __init__(self, X, confidence_level=0.95):
+        self.clr = CLRtransformation()
+        Z = self.clr(X)
+        self.mean_ = np.mean(X, axis=0)
+        self.conf_region_clr = ConfidenceEllipseSimplex(Z, confidence_level=confidence_level)
+
+    def point_estimate(self):
+        """
+        Returns the point estimate, the center of the ellipse.
+
+        :return: np.ndarray of shape (n_classes,)
+        """
+        # The inverse of the CLR does not coincide with the true mean, because the geometric mean
+        # requires smoothing the prevalence vectors and this affects the softmax (inverse);
+        # return self.clr.inverse(self.mean_) # <- does not coincide
+        return self.mean_
+
+    def coverage(self, true_value):
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        transformed_values = self.clr(true_value)
+        return self.conf_region_clr.coverage(transformed_values)
+
+
+class ConfidenceIntervals(ConfidenceRegionABC):
+    """
+    Instantiates a region based on (independent) Confidence Intervals.
+
+    :param X: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    """
+    def __init__(self, X, confidence_level=0.95):
+        assert 0 < confidence_level < 1, f'{confidence_level=} must be in range(0,1)'
+
+        X = np.asarray(X)
+
+        self.means_ = X.mean(axis=0)
+        alpha = 1-confidence_level
+        low_perc = (alpha/2.)*100
+        high_perc = (1-alpha/2.)*100
+        self.I_low, self.I_high = np.percentile(X, q=[low_perc, high_perc], axis=0)
+
+    def point_estimate(self):
+        """
+        Returns the point estimate, the class-wise average of the bootstrapped estimates
+
+        :return: np.ndarray of shape (n_classes,)
+        """
+        return self.means_
+
+    def coverage(self, true_value):
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        within_intervals = np.logical_and(self.I_low <= true_value, true_value <= self.I_high)
+        within_all_intervals = np.all(within_intervals, axis=-1, keepdims=True)
+        proportion = within_all_intervals.mean()
+
+        return proportion
+
+
+class CLRtransformation:
+    """
+    Centered log-ratio, from component analysis
+    """
+    def __call__(self, X, epsilon=1e-6):
+        """
+        Applies the CLR function to X thus mapping the instances, which are contained in `\\mathcal{R}^{n}` but
+        actually lie on a `\\mathcal{R}^{n-1}` simplex, onto an unrestricted space in :math:`\\mathcal{R}^{n}`
+
+        :param X: np.ndarray of (n_instances, n_dimensions) to be transformed
+        :param epsilon: small float for prevalence smoothing
+        :return: np.ndarray of (n_instances, n_dimensions), the CLR-transformed points
+        """
+        X = np.asarray(X)
+        X = qp.error.smooth(X, epsilon)
+        G = np.exp(np.mean(np.log(X), axis=-1, keepdims=True))  # geometric mean
+        return np.log(X / G)
+
+    def inverse(self, X):
+        """
+        Inverse function. However, clr.inverse(clr(X)) does not exactly coincide with X due to smoothing.
+
+        :param X: np.ndarray of (n_instances, n_dimensions) to be transformed
+        :return: np.ndarray of (n_instances, n_dimensions), the CLR-transformed points
+        """
+        return softmax(X, axis=-1)
+
+
+class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
+    """
+    Aggregative Bootstrap allows any AggregativeQuantifier to get confidence regions around
+    point-estimates of class prevalence values. This method implements some optimizations for
+    speeding up the computations, which are only possible due to the two phases of the aggregative
+    quantifiers.
+
+    During training, the bootstrap repetitions are only carried out over pre-classified training instances,
+    after the classifier has been trained (only once), in order to train a series of aggregation
+    functions (model-based approach).
+
+    During inference, the bootstrap repetitions are applied to the pre-classified test instances.
+
+    :param quantifier: an aggregative quantifier
+    :para n_train_samples: int, the number of training resamplings (defaults to 1, set to > 1 to activate a
+        model-based bootstrap approach)
+    :para n_test_samples: int, the number of test resamplings (defaults to 500, set to > 1 to activate a
+        population-based bootstrap approach)
+    :param confidence_level: float, the confidence level for the confidence region (default 0.95)
+    :param region: string, set to `intervals` for constructing confidence intervals (default), or to
+        `ellipse` for constructing an ellipse in the probability simplex, or to `ellipse-clr` for
+        constructing an ellipse in the Centered-Log Ratio (CLR) unconstrained space.
+    :param random_state: int for replicating samples, None (default) for non-replicable samples
+    """
+
+    def __init__(self,
+                 quantifier: AggregativeQuantifier,
+                 n_train_samples=1,
+                 n_test_samples=500,
+                 confidence_level=0.95,
+                 region='intervals',
+                 random_state=None):
+
+        assert isinstance(quantifier, AggregativeQuantifier), \
+            f'base quantifier does not seem to be an instance of {AggregativeQuantifier.__name__}'
+        assert n_train_samples >= 1, \
+            f'{n_train_samples=} must be >= 1'
+        assert n_test_samples >= 1, \
+            f'{n_test_samples=} must be >= 1'
+        assert n_test_samples>1 or n_train_samples>1, \
+            f'either {n_test_samples=} or {n_train_samples=} must be >1'
+
+        self.quantifier = quantifier
+        self.n_train_samples = n_train_samples
+        self.n_test_samples = n_test_samples
+        self.confidence_level = confidence_level
+        self.region = region
+        self.random_state = random_state
+
+    def aggregation_fit(self, classif_predictions, labels):
+        data = LabelledCollection(classif_predictions, labels, classes=self.classes_)
+        self.quantifiers = []
+        if self.n_train_samples==1:
+            self.quantifier.aggregation_fit(classif_predictions, labels)
+            self.quantifiers.append(self.quantifier)
+        else:
+            # model-based bootstrap (only on the aggregative part)
+            n_examples = len(data)
+            full_index = np.arange(n_examples)
+            with qp.util.temp_seed(self.random_state):
+                for i in range(self.n_train_samples):
+                    quantifier = copy.deepcopy(self.quantifier)
+                    index = resample(full_index, n_samples=n_examples)
+                    classif_predictions_i = classif_predictions.sampling_from_index(index)
+                    data_i = data.sampling_from_index(index)
+                    quantifier.aggregation_fit(classif_predictions_i, data_i)
+                    self.quantifiers.append(quantifier)
+        return self
+
+    def aggregate(self, classif_predictions: np.ndarray):
+        prev_mean, self.confidence = self.aggregate_conf(classif_predictions)
+        return prev_mean
+
+    def aggregate_conf(self, classif_predictions: np.ndarray, confidence_level=None):
+        if confidence_level is None:
+            confidence_level = self.confidence_level
+
+        n_samples = classif_predictions.shape[0]
+        prevs = []
+        with qp.util.temp_seed(self.random_state):
+            for quantifier in self.quantifiers:
+                for i in range(self.n_test_samples):
+                    sample_i = resample(classif_predictions, n_samples=n_samples)
+                    prev_i = quantifier.aggregate(sample_i)
+                    prevs.append(prev_i)
+
+        conf = WithConfidenceABC.construct_region(prevs, confidence_level, method=self.region)
+        prev_estim = conf.point_estimate()
+
+        return prev_estim, conf
+
+    def fit(self, X, y):
+        self.quantifier._check_init_parameters()
+        classif_predictions, labels = self.quantifier.classifier_fit_predict(X, y)
+        self.aggregation_fit(classif_predictions, labels)
+        return self
+
+    def quantify_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
+        predictions = self.quantifier.classify(instances)
+        return self.aggregate_conf(predictions, confidence_level=confidence_level)
+
+    @property
+    def classifier(self):
+        return self.quantifier.classifier
+
+    def _classifier_method(self):
+        return self.quantifier._classifier_method()
+
+
+class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
+    """
+    `Bayesian quantification <https://arxiv.org/abs/2302.09159>`_ method,
+    which is a variant of :class:`ACC` that calculates the posterior probability distribution
+    over the prevalence vectors, rather than providing a point estimate obtained
+    by matrix inversion.
+
+    Can be used to diagnose degeneracy in the predictions visible when the confusion
+    matrix has high condition number or to quantify uncertainty around the point estimate.
+
+    This method relies on extra dependencies, which have to be installed via:
+    `$ pip install quapy[bayes]`
+
+    :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 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. Set to
+        None when the method does not require any validation data, in order to avoid that some portion of
+        the training data be wasted.
+    :param num_warmup: number of warmup iterations for the MCMC sampler (default 500)
+    :param num_samples: number of samples to draw from the posterior (default 1000)
+    :param mcmc_seed: random seed for the MCMC sampler (default 0)
+    :param confidence_level: float in [0,1] to construct a confidence region around the point estimate (default 0.95)
+    :param region: string, set to `intervals` for constructing confidence intervals (default), or to
+        `ellipse` for constructing an ellipse in the probability simplex, or to `ellipse-clr` for
+        constructing an ellipse in the Centered-Log Ratio (CLR) unconstrained space.
+    """
+    def __init__(self,
+                 classifier: BaseEstimator=None,
+                 fit_classifier=True,
+                 val_split: int = 5,
+                 num_warmup: int = 500,
+                 num_samples: int = 1_000,
+                 mcmc_seed: int = 0,
+                 confidence_level: float = 0.95,
+                 region: str = 'intervals'):
+
+        if num_warmup <= 0:
+            raise ValueError(f'parameter {num_warmup=} must be a positive integer')
+        if num_samples <= 0:
+            raise ValueError(f'parameter {num_samples=} must be a positive integer')
+
+        if _bayesian.DEPENDENCIES_INSTALLED is False:
+            raise ImportError("Auxiliary dependencies are required. "
+                              "Run `$ pip install quapy[bayes]` to install them.")
+
+        super().__init__(classifier, fit_classifier, val_split)
+        self.num_warmup = num_warmup
+        self.num_samples = num_samples
+        self.mcmc_seed = mcmc_seed
+        self.confidence_level = confidence_level
+        self.region = region
+
+        # Array of shape (n_classes, n_predicted_classes,) where entry (y, c) is the number of instances
+        # labeled as class y and predicted as class c.
+        # By default, this array is set to None and later defined as part of the `aggregation_fit` phase
+        self._n_and_c_labeled = None
+
+        # Dictionary with posterior samples, set when `aggregate` is provided.
+        self._samples = None
+
+    def aggregation_fit(self, classif_predictions, labels):
+        """
+        Estimates the misclassification rates.
+
+        :param classif_predictions: array-like with the label predictions returned by the classifier
+        :param labels: array-like with the true labels associated to each classifier prediction
+        """
+        pred_labels = classif_predictions
+        true_labels = labels
+        self._n_and_c_labeled = confusion_matrix(
+            y_true=true_labels,
+            y_pred=pred_labels,
+            labels=self.classifier.classes_
+        ).astype(float)
+
+    def sample_from_posterior(self, classif_predictions):
+        if self._n_and_c_labeled is None:
+            raise ValueError("aggregation_fit must be called before sample_from_posterior")
+
+        n_c_unlabeled = F.counts_from_labels(classif_predictions, self.classifier.classes_).astype(float)
+
+        self._samples = _bayesian.sample_posterior(
+            n_c_unlabeled=n_c_unlabeled,
+            n_y_and_c_labeled=self._n_and_c_labeled,
+            num_warmup=self.num_warmup,
+            num_samples=self.num_samples,
+            seed=self.mcmc_seed,
+        )
+        return self._samples
+
+    def get_prevalence_samples(self):
+        if self._samples is None:
+            raise ValueError("sample_from_posterior must be called before get_prevalence_samples")
+        return self._samples[_bayesian.P_TEST_Y]
+
+    def get_conditional_probability_samples(self):
+        if self._samples is None:
+            raise ValueError("sample_from_posterior must be called before get_conditional_probability_samples")
+        return self._samples[_bayesian.P_C_COND_Y]
+
+    def aggregate(self, classif_predictions):
+        samples = self.sample_from_posterior(classif_predictions)[_bayesian.P_TEST_Y]
+        return np.asarray(samples.mean(axis=0), dtype=float)
+
+    def quantify_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
+        classif_predictions = self.classify(instances)
+        point_estimate = self.aggregate(classif_predictions)
+        samples = self.get_prevalence_samples()  # available after calling "aggregate" function
+        region = WithConfidenceABC.construct_region(samples, confidence_level=self.confidence_level, method=self.region)
+        return point_estimate, region

quapy/method/meta.py

@@ -1,6 +1,6 @@
 import itertools
 from copy import deepcopy
-from typing import Union
+from typing import Union, List
 import numpy as np
 from sklearn.linear_model import LogisticRegression
 from sklearn.metrics import f1_score, make_scorer, accuracy_score
@@ -12,7 +12,7 @@ from quapy import functional as F
 from quapy.data import LabelledCollection
 from quapy.model_selection import GridSearchQ
 from quapy.method.base import BaseQuantifier, BinaryQuantifier
-from quapy.method.aggregative import CC, ACC, PACC, HDy, EMQ, AggregativeQuantifier
+from quapy.method.aggregative import CC, ACC, PACC, HDy, EMQ, AggregativeQuantifier, AggregativeSoftQuantifier
 
 try:
     from . import _neural
@@ -52,19 +52,19 @@ class MedianEstimator2(BinaryQuantifier):
 
     def _delayed_fit(self, args):
         with qp.util.temp_seed(self.random_state):
-            params, training = args
+            params, X, y = args
             model = deepcopy(self.base_quantifier)
             model.set_params(**params)
-            model.fit(training)
+            model.fit(X, y)
             return model
 
-    def fit(self, training: LabelledCollection):
+    def fit(self, X, y):
-        self._check_binary(training, self.__class__.__name__)
+        self._check_binary(y, self.__class__.__name__)
 
         configs = qp.model_selection.expand_grid(self.param_grid)
         self.models = qp.util.parallel(
             self._delayed_fit,
-            ((params, training) for params in configs),
+            ((params, X, y) for params in configs),
             seed=qp.environ.get('_R_SEED', None),
             n_jobs=self.n_jobs
         )
@@ -72,12 +72,12 @@ class MedianEstimator2(BinaryQuantifier):
 
     def _delayed_predict(self, args):
         model, instances = args
-        return model.quantify(instances)
+        return model.predict(instances)
 
-    def quantify(self, instances):
+    def predict(self, X):
         prev_preds = qp.util.parallel(
             self._delayed_predict,
-            ((model, instances) for model in self.models),
+            ((model, X) for model in self.models),
             seed=qp.environ.get('_R_SEED', None),
             n_jobs=self.n_jobs
         )
@@ -95,7 +95,7 @@ class MedianEstimator(BinaryQuantifier):
     :param base_quantifier: the base, binary quantifier
     :param random_state: a seed to be set before fitting any base quantifier (default None)
     :param param_grid: the grid or parameters towards which the median will be computed
-    :param n_jobs: number of parllel workes
+    :param n_jobs: number of parallel workers
     """
     def __init__(self, base_quantifier: BinaryQuantifier, param_grid: dict, random_state=None, n_jobs=None):
         self.base_quantifier = base_quantifier
@@ -111,75 +111,33 @@ class MedianEstimator(BinaryQuantifier):
 
     def _delayed_fit(self, args):
         with qp.util.temp_seed(self.random_state):
-            params, training = args
+            params, X, y = args
             model = deepcopy(self.base_quantifier)
             model.set_params(**params)
-            model.fit(training)
+            model.fit(X, y)
             return model
 
-    def _delayed_fit_classifier(self, args):
+    def fit(self, X, y):
-        with qp.util.temp_seed(self.random_state):
+        self._check_binary(y, self.__class__.__name__)
-            cls_params, training = args
-            model = deepcopy(self.base_quantifier)
-            model.set_params(**cls_params)
-            predictions = model.classifier_fit_predict(training, predict_on=model.val_split)
-            return (model, predictions)
 
-    def _delayed_fit_aggregation(self, args):
+        configs = qp.model_selection.expand_grid(self.param_grid)
-        with qp.util.temp_seed(self.random_state):
+        self.models = qp.util.parallel(
-            ((model, predictions), q_params), training = args
+            self._delayed_fit,
-            model = deepcopy(model)
+            ((params, X, y) for params in configs),
-            model.set_params(**q_params)
+            seed=qp.environ.get('_R_SEED', None),
-            model.aggregation_fit(predictions, training)
+            n_jobs=self.n_jobs,
-            return model
+            asarray=False
-
+        )
-
-    def fit(self, training: LabelledCollection):
-        self._check_binary(training, self.__class__.__name__)
-
-        if isinstance(self.base_quantifier, AggregativeQuantifier):
-            cls_configs, q_configs = qp.model_selection.group_params(self.param_grid)
-
-            if len(cls_configs) > 1:
-                models_preds = qp.util.parallel(
-                    self._delayed_fit_classifier,
-                    ((params, training) for params in cls_configs),
-                    seed=qp.environ.get('_R_SEED', None),
-                    n_jobs=self.n_jobs,
-                    asarray=False
-                )
-            else:
-                model = self.base_quantifier
-                model.set_params(**cls_configs[0])
-                predictions = model.classifier_fit_predict(training, predict_on=model.val_split)
-                models_preds = [(model, predictions)]
-
-            self.models = qp.util.parallel(
-                self._delayed_fit_aggregation,
-                ((setup, training) for setup in itertools.product(models_preds, q_configs)),
-                seed=qp.environ.get('_R_SEED', None),
-                n_jobs=self.n_jobs,
-                asarray=False
-            )
-        else:
-            configs = qp.model_selection.expand_grid(self.param_grid)
-            self.models = qp.util.parallel(
-                self._delayed_fit,
-                ((params, training) for params in configs),
-                seed=qp.environ.get('_R_SEED', None),
-                n_jobs=self.n_jobs,
-                asarray=False
-            )
         return self
 
     def _delayed_predict(self, args):
         model, instances = args
-        return model.quantify(instances)
+        return model.predict(instances)
 
-    def quantify(self, instances):
+    def predict(self, X):
         prev_preds = qp.util.parallel(
             self._delayed_predict,
-            ((model, instances) for model in self.models),
+            ((model, X) for model in self.models),
             seed=qp.environ.get('_R_SEED', None),
             n_jobs=self.n_jobs,
             asarray=False
@@ -257,13 +215,14 @@ class Ensemble(BaseQuantifier):
         if self.verbose:
             print('[Ensemble]' + msg)
 
-    def fit(self, data: qp.data.LabelledCollection, val_split: Union[qp.data.LabelledCollection, float] = None):
+    def fit(self, X, y):
+
+        data = LabelledCollection(X, y)
 
         if self.policy == 'ds' and not data.binary:
             raise ValueError(f'ds policy is only defined for binary quantification, but this dataset is not binary')
 
-        if val_split is None:
+        val_split = self.val_split
-            val_split = self.val_split
 
         # randomly chooses the prevalences for each member of the ensemble (preventing classes with less than
         # min_pos positive examples)
@@ -294,15 +253,15 @@ class Ensemble(BaseQuantifier):
         self._sout('Fit [Done]')
         return self
 
-    def quantify(self, instances):
+    def predict(self, X):
         predictions = np.asarray(
-            qp.util.parallel(_delayed_quantify, ((Qi, instances) for Qi in self.ensemble), n_jobs=self.n_jobs)
+            qp.util.parallel(_delayed_quantify, ((Qi, X) for Qi in self.ensemble), n_jobs=self.n_jobs)
         )
 
         if self.policy == 'ptr':
             predictions = self._ptr_policy(predictions)
         elif self.policy == 'ds':
-            predictions = self._ds_policy(predictions, instances)
+            predictions = self._ds_policy(predictions, X)
 
         predictions = np.mean(predictions, axis=0)
         return F.normalize_prevalence(predictions)
@@ -455,22 +414,22 @@ def _delayed_new_instance(args):
     sample = data.sampling_from_index(sample_index)
 
     if val_split is not None:
-        model.fit(sample, val_split=val_split)
+        model.fit(*sample.Xy, val_split=val_split)
     else:
-        model.fit(sample)
+        model.fit(*sample.Xy)
 
     tr_prevalence = sample.prevalence()
     tr_distribution = get_probability_distribution(posteriors[sample_index]) if (posteriors is not None) else None
 
     if verbose:
-        print(f'\t\--fit-ended for prev {F.strprev(prev)}')
+        print(f'\t--fit-ended for prev {F.strprev(prev)}')
 
     return (model, tr_prevalence, tr_distribution, sample if keep_samples else None)
 
 
 def _delayed_quantify(args):
     quantifier, instances = args
-    return quantifier[0].quantify(instances)
+    return quantifier[0].predict(instances)
 
 
 def _draw_simplex(ndim, min_val, max_trials=100):
@@ -691,3 +650,107 @@ def EEMQ(classifier, param_grid=None, optim=None, param_mod_sel=None, **kwargs):
     """
 
     return ensembleFactory(classifier, EMQ, param_grid, optim, param_mod_sel, **kwargs)
+
+
+def merge(prev_predictions, merge_fun):
+    prev_predictions = np.asarray(prev_predictions)
+    if merge_fun == 'median':
+        prevalences = np.median(prev_predictions, axis=0)
+        prevalences = F.normalize_prevalence(prevalences, method='l1')
+    elif merge_fun == 'mean':
+        prevalences = np.mean(prev_predictions, axis=0)
+    else:
+        raise NotImplementedError(f'merge function {merge_fun} not implemented!')
+    return prevalences
+
+
+class SCMQ(AggregativeSoftQuantifier):
+
+    MERGE_FUNCTIONS = ['median', 'mean']
+
+    def __init__(self, classifier, quantifiers: List[AggregativeSoftQuantifier], merge_fun='median', val_split=5):
+        self.classifier = classifier
+        self.quantifiers = [deepcopy(q) for q in quantifiers]
+        assert merge_fun in self.MERGE_FUNCTIONS, f'unknown {merge_fun=}, valid ones are {self.MERGE_FUNCTIONS}'
+        self.merge_fun = merge_fun
+        self.val_split = val_split
+
+    def aggregation_fit(self, classif_predictions, labels):
+        for quantifier in self.quantifiers:
+            quantifier.classifier = self.classifier
+            quantifier.aggregation_fit(classif_predictions, labels)
+        return self
+
+    def aggregate(self, classif_predictions: np.ndarray):
+        prev_predictions = []
+        for quantifier_i in self.quantifiers:
+            prevalence_i = quantifier_i.aggregate(classif_predictions)
+            prev_predictions.append(prevalence_i)
+        return merge(prev_predictions, merge_fun=self.merge_fun)
+
+
+class MCSQ(BaseQuantifier):
+    def __init__(self, classifiers, quantifier: AggregativeSoftQuantifier, merge_fun='median', val_split=5):
+        self.merge_fun = merge_fun
+        self.val_split = val_split
+        self.mcsqs = []
+        for classifier in classifiers:
+            quantifier = deepcopy(quantifier)
+            quantifier.classifier = classifier
+            self.mcsqs.append(quantifier)
+
+    def fit(self, data: LabelledCollection):
+        for q in self.mcsqs:
+            q.fit(data, val_split=self.val_split)
+        return self
+
+    def quantify(self, instances):
+        prev_predictions = []
+        for q in self.mcsqs:
+            prevalence_i = q.quantify(instances)
+            prev_predictions.append(prevalence_i)
+        return merge(prev_predictions, merge_fun=self.merge_fun)
+
+
+class MCMQ(BaseQuantifier):
+    def __init__(self, classifiers, quantifiers: List[AggregativeSoftQuantifier], merge_fun='median', val_split=5):
+        self.merge_fun = merge_fun
+        self.scmqs = []
+        for classifier in classifiers:
+            self.scmqs.append(SCMQ(classifier, quantifiers, val_split=val_split))
+
+    def fit(self, data: LabelledCollection):
+        for q in self.scmqs:
+            q.fit(data)
+        return self
+
+    def quantify(self, instances):
+        prev_predictions = []
+        for q in self.scmqs:
+            prevalence_i = q.quantify(instances)
+            prev_predictions.append(prevalence_i)
+        return merge(prev_predictions, merge_fun=self.merge_fun)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

quapy/method/non_aggregative.py

@@ -20,21 +20,23 @@ class MaximumLikelihoodPrevalenceEstimation(BaseQuantifier):
     def __init__(self):
         self._classes_ = None
 
-    def fit(self, data: LabelledCollection):
+    def fit(self, X, y):
         """
         Computes the training prevalence and stores it.
 
-        :param data: the training sample
+        :param X: array-like of shape `(n_samples, n_features)`, the training instances
+        :param y: array-like of shape `(n_samples,)`, the labels
         :return: self
         """
-        self.estimated_prevalence = data.prevalence()
+        self._classes_ = F.classes_from_labels(labels=y)
+        self.estimated_prevalence = F.prevalence_from_labels(y, classes=self._classes_)
         return self
 
-    def quantify(self, instances):
+    def predict(self, X):
         """
         Ignores the input instances and returns, as the class prevalence estimantes, the training prevalence.
 
-        :param instances: array-like (ignored)
+        :param X: array-like (ignored)
         :return: the class prevalence seen during training
         """
         return self.estimated_prevalence
@@ -100,7 +102,7 @@ class DMx(BaseQuantifier):
 
         return distributions
 
-    def fit(self, data: LabelledCollection):
+    def fit(self, X, y):
         """
         Generates the validation distributions out of the training data (covariates).
         The validation distributions have shape `(n, nfeats, nbins)`, with `n` the number of classes, `nfeats`
@@ -109,33 +111,33 @@ class DMx(BaseQuantifier):
         training data labelled with class `i`; while `dij = di[j]` is the discrete distribution for feature j in
         training data labelled with class `i`, and `dij[k]` is the fraction of instances with a value in the `k`-th bin.
 
-        :param data: the training set
+        :param X: array-like of shape `(n_samples, n_features)`, the training instances
+        :param y: array-like of shape `(n_samples,)`, the labels
         """
-        X, y = data.Xy
-
         self.nfeats = X.shape[1]
         self.feat_ranges = _get_features_range(X)
+        n_classes = len(np.unique(y))
 
         self.validation_distribution = np.asarray(
-            [self.__get_distributions(X[y==cat]) for cat in range(data.n_classes)]
+            [self.__get_distributions(X[y==cat]) for cat in range(n_classes)]
         )
 
         return self
 
-    def quantify(self, instances):
+    def predict(self, X):
         """
         Searches for the mixture model parameter (the sought prevalence values) that yields a validation distribution
         (the mixture) that best matches the test distribution, in terms of the divergence measure of choice.
         The matching is computed as the average dissimilarity (in terms of the dissimilarity measure of choice)
         between all feature-specific discrete distributions.
 
-        :param instances: instances in the sample
+        :param X: instances in the sample
         :return: a vector of class prevalence estimates
         """
 
-        assert instances.shape[1] == self.nfeats, f'wrong shape; expected {self.nfeats}, found {instances.shape[1]}'
+        assert X.shape[1] == self.nfeats, f'wrong shape; expected {self.nfeats}, found {X.shape[1]}'
 
-        test_distribution = self.__get_distributions(instances)
+        test_distribution = self.__get_distributions(X)
         divergence = get_divergence(self.divergence)
         n_classes, n_feats, nbins = self.validation_distribution.shape
         def loss(prev):
@@ -147,53 +149,53 @@ class DMx(BaseQuantifier):
         return F.argmin_prevalence(loss, n_classes, method=self.search)
 
 
-class ReadMe(BaseQuantifier):
+# class ReadMe(BaseQuantifier):
-
+#
-    def __init__(self, bootstrap_trials=100, bootstrap_range=100, bagging_trials=100, bagging_range=25, **vectorizer_kwargs):
+#     def __init__(self, bootstrap_trials=100, bootstrap_range=100, bagging_trials=100, bagging_range=25, **vectorizer_kwargs):
-        raise NotImplementedError('under development ...')
+#         raise NotImplementedError('under development ...')
-        self.bootstrap_trials = bootstrap_trials
+#         self.bootstrap_trials = bootstrap_trials
-        self.bootstrap_range = bootstrap_range
+#         self.bootstrap_range = bootstrap_range
-        self.bagging_trials = bagging_trials
+#         self.bagging_trials = bagging_trials
-        self.bagging_range = bagging_range
+#         self.bagging_range = bagging_range
-        self.vectorizer_kwargs = vectorizer_kwargs
+#         self.vectorizer_kwargs = vectorizer_kwargs
-
+#
-    def fit(self, data: LabelledCollection):
+#     def fit(self, data: LabelledCollection):
-        X, y = data.Xy
+#         X, y = data.Xy
-        self.vectorizer = CountVectorizer(binary=True, **self.vectorizer_kwargs)
+#         self.vectorizer = CountVectorizer(binary=True, **self.vectorizer_kwargs)
-        X = self.vectorizer.fit_transform(X)
+#         X = self.vectorizer.fit_transform(X)
-        self.class_conditional_X = {i: X[y==i] for i in range(data.classes_)}
+#         self.class_conditional_X = {i: X[y==i] for i in range(data.classes_)}
-
+#
-    def quantify(self, instances):
+#     def predict(self, X):
-        X = self.vectorizer.transform(instances)
+#         X = self.vectorizer.transform(X)
-
+#
-        # number of features
+#         # number of features
-        num_docs, num_feats = X.shape
+#         num_docs, num_feats = X.shape
-
+#
-        # bootstrap
+#         # bootstrap
-        p_boots = []
+#         p_boots = []
-        for _ in range(self.bootstrap_trials):
+#         for _ in range(self.bootstrap_trials):
-            docs_idx = np.random.choice(num_docs, size=self.bootstra_range, replace=False)
+#             docs_idx = np.random.choice(num_docs, size=self.bootstra_range, replace=False)
-            class_conditional_X = {i: X[docs_idx] for i, X in self.class_conditional_X.items()}
+#             class_conditional_X = {i: X[docs_idx] for i, X in self.class_conditional_X.items()}
-            Xboot = X[docs_idx]
+#             Xboot = X[docs_idx]
-
+#
-            # bagging
+#             # bagging
-            p_bags = []
+#             p_bags = []
-            for _ in range(self.bagging_trials):
+#             for _ in range(self.bagging_trials):
-                feat_idx = np.random.choice(num_feats, size=self.bagging_range, replace=False)
+#                 feat_idx = np.random.choice(num_feats, size=self.bagging_range, replace=False)
-                class_conditional_Xbag = {i: X[:, feat_idx] for i, X in class_conditional_X.items()}
+#                 class_conditional_Xbag = {i: X[:, feat_idx] for i, X in class_conditional_X.items()}
-                Xbag = Xboot[:,feat_idx]
+#                 Xbag = Xboot[:,feat_idx]
-                p = self.std_constrained_linear_ls(Xbag, class_conditional_Xbag)
+#                 p = self.std_constrained_linear_ls(Xbag, class_conditional_Xbag)
-                p_bags.append(p)
+#                 p_bags.append(p)
-            p_boots.append(np.mean(p_bags, axis=0))
+#             p_boots.append(np.mean(p_bags, axis=0))
-
+#
-        p_mean = np.mean(p_boots, axis=0)
+#         p_mean = np.mean(p_boots, axis=0)
-        p_std  = np.std(p_bags, axis=0)
+#         p_std  = np.std(p_bags, axis=0)
-
+#
-        return p_mean
+#         return p_mean
-
+#
-
+#
-    def std_constrained_linear_ls(self, X, class_cond_X: dict):
+#     def std_constrained_linear_ls(self, X, class_cond_X: dict):
-        pass
+#         pass
 
 
 def _get_features_range(X):

quapy/model_selection.py

@@ -86,14 +86,14 @@ class GridSearchQ(BaseQuantifier):
         self.n_jobs = qp._get_njobs(n_jobs)
         self.raise_errors = raise_errors
         self.verbose = verbose
-        self.__check_error(error)
+        self.__check_error_measure(error)
         assert isinstance(protocol, AbstractProtocol), 'unknown protocol'
 
     def _sout(self, msg):
         if self.verbose:
             print(f'[{self.__class__.__name__}:{self.model.__class__.__name__}]: {msg}')
 
-    def __check_error(self, error):
+    def __check_error_measure(self, error):
         if error in qp.error.QUANTIFICATION_ERROR:
             self.error = error
         elif isinstance(error, str):
@@ -109,7 +109,7 @@ class GridSearchQ(BaseQuantifier):
 
         def job(cls_params):
             model.set_params(**cls_params)
-            predictions = model.classifier_fit_predict(self._training)
+            predictions = model.classifier_fit_predict(self._training_X, self._training_y)
             return predictions
 
         predictions, status, took = self._error_handler(job, cls_params)
@@ -123,7 +123,8 @@ class GridSearchQ(BaseQuantifier):
 
         def job(q_params):
             model.set_params(**q_params)
-            model.aggregation_fit(predictions, self._training)
+            P, y = predictions
+            model.aggregation_fit(P, y)
             score = evaluation.evaluate(model, protocol=self.protocol, error_metric=self.error)
             return score
 
@@ -136,7 +137,7 @@ class GridSearchQ(BaseQuantifier):
 
         def job(params):
             model.set_params(**params)
-            model.fit(self._training)
+            model.fit(self._training_X, self._training_y)
             score = evaluation.evaluate(model, protocol=self.protocol, error_metric=self.error)
             return score
 
@@ -159,17 +160,19 @@ class GridSearchQ(BaseQuantifier):
             return False
         return True
 
-    def _compute_scores_aggregative(self, training):
+    def _compute_scores_aggregative(self, X, y):
         # break down the set of hyperparameters into two: classifier-specific, quantifier-specific
         cls_configs, q_configs = group_params(self.param_grid)
 
         # train all classifiers and get the predictions
-        self._training = training
+        self._training_X = X
+        self._training_y = y
         cls_outs = qp.util.parallel(
             self._prepare_classifier,
             cls_configs,
             seed=qp.environ.get('_R_SEED', None),
-            n_jobs=self.n_jobs
+            n_jobs=self.n_jobs,
+            asarray=False
         )
 
         # filter out classifier configurations that yielded any error
@@ -194,9 +197,10 @@ class GridSearchQ(BaseQuantifier):
 
         return aggr_outs
 
-    def _compute_scores_nonaggregative(self, training):
+    def _compute_scores_nonaggregative(self, X, y):
         configs = expand_grid(self.param_grid)
-        self._training = training
+        self._training_X = X
+        self._training_y = y
         scores = qp.util.parallel(
             self._prepare_nonaggr_model,
             configs,
@@ -211,11 +215,12 @@ class GridSearchQ(BaseQuantifier):
         else:
             self._sout(f'error={status}')
 
-    def fit(self, training: LabelledCollection):
+    def fit(self, X, y):
         """ Learning routine. Fits methods with all combinations of hyperparameters and selects the one minimizing
             the error metric.
 
-        :param training: the training set on which to optimize the hyperparameters
+        :param X: array-like, training covariates
+        :param y: array-like, labels of training data
         :return: self
         """
 
@@ -231,9 +236,9 @@ class GridSearchQ(BaseQuantifier):
 
         self._sout(f'starting model selection with n_jobs={self.n_jobs}')
         if self._break_down_fit():
-            results = self._compute_scores_aggregative(training)
+            results = self._compute_scores_aggregative(X, y)
         else:
-            results = self._compute_scores_nonaggregative(training)
+            results = self._compute_scores_nonaggregative(X, y)
 
         self.param_scores_ = {}
         self.best_score_ = None
@@ -248,13 +253,13 @@ class GridSearchQ(BaseQuantifier):
                 self.param_scores_[str(params)] = status.status
                 self.error_collector.append(status)
 
-        tend = time()-tinit
+        self.fit_time_ = time()-tinit
 
         if self.best_score_ is None:
             raise ValueError('no combination of hyperparameters seemed to work')
 
         self._sout(f'optimization finished: best params {self.best_params_} (score={self.best_score_:.5f}) '
-                   f'[took {tend:.4f}s]')
+                   f'[took {self.fit_time_:.4f}s]')
 
         no_errors = len(self.error_collector)
         if no_errors>0:
@@ -266,7 +271,10 @@ class GridSearchQ(BaseQuantifier):
             if isinstance(self.protocol, OnLabelledCollectionProtocol):
                 tinit = time()
                 self._sout(f'refitting on the whole development set')
-                self.best_model_.fit(training + self.protocol.get_labelled_collection())
+                validation_collection = self.protocol.get_labelled_collection()
+                training_collection = LabelledCollection(X, y, classes=validation_collection.classes)
+                devel_collection = training_collection + validation_collection
+                self.best_model_.fit(*devel_collection.Xy)
                 tend = time() - tinit
                 self.refit_time_ = tend
             else:
@@ -275,15 +283,15 @@ class GridSearchQ(BaseQuantifier):
 
         return self
 
-    def quantify(self, instances):
+    def predict(self, X):
         """Estimate class prevalence values using the best model found after calling the :meth:`fit` method.
 
-        :param instances: sample contanining the instances
+        :param X: sample contanining the instances
         :return: a ndarray of shape `(n_classes)` with class prevalence estimates as according to the best model found
             by the model selection process.
         """
         assert hasattr(self, 'best_model_'), 'quantify called before fit'
-        return self.best_model().quantify(instances)
+        return self.best_model().predict(X)
 
     def set_params(self, **parameters):
         """Sets the hyper-parameters to explore.
@@ -364,8 +372,8 @@ def cross_val_predict(quantifier: BaseQuantifier, data: LabelledCollection, nfol
     total_prev = np.zeros(shape=data.n_classes)
 
     for train, test in data.kFCV(nfolds=nfolds, random_state=random_state):
-        quantifier.fit(train)
+        quantifier.fit(*train.Xy)
-        fold_prev = quantifier.quantify(test.X)
+        fold_prev = quantifier.predict(test.X)
         rel_size = 1. * len(test) / len(data)
         total_prev += fold_prev*rel_size
 

quapy/plot.py

@@ -23,21 +23,29 @@ def binary_diagonal(method_names, true_prevs, estim_prevs, pos_class=1, title=No
     indicating which class is to be taken as the positive class. (For multiclass quantification problems, other plots
     like the :meth:`error_by_drift` might be preferable though).
 
+    The format convention is as follows: `method_names`, `true_prevs`, and `estim_prevs` are array-like of the same
+    length, with the ith element describing the output of an independent experiment. The elements of `true_prevs`, and
+    `estim_prevs` are `ndarrays` with coherent shape for the same experiment. Experiments for the same method on
+    different datasets can be used, in which case the method name can appear more than once in `method_names`.
+
     :param method_names: array-like with the method names for each experiment
-    :param true_prevs: array-like with the true prevalence values (each being a ndarray with n_classes components) for
+    :param true_prevs: array-like with the true prevalence values for each experiment. Each entry is a ndarray of
-        each experiment
+        shape `(n_samples, n_classes)` components.
-    :param estim_prevs: array-like with the estimated prevalence values (each being a ndarray with n_classes components)
+    :param estim_prevs: array-like with the estimated prevalence values for each experiment. Each entry is a ndarray of
-        for each experiment
+        shape `(n_samples, n_classes)` components and `n_samples` must coincide with the corresponding entry in
-    :param pos_class: index of the positive class
+        `true_prevs`.
-    :param title: the title to be displayed in the plot
+    :param pos_class: index of the positive class (default 1)
-    :param show_std: whether or not to show standard deviations (represented by color bands). This might be inconvenient
+    :param title: the title to be displayed in the plot (default None)
+    :param show_std: whether to show standard deviations (represented by color bands). This might be inconvenient
         for cases in which many methods are compared, or when the standard deviations are high -- default True)
-    :param legend: whether or not to display the leyend (default True)
+    :param legend: whether to display the legend (default True)
-    :param train_prev: if indicated (default is None), the training prevalence (for the positive class) is hightlighted
+    :param train_prev: if indicated (default is None), the training prevalence (for the positive class) is highlighted
-        in the plot. This is convenient when all the experiments have been conducted in the same dataset.
+        in the plot. This is convenient when all the experiments have been conducted in the same dataset, or in
+        datasets with the same training prevalence.
     :param savepath: path where to save the plot. If not indicated (as default), the plot is shown.
     :param method_order: if indicated (default is None), imposes the order in which the methods are processed (i.e.,
         listed in the legend and associated with matplotlib colors).
+    :return: returns (fig, ax) matplotlib objects for eventual customisation
     """
     fig, ax = plt.subplots()
     ax.set_aspect('equal')
@@ -78,13 +86,9 @@ def binary_diagonal(method_names, true_prevs, estim_prevs, pos_class=1, title=No
 
     if legend:
         ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
-        # box = ax.get_position()
-        # ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
-        # ax.legend(loc='lower center',
-        #           bbox_to_anchor=(1, -0.5),
-        #           ncol=(len(method_names)+1)//2)
 
     _save_or_show(savepath)
+    return fig, ax
 
 
 def binary_bias_global(method_names, true_prevs, estim_prevs, pos_class=1, title=None, savepath=None):
@@ -92,14 +96,21 @@ def binary_bias_global(method_names, true_prevs, estim_prevs, pos_class=1, title
     Box-plots displaying the global bias (i.e., signed error computed as the estimated value minus the true value)
     for each quantification method with respect to a given positive class.
 
+    The format convention is as follows: `method_names`, `true_prevs`, and `estim_prevs` are array-like of the same
+    length, with the ith element describing the output of an independent experiment. The elements of `true_prevs`, and
+    `estim_prevs` are `ndarrays` with coherent shape for the same experiment. Experiments for the same method on
+    different datasets can be used, in which case the method name can appear more than once in `method_names`.
+
     :param method_names: array-like with the method names for each experiment
-    :param true_prevs: array-like with the true prevalence values (each being a ndarray with n_classes components) for
+    :param true_prevs: array-like with the true prevalence values for each experiment. Each entry is a ndarray of
-        each experiment
+        shape `(n_samples, n_classes)` components.
-    :param estim_prevs: array-like with the estimated prevalence values (each being a ndarray with n_classes components)
+    :param estim_prevs: array-like with the estimated prevalence values for each experiment. Each entry is a ndarray of
-        for each experiment
+        shape `(n_samples, n_classes)` components and `n_samples` must coincide with the corresponding entry in
+        `true_prevs`.
     :param pos_class: index of the positive class
-    :param title: the title to be displayed in the plot
+    :param title: the title to be displayed in the plot (default None)
     :param savepath: path where to save the plot. If not indicated (as default), the plot is shown.
+    :return: returns (fig, ax) matplotlib objects for eventual customisation
     """
 
     method_names, true_prevs, estim_prevs = _merge(method_names, true_prevs, estim_prevs)
@@ -120,25 +131,34 @@ def binary_bias_global(method_names, true_prevs, estim_prevs, pos_class=1, title
 
     _save_or_show(savepath)
 
+    return fig, ax
+
 
 def binary_bias_bins(method_names, true_prevs, estim_prevs, pos_class=1, title=None, nbins=5, colormap=cm.tab10,
                      vertical_xticks=False, legend=True, savepath=None):
     """
     Box-plots displaying the local bias (i.e., signed error computed as the estimated value minus the true value)
-    for different bins of (true) prevalence of the positive classs, for each quantification method.
+    for different bins of (true) prevalence of the positive class, for each quantification method.
+
+    The format convention is as follows: `method_names`, `true_prevs`, and `estim_prevs` are array-like of the same
+    length, with the ith element describing the output of an independent experiment. The elements of `true_prevs`, and
+    `estim_prevs` are `ndarrays` with coherent shape for the same experiment. Experiments for the same method on
+    different datasets can be used, in which case the method name can appear more than once in `method_names`.
 
     :param method_names: array-like with the method names for each experiment
-    :param true_prevs: array-like with the true prevalence values (each being a ndarray with n_classes components) for
+    :param true_prevs: array-like with the true prevalence values for each experiment. Each entry is a ndarray of
-        each experiment
+        shape `(n_samples, n_classes)` components.
-    :param estim_prevs: array-like with the estimated prevalence values (each being a ndarray with n_classes components)
+    :param estim_prevs: array-like with the estimated prevalence values for each experiment. Each entry is a ndarray of
-        for each experiment
+        shape `(n_samples, n_classes)` components and `n_samples` must coincide with the corresponding entry in
+        `true_prevs`.
     :param pos_class: index of the positive class
-    :param title: the title to be displayed in the plot
+    :param title: the title to be displayed in the plot (default None)
-    :param nbins: number of bins
+    :param nbins: number of bins (default 5)
     :param colormap: the matplotlib colormap to use (default cm.tab10)
     :param vertical_xticks: whether or not to add secondary grid (default is False)
     :param legend: whether or not to display the legend (default is True)
     :param savepath: path where to save the plot. If not indicated (as default), the plot is shown.
+    :return: returns (fig, ax) matplotlib objects for eventual customisation
     """
     from pylab import boxplot, plot, setp
 
@@ -210,13 +230,15 @@ def binary_bias_bins(method_names, true_prevs, estim_prevs, pos_class=1, title=N
 
     _save_or_show(savepath)
 
+    return fig, ax
+
 
 def error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
                    n_bins=20, error_name='ae', show_std=False,
                    show_density=True,
                    show_legend=True,
                    logscale=False,
-                   title=f'Quantification error as a function of distribution shift',
+                   title=None,
                    vlines=None,
                    method_order=None,
                    savepath=None):
@@ -227,11 +249,17 @@ def error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
     fare in different regions of the prior probability shift spectrum (e.g., in the low-shift regime vs. in the
     high-shift regime).
 
+    The format convention is as follows: `method_names`, `true_prevs`, and `estim_prevs` are array-like of the same
+    length, with the ith element describing the output of an independent experiment. The elements of `true_prevs`, and
+    `estim_prevs` are `ndarrays` with coherent shape for the same experiment. Experiments for the same method on
+    different datasets can be used, in which case the method name can appear more than once in `method_names`.
+
     :param method_names: array-like with the method names for each experiment
-    :param true_prevs: array-like with the true prevalence values (each being a ndarray with n_classes components) for
+    :param true_prevs: array-like with the true prevalence values for each experiment. Each entry is a ndarray of
-        each experiment
+        shape `(n_samples, n_classes)` components.
-    :param estim_prevs: array-like with the estimated prevalence values (each being a ndarray with n_classes components)
+    :param estim_prevs: array-like with the estimated prevalence values for each experiment. Each entry is a ndarray of
-        for each experiment
+        shape `(n_samples, n_classes)` components and `n_samples` must coincide with the corresponding entry in
+        `true_prevs`.
     :param tr_prevs: training prevalence of each experiment
     :param n_bins: number of bins in which the y-axis is to be divided (default is 20)
     :param error_name: a string representing the name of an error function (as defined in `quapy.error`, default is "ae")
@@ -239,12 +267,13 @@ def error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
     :param show_density: whether or not to display the distribution of experiments for each bin (default is True)
     :param show_density: whether or not to display the legend of the chart (default is True)
     :param logscale: whether or not to log-scale the y-error measure (default is False)
-    :param title: title of the plot (default is "Quantification error as a function of distribution shift")
+    :param title: title of the plot (default is None)
     :param vlines: array-like list of values (default is None). If indicated, highlights some regions of the space
         using vertical dotted lines.
     :param method_order: if indicated (default is None), imposes the order in which the methods are processed (i.e.,
         listed in the legend and associated with matplotlib colors).
     :param savepath: path where to save the plot. If not indicated (as default), the plot is shown.
+    :return: returns (fig, ax) matplotlib objects for eventual customisation
     """
 
     fig, ax = plt.subplots()
@@ -253,14 +282,14 @@ def error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
     x_error = qp.error.ae
     y_error = getattr(qp.error, error_name)
 
+    if method_order is None:
+        method_order = []
+
     # get all data as a dictionary {'m':{'x':ndarray, 'y':ndarray}} where 'm' is a method name (in the same
     # order as in method_order (if specified), and where 'x' are the train-test shifts (computed as according to
     # x_error function) and 'y' is the estim-test shift (computed as according to y_error)
     data = _join_data_by_drift(method_names, true_prevs, estim_prevs, tr_prevs, x_error, y_error, method_order)
 
-    if method_order is None:
-        method_order = method_names
-
     _set_colors(ax, n_methods=len(method_order))
 
     bins = np.linspace(0, 1, n_bins+1)
@@ -313,11 +342,11 @@ def error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
         ax2.spines['right'].set_color('g')
         ax2.tick_params(axis='y', colors='g')
     
-    ax.set(xlabel=f'Distribution shift between training set and test sample',
+    ax.set(xlabel=f'Prior shift between training set and test sample',
-           ylabel=f'{error_name.upper()} (true distribution, predicted distribution)',
+           ylabel=f'{error_name.upper()} (true prev, predicted prev)',
            title=title)
-    box = ax.get_position()
+    # box = ax.get_position()
-    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
+    # ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
     if vlines:
         for vline in vlines:
             ax.axvline(vline, 0, 1, linestyle='--', color='k')
@@ -327,14 +356,15 @@ def error_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
         #nice scale for the logaritmic axis
         ax.set_ylim(0,10 ** math.ceil(math.log10(max_y)))
     
-    
     if show_legend:
-        fig.legend(loc='lower center',
+        fig.legend(loc='center left',
                   bbox_to_anchor=(1, 0.5),
-                  ncol=(len(method_names)+1)//2)
+                  ncol=1)
 
     _save_or_show(savepath)
 
+    return fig, ax
+
 
 def brokenbar_supremacy_by_drift(method_names, true_prevs, estim_prevs, tr_prevs,
                                  n_bins=20, binning='isomerous',
@@ -350,11 +380,17 @@ def brokenbar_supremacy_by_drift(method_names, true_prevs, estim_prevs, tr_prevs
     plot is displayed on top, that displays the distribution of experiments for each bin (when binning="isometric") or
     the percentiles points of the distribution (when binning="isomerous").
 
+    The format convention is as follows: `method_names`, `true_prevs`, and `estim_prevs` are array-like of the same
+    length, with the ith element describing the output of an independent experiment. The elements of `true_prevs`, and
+    `estim_prevs` are `ndarrays` with coherent shape for the same experiment. Experiments for the same method on
+    different datasets can be used, in which case the method name can appear more than once in `method_names`.
+
     :param method_names: array-like with the method names for each experiment
-    :param true_prevs: array-like with the true prevalence values (each being a ndarray with n_classes components) for
+    :param true_prevs: array-like with the true prevalence values for each experiment. Each entry is a ndarray of
-        each experiment
+        shape `(n_samples, n_classes)` components.
-    :param estim_prevs: array-like with the estimated prevalence values (each being a ndarray with n_classes components)
+    :param estim_prevs: array-like with the estimated prevalence values for each experiment. Each entry is a ndarray of
-        for each experiment
+        shape `(n_samples, n_classes)` components and `n_samples` must coincide with the corresponding entry in
+        `true_prevs`.
     :param tr_prevs: training prevalence of each experiment
     :param n_bins: number of bins in which the y-axis is to be divided (default is 20)
     :param binning: type of binning, either "isomerous" (default) or "isometric"
@@ -371,13 +407,16 @@ def brokenbar_supremacy_by_drift(method_names, true_prevs, estim_prevs, tr_prevs
     :param method_order: if indicated (default is None), imposes the order in which the methods are processed (i.e.,
         listed in the legend and associated with matplotlib colors).
     :param savepath: path where to save the plot. If not indicated (as default), the plot is shown.
-    :return:
+    :return: returns (fig, ax) matplotlib objects for eventual customisation
     """
     assert binning in ['isomerous', 'isometric'], 'unknown binning type; valid types are "isomerous" and "isometric"'
 
     x_error = getattr(qp.error, x_error)
     y_error = getattr(qp.error, y_error)
 
+    if method_order is None:
+        method_order = []
+
     # get all data as a dictionary {'m':{'x':ndarray, 'y':ndarray}} where 'm' is a method name (in the same
     # order as in method_order (if specified), and where 'x' are the train-test shifts (computed as according to
     # x_error function) and 'y' is the estim-test shift (computed as according to y_error)
@@ -518,6 +557,8 @@ def brokenbar_supremacy_by_drift(method_names, true_prevs, estim_prevs, tr_prevs
 
     _save_or_show(savepath)
 
+    return fig, ax
+
 
 def _merge(method_names, true_prevs, estim_prevs):
     ndims = true_prevs[0].shape[1]
@@ -535,8 +576,9 @@ def _merge(method_names, true_prevs, estim_prevs):
 
 def _set_colors(ax, n_methods):
     NUM_COLORS = n_methods
-    cm = plt.get_cmap('tab20')
+    if NUM_COLORS>10:
-    ax.set_prop_cycle(color=[cm(1. * i / NUM_COLORS) for i in range(NUM_COLORS)])
+        cm = plt.get_cmap('tab20')
+        ax.set_prop_cycle(color=[cm(1. * i / NUM_COLORS) for i in range(NUM_COLORS)])
 
 
 def _save_or_show(savepath):
@@ -568,3 +610,39 @@ def _join_data_by_drift(method_names, true_prevs, estim_prevs, tr_prevs, x_error
             method_order.append(method)
 
     return data
+
+
+def calibration_plot(prob_classifier, X, y, nbins=10, savepath=None):
+    posteriors = prob_classifier.predict_proba(X)
+    assert posteriors.ndim==2, 'calibration plot only works for binary problems'
+    posteriors = posteriors[:,1]
+    pred_y = posteriors>=0.5
+    bins = np.linspace(0, 1, nbins + 1)
+    binned_values = np.digitize(posteriors, bins, right=False)
+    print(np.unique(binned_values))
+    correct = pred_y == y
+    bin_centers = (bins[:-1] + bins[1:]) / 2
+    bins_names = np.arange(nbins)
+    y_axis = [correct[binned_values==bin].mean() for bin in bins_names]
+    y_axis = [v if not np.isnan(v) else 0 for v in y_axis]
+    # Crear el gráfico de barras
+    plt.bar(bin_centers, y_axis, width=bins[1]-bins[0], edgecolor='black', alpha=0.7)
+
+    # Etiquetas y título
+    plt.xlabel("Bin")
+    plt.ylabel("Value")
+    plt.title("Bar plot of calculated values per bin")
+    plt.xticks(bin_centers, [f"{b:.2f}" for b in bin_centers], rotation=45)
+
+    # Mostrar el gráfico
+    plt.tight_layout()
+    plt.show()
+
+if __name__ == '__main__':
+    import quapy as qp
+    from sklearn.linear_model import LogisticRegression
+    data = qp.datasets.fetch_UCIBinaryDataset(qp.datasets.UCI_BINARY_DATASETS[6])
+    train, test = data.train_test
+    classifier = LogisticRegression()
+    classifier.fit(*train.Xy)
+    calibration_plot(classifier, *test.Xy)

quapy/protocol.py

@@ -1,4 +1,6 @@
 from copy import deepcopy
+from typing import Iterable
+
 import quapy as qp
 import numpy as np
 import itertools
@@ -62,6 +64,36 @@ class IterateProtocol(AbstractProtocol):
         return len(self.samples)
 
 
+class ProtocolFromIndex(AbstractProtocol):
+    """
+    A protocol from a list of indexes
+
+    :param data: a :class:`quapy.data.base.LabelledCollection`
+    :param indexes: a list of indexes
+    """
+    def __init__(self, data: LabelledCollection, indexes: Iterable):
+        self.data = data
+        self.indexes = indexes
+
+    def __call__(self):
+        """
+        Yields one sample at a time extracted using the indexes
+
+        :return: yields a tuple `(sample, prev) at a time, where `sample` is a set of instances
+            and in which `prev` is an `nd.array` with the class prevalence values
+        """
+        for index in self.indexes:
+            yield self.data.sampling_from_index(index).Xp
+
+    def total(self):
+        """
+        Returns the number of samples in this protocol
+
+        :return: int
+        """
+        return len(self.indexes)
+
+
 class AbstractStochasticSeededProtocol(AbstractProtocol):
     """
     An `AbstractStochasticSeededProtocol` is a protocol that generates, via any random procedure (e.g.,
@@ -124,9 +156,9 @@ class AbstractStochasticSeededProtocol(AbstractProtocol):
             if self.random_state is not None:
                 stack.enter_context(qp.util.temp_seed(self.random_state))
             for params in self.samples_parameters():
-                yield self.collator(self.sample(params))
+                yield self.collator(self.sample(params), params)
 
-    def collator(self, sample, *args):
+    def collator(self, sample, params):
         """
         The collator prepares the sample to accommodate the desired output format before returning the output.
         This collator simply returns the sample as it is. Classes inheriting from this abstract class can
@@ -191,9 +223,11 @@ class OnLabelledCollectionProtocol:
         assert return_type in cls.RETURN_TYPES, \
             f'unknown return type passed as argument; valid ones are {cls.RETURN_TYPES}'
         if return_type=='sample_prev':
-            return lambda lc:lc.Xp
+            return lambda lc,params:lc.Xp
         elif return_type=='labelled_collection':
-            return lambda lc:lc
+            return lambda lc,params:lc
+        elif return_type=='index':
+            return lambda lc,params:params
 
 
 class APP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):

quapy/tests/test_datasets.py

@@ -15,10 +15,13 @@ class TestDatasets(unittest.TestCase):
         return PCC(LogisticRegression(C=0.001, max_iter=100))
 
     def _check_dataset(self, dataset):
+        train, test = dataset.reduce().train_test
         q = self.new_quantifier()
         print(f'testing method {q} in {dataset.name}...', end='')
-        q.fit(dataset.training)
+        if len(train)>500:
-        estim_prevalences = q.quantify(dataset.test.instances)
+            train = train.sampling(500)
+        q.fit(*dataset.training.Xy)
+        estim_prevalences = q.predict(dataset.test.instances)
         self.assertTrue(F.check_prevalence_vector(estim_prevalences))
         print(f'[done]')
 
@@ -26,7 +29,7 @@ class TestDatasets(unittest.TestCase):
         for X, p in gen():
             if vectorizer is not None:
                 X = vectorizer.transform(X)
-            estim_prevalences = q.quantify(X)
+            estim_prevalences = q.predict(X)
             self.assertTrue(F.check_prevalence_vector(estim_prevalences))
             max_samples_test -= 1
             if max_samples_test == 0:
@@ -42,7 +45,9 @@ class TestDatasets(unittest.TestCase):
             self._check_dataset(dataset)
 
     def test_twitter(self):
-        for dataset_name in TWITTER_SENTIMENT_DATASETS_TEST:
+        # all the datasets are contained in the same resource; if the first one
+        # works, there is no need to test for the rest
+        for dataset_name in TWITTER_SENTIMENT_DATASETS_TEST[:1]:
             print(f'loading dataset {dataset_name}...', end='')
             dataset = fetch_twitter(dataset_name, min_df=10)
             dataset.stats()
@@ -52,18 +57,12 @@ class TestDatasets(unittest.TestCase):
 
     def test_UCIBinaryDataset(self):
         for dataset_name in UCI_BINARY_DATASETS:
-            try:
+            print(f'loading dataset {dataset_name}...', end='')
-                print(f'loading dataset {dataset_name}...', end='')
+            dataset = fetch_UCIBinaryDataset(dataset_name)
-                dataset = fetch_UCIBinaryDataset(dataset_name)
+            dataset.stats()
-                dataset.stats()
+            dataset.reduce()
-                dataset.reduce()
+            print(f'[done]')
-                print(f'[done]')
+            self._check_dataset(dataset)
-                self._check_dataset(dataset)
-            except FileNotFoundError as fnfe:
-                if dataset_name == 'pageblocks.5' and fnfe.args[0].find(
-                        'If this is the first time you attempt to load this dataset') > 0:
-                    print('The pageblocks.5 dataset requires some hand processing to be usable; skipping this test.')
-                    continue
 
     def test_UCIMultiDataset(self):
         for dataset_name in UCI_MULTICLASS_DATASETS:
@@ -83,18 +82,18 @@ class TestDatasets(unittest.TestCase):
             return
 
         for dataset_name in LEQUA2022_VECTOR_TASKS:
-            print(f'loading dataset {dataset_name}...', end='')
+            print(f'LeQu2022: loading dataset {dataset_name}...', end='')
             train, gen_val, gen_test = fetch_lequa2022(dataset_name)
             train.stats()
             n_classes = train.n_classes
             train = train.sampling(100, *F.uniform_prevalence(n_classes))
             q = self.new_quantifier()
-            q.fit(train)
+            q.fit(*train.Xy)
             self._check_samples(gen_val, q, max_samples_test=5)
             self._check_samples(gen_test, q, max_samples_test=5)
 
         for dataset_name in LEQUA2022_TEXT_TASKS:
-            print(f'loading dataset {dataset_name}...', end='')
+            print(f'LeQu2022: loading dataset {dataset_name}...', end='')
             train, gen_val, gen_test = fetch_lequa2022(dataset_name)
             train.stats()
             n_classes = train.n_classes
@@ -102,10 +101,26 @@ class TestDatasets(unittest.TestCase):
             tfidf = TfidfVectorizer()
             train.instances = tfidf.fit_transform(train.instances)
             q = self.new_quantifier()
-            q.fit(train)
+            q.fit(*train.Xy)
             self._check_samples(gen_val, q, max_samples_test=5, vectorizer=tfidf)
             self._check_samples(gen_test, q, max_samples_test=5, vectorizer=tfidf)
 
+    def test_lequa2024(self):
+        if os.environ.get('QUAPY_TESTS_OMIT_LARGE_DATASETS'):
+            print("omitting test_lequa2024 because QUAPY_TESTS_OMIT_LARGE_DATASETS is set")
+            return
+
+        for task in LEQUA2024_TASKS:
+            print(f'LeQu2024: loading task {task}...', end='')
+            train, gen_val, gen_test = fetch_lequa2024(task, merge_T3=True)
+            train.stats()
+            n_classes = train.n_classes
+            train = train.sampling(100, *F.uniform_prevalence(n_classes))
+            q = self.new_quantifier()
+            q.fit(*train.Xy)
+            self._check_samples(gen_val, q, max_samples_test=5)
+            self._check_samples(gen_test, q, max_samples_test=5)
+
 
     def test_IFCB(self):
         if os.environ.get('QUAPY_TESTS_OMIT_LARGE_DATASETS'):
@@ -119,7 +134,7 @@ class TestDatasets(unittest.TestCase):
             n_classes = train.n_classes
             train = train.sampling(100, *F.uniform_prevalence(n_classes))
             q = self.new_quantifier()
-            q.fit(train)
+            q.fit(*train.Xy)
             self._check_samples(gen, q, max_samples_test=5)
 
 

quapy/tests/test_evaluation.py

@@ -29,7 +29,7 @@ class EvalTestCase(unittest.TestCase):
                 time.sleep(1)
                 return super().predict_proba(X)
 
-        emq = EMQ(SlowLR()).fit(train)
+        emq = EMQ(SlowLR()).fit(*train.Xy)
 
         tinit = time()
         score = qp.evaluation.evaluate(emq, protocol, error_metric='mae', verbose=True, aggr_speedup='force')
@@ -41,14 +41,14 @@ class EvalTestCase(unittest.TestCase):
             def __init__(self, cls):
                 self.emq = EMQ(cls)
 
-            def quantify(self, instances):
+            def predict(self, X):
-                return self.emq.quantify(instances)
+                return self.emq.predict(X)
 
-            def fit(self, data):
+            def fit(self, X, y):
-                self.emq.fit(data)
+                self.emq.fit(X, y)
                 return self
 
-        emq = NonAggregativeEMQ(SlowLR()).fit(train)
+        emq = NonAggregativeEMQ(SlowLR()).fit(*train.Xy)
 
         tinit = time()
         score = qp.evaluation.evaluate(emq, protocol, error_metric='mae', verbose=True)
@@ -69,7 +69,7 @@ class EvalTestCase(unittest.TestCase):
 
         protocol = qp.protocol.APP(test, random_state=0)
 
-        q = PCC(LogisticRegression()).fit(train)
+        q = PCC(LogisticRegression()).fit(*train.Xy)
 
         single_errors = list(QUANTIFICATION_ERROR_SINGLE_NAMES)
         averaged_errors = ['m'+e for e in single_errors]

quapy/tests/test_hierarchy.py

@@ -9,9 +9,8 @@ import inspect
 class HierarchyTestCase(unittest.TestCase):
 
     def test_aggregative(self):
-        lr = LogisticRegression()
         for m in AGGREGATIVE_METHODS:
-            self.assertEqual(isinstance(m(lr), AggregativeQuantifier), True)
+            self.assertEqual(isinstance(m(), AggregativeQuantifier), True)
 
     def test_inspect_aggregative(self):
 
@@ -22,6 +21,7 @@ class HierarchyTestCase(unittest.TestCase):
         quantifiers = [cls for cls in classes if issubclass(cls, BaseQuantifier)]
         quantifiers = [cls for cls in quantifiers if issubclass(cls, AggregativeQuantifier)]
         quantifiers = [cls for cls in quantifiers if not inspect.isabstract(cls) ]
+        quantifiers = [cls for cls in quantifiers if cls is not OneVsAllAggregative]
 
         for cls in quantifiers:
             self.assertIn(cls, AGGREGATIVE_METHODS)

quapy/tests/test_methods.py

@@ -10,24 +10,27 @@ from quapy.method import AGGREGATIVE_METHODS, BINARY_METHODS, NON_AGGREGATIVE_ME
 from quapy.functional import check_prevalence_vector
 
 # a random selection of composed methods to test the qunfold integration
+from quapy.method.composable import check_compatible_qunfold_version
+
 from quapy.method.composable import (
     ComposableQuantifier,
     LeastSquaresLoss,
     HellingerSurrogateLoss,
-    ClassTransformer,
+    ClassRepresentation,
-    HistogramTransformer,
+    HistogramRepresentation,
-    CVClassifier,
+    CVClassifier
 )
+
 COMPOSABLE_METHODS = [
     ComposableQuantifier( # ACC
         LeastSquaresLoss(),
-        ClassTransformer(CVClassifier(LogisticRegression()))
+        ClassRepresentation(CVClassifier(LogisticRegression()))
     ),
     ComposableQuantifier( # HDy
         HellingerSurrogateLoss(),
-        HistogramTransformer(
+        HistogramRepresentation(
             3, # 3 bins per class
-            preprocessor = ClassTransformer(CVClassifier(LogisticRegression()))
+            preprocessor = ClassRepresentation(CVClassifier(LogisticRegression()))
         )
     ),
 ]
@@ -48,10 +51,10 @@ class TestMethods(unittest.TestCase):
                     print(f'skipping the test of binary model {model.__name__} on multiclass dataset {dataset.name}')
                     continue
 
-                q = model(learner)
+                q = model(learner, fit_classifier=False)
                 print('testing', q)
-                q.fit(dataset.training, fit_classifier=False)
+                q.fit(*dataset.training.Xy)
-                estim_prevalences = q.quantify(dataset.test.X)
+                estim_prevalences = q.predict(dataset.test.X)
                 self.assertTrue(check_prevalence_vector(estim_prevalences))
 
     def test_non_aggregative(self):
@@ -64,12 +67,11 @@ class TestMethods(unittest.TestCase):
 
                 q = model()
                 print(f'testing {q} on dataset {dataset.name}')
-                q.fit(dataset.training)
+                q.fit(*dataset.training.Xy)
-                estim_prevalences = q.quantify(dataset.test.X)
+                estim_prevalences = q.predict(dataset.test.X)
                 self.assertTrue(check_prevalence_vector(estim_prevalences))
 
     def test_ensembles(self):
-
         qp.environ['SAMPLE_SIZE'] = 10
 
         base_quantifier = ACC(LogisticRegression())
@@ -80,8 +82,8 @@ class TestMethods(unittest.TestCase):
 
             print(f'testing {base_quantifier} on dataset {dataset.name} with {policy=}')
             ensemble = Ensemble(quantifier=base_quantifier, size=3, policy=policy, n_jobs=-1)
-            ensemble.fit(dataset.training)
+            ensemble.fit(*dataset.training.Xy)
-            estim_prevalences = ensemble.quantify(dataset.test.instances)
+            estim_prevalences = ensemble.predict(dataset.test.instances)
             self.assertTrue(check_prevalence_vector(estim_prevalences))
 
     def test_quanet(self):
@@ -106,17 +108,22 @@ class TestMethods(unittest.TestCase):
         from quapy.method.meta import QuaNet
         model = QuaNet(learner, device='cpu', n_epochs=2, tr_iter_per_poch=10, va_iter_per_poch=10, patience=2)
 
-        model.fit(dataset.training)
+        model.fit(*dataset.training.Xy)
-        estim_prevalences = model.quantify(dataset.test.instances)
+        estim_prevalences = model.predict(dataset.test.instances)
         self.assertTrue(check_prevalence_vector(estim_prevalences))
 
     def test_composable(self):
-        for dataset in TestMethods.datasets:
+        if check_compatible_qunfold_version():
-            for q in COMPOSABLE_METHODS:
+            for dataset in TestMethods.datasets:
-                print('testing', q)
+                for q in COMPOSABLE_METHODS:
-                q.fit(dataset.training)
+                    print('testing', q)
-                estim_prevalences = q.quantify(dataset.test.X)
+                    q.fit(*dataset.training.Xy)
-                self.assertTrue(check_prevalence_vector(estim_prevalences))
+                    estim_prevalences = q.predict(dataset.test.X)
+                    print(estim_prevalences)
+                    self.assertTrue(check_prevalence_vector(estim_prevalences))
+        else:
+            from quapy.method.composable import __old_version_message
+            print(__old_version_message)
 
 
 if __name__ == '__main__':

quapy/tests/test_modsel.py

@@ -25,8 +25,8 @@ class ModselTestCase(unittest.TestCase):
         param_grid = {'classifier__C': [0.000001, 10.]}
         app = APP(validation, sample_size=100, random_state=1)
         q = GridSearchQ(
-            q, param_grid, protocol=app, error='mae', refit=True, timeout=-1, verbose=True
+            q, param_grid, protocol=app, error='mae', refit=False, timeout=-1, verbose=True, n_jobs=-1
-        ).fit(training)
+        ).fit(*training.Xy)
         print('best params', q.best_params_)
         print('best score', q.best_score_)
 
@@ -39,31 +39,30 @@ class ModselTestCase(unittest.TestCase):
         obtains the same optimal parameters
         """
 
-        q = PACC(LogisticRegression(random_state=1, max_iter=5000))
+        q = PACC(LogisticRegression(random_state=1, max_iter=3000))
 
-        data = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=10).reduce(n_train=500, random_state=1)
+        data = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=50)
         training, validation = data.training.split_stratified(0.7, random_state=1)
 
-        param_grid = {'classifier__C': np.logspace(-3,3,7)}
+        param_grid = {'classifier__C': np.logspace(-3,3,7), 'classifier__class_weight': ['balanced', None]}
         app = APP(validation, sample_size=100, random_state=1)
 
-        print('starting model selection in sequential exploration')
+        def do_gridsearch(n_jobs):
-        tinit = time.time()
+            print('starting model selection in sequential exploration')
-        modsel = GridSearchQ(
+            t_init = time.time()
-            q, param_grid, protocol=app, error='mae', refit=True, timeout=-1, n_jobs=1, verbose=True
+            modsel = GridSearchQ(
-        ).fit(training)
+                q, param_grid, protocol=app, error='mae', refit=False, timeout=-1, n_jobs=n_jobs, verbose=True
-        tend_seq = time.time()-tinit
+            ).fit(*training.Xy)
-        best_c_seq = modsel.best_params_['classifier__C']
+            t_end = time.time()-t_init
-        print(f'[done] took {tend_seq:.2f}s best C = {best_c_seq}')
+            best_c = modsel.best_params_['classifier__C']
+            print(f'[done] took {t_end:.2f}s best C = {best_c}')
+            return t_end, best_c
 
-        print('starting model selection in parallel exploration')
+        tend_seq, best_c_seq = do_gridsearch(n_jobs=1)
-        tinit = time.time()
+        tend_par, best_c_par = do_gridsearch(n_jobs=-1)
-        modsel = GridSearchQ(
+
-            q, param_grid, protocol=app, error='mae', refit=True, timeout=-1, n_jobs=-1, verbose=True
+        print(tend_seq, best_c_seq)
-        ).fit(training)
+        print(tend_par, best_c_par)
-        tend_par = time.time() - tinit
-        best_c_par = modsel.best_params_['classifier__C']
-        print(f'[done] took {tend_par:.2f}s best C = {best_c_par}')
 
         self.assertEqual(best_c_seq, best_c_par)
         self.assertLess(tend_par, tend_seq)
@@ -90,7 +89,7 @@ class ModselTestCase(unittest.TestCase):
             q, param_grid, protocol=app, timeout=3, n_jobs=-1, verbose=True, raise_errors=True
         )
         with self.assertRaises(TimeoutError):
-            modsel.fit(training)
+            modsel.fit(*training.Xy)
 
         print('Expecting ValueError to be raised')
         modsel = GridSearchQ(
@@ -99,7 +98,7 @@ class ModselTestCase(unittest.TestCase):
         with self.assertRaises(ValueError):
             # this exception is not raised because of the timeout, but because no combination of hyperparams
             # succedded (in this case, a ValueError is raised, regardless of "raise_errors"
-            modsel.fit(training)
+            modsel.fit(*training.Xy)
 
 
 if __name__ == '__main__':

quapy/tests/test_protocols.py

@@ -71,7 +71,7 @@ class TestProtocols(unittest.TestCase):
         # surprisingly enough, for some n_prevalences the test fails, notwithstanding
         # everything is correct. The problem is that in function APP.prevalence_grid()
         # there is sometimes one rounding error that gets cumulated and
-        # surpasses 1.0 (by a very small float value, 0.0000000000002 or sthe like)
+        # surpasses 1.0 (by a very small float value, 0.0000000000002 or the like)
         # so these tuples are mistakenly removed... I have tried with np.close, and
         # other workarounds, but eventually happens that there is some negative probability
         # in the sampling function...

quapy/tests/test_replicability.py

@@ -13,17 +13,18 @@ class TestReplicability(unittest.TestCase):
     def test_prediction_replicability(self):
 
         dataset = qp.datasets.fetch_UCIBinaryDataset('yeast')
+        train, test = dataset.train_test
 
         with qp.util.temp_seed(0):
             lr = LogisticRegression(random_state=0, max_iter=10000)
             pacc = PACC(lr)
-            prev = pacc.fit(dataset.training).quantify(dataset.test.X)
+            prev = pacc.fit(*train.Xy).predict(test.X)
             str_prev1 = strprev(prev, prec=5)
 
         with qp.util.temp_seed(0):
             lr = LogisticRegression(random_state=0, max_iter=10000)
             pacc = PACC(lr)
-            prev2 = pacc.fit(dataset.training).quantify(dataset.test.X)
+            prev2 = pacc.fit(*train.Xy).predict(test.X)
             str_prev2 = strprev(prev2, prec=5)
 
         self.assertEqual(str_prev1, str_prev2)
@@ -83,19 +84,19 @@ class TestReplicability(unittest.TestCase):
         test = test.sampling(500, *[0.1, 0.0, 0.1, 0.1, 0.2, 0.5, 0.0])
 
         with qp.util.temp_seed(10):
-            pacc = PACC(LogisticRegression(), val_split=2, n_jobs=2)
+            pacc = PACC(LogisticRegression(), val_split=.5, n_jobs=2)
-            pacc.fit(train, val_split=0.5)
+            pacc.fit(*train.Xy)
-            prev1 = F.strprev(pacc.quantify(test.instances))
+            prev1 = F.strprev(pacc.predict(test.instances))
 
         with qp.util.temp_seed(0):
-            pacc = PACC(LogisticRegression(), val_split=2, n_jobs=2)
+            pacc = PACC(LogisticRegression(), val_split=.5, n_jobs=2)
-            pacc.fit(train, val_split=0.5)
+            pacc.fit(*train.Xy)
-            prev2 = F.strprev(pacc.quantify(test.instances))
+            prev2 = F.strprev(pacc.predict(test.instances))
 
         with qp.util.temp_seed(0):
-            pacc = PACC(LogisticRegression(), val_split=2, n_jobs=2)
+            pacc = PACC(LogisticRegression(), val_split=.5, n_jobs=2)
-            pacc.fit(train, val_split=0.5)
+            pacc.fit(*train.Xy)
-            prev3 = F.strprev(pacc.quantify(test.instances))
+            prev3 = F.strprev(pacc.predict(test.instances))
 
         print(prev1)
         print(prev2)

quapy/util.py

@@ -208,7 +208,7 @@ def save_text_file(path, text):
     :param text: text to save.
     """
     create_parent_dir(path)
-    with open(text, 'wt') as fout:
+    with open(path, 'wt') as fout:
         fout.write(text)