diff --git a/README.md b/README.md
index 730a433..fcb820d 100644
--- a/README.md
+++ b/README.md
@@ -13,7 +13,7 @@ for facilitating the analysis and interpretation of the experimental results.
 
 ### Last updates:
 
-* Version 0.2.0 is released! major changes can be consulted [here](CHANGE_LOG.txt).
+* Version 0.2.1 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)
 * Manuals are available [here](https://hlt-isti.github.io/QuaPy/manuals.html)
 
diff --git a/docs/source/_static/custom.css b/docs/source/_static/custom.css
new file mode 100644
index 0000000..600ce46
--- /dev/null
+++ b/docs/source/_static/custom.css
@@ -0,0 +1,53 @@
+.navbar-brand.logo .title.logo__title {
+    display: none;
+}
+
+.navbar-brand img {
+    max-height: 2.2rem;
+    width: auto;
+}
+
+.navbar-brand.logo .title.logo__title {
+    font-size: 1.05rem;
+    line-height: 1.15;
+    white-space: nowrap;
+}
+
+@media (max-width: 1200px) {
+    .bd-header .navbar-header-items {
+        min-width: 0;
+    }
+
+    .bd-header .navbar-header-items__center {
+        overflow-x: auto;
+    }
+
+    .bd-header .bd-navbar-elements {
+        flex-wrap: nowrap;
+    }
+}
+
+.hero-copy {
+    font-size: 1.15rem;
+    line-height: 1.7;
+    max-width: 56rem;
+    margin: 0 0 1.5rem 0;
+}
+
+.landing-grid {
+    margin: 1.2rem 0 2rem 0;
+}
+
+.landing-card {
+    border-radius: 1rem;
+    border: 1px solid var(--pst-color-border, #d0d7de);
+    box-shadow: 0 10px 24px rgba(15, 23, 42, 0.08);
+}
+
+.landing-card .sd-card-title {
+    font-size: 1.1rem;
+}
+
+[data-theme="dark"] .landing-card {
+    box-shadow: 0 10px 24px rgba(0, 0, 0, 0.22);
+}
diff --git a/docs/source/_static/quapy_logo.png b/docs/source/_static/quapy_logo.png
new file mode 100644
index 0000000..3798993
Binary files /dev/null and b/docs/source/_static/quapy_logo.png differ
diff --git a/docs/source/_static/quapy_logo_dark.png b/docs/source/_static/quapy_logo_dark.png
new file mode 100644
index 0000000..18ffbad
Binary files /dev/null and b/docs/source/_static/quapy_logo_dark.png differ
diff --git a/examples/19.visualizing_simplex.py b/examples/19.visualizing_simplex.py
new file mode 100644
index 0000000..1ec5f9c
--- /dev/null
+++ b/examples/19.visualizing_simplex.py
@@ -0,0 +1,68 @@
+import numpy as np
+import quapy as qp
+from quapy.method.confidence import ConfidenceIntervals
+
+
+"""
+A minimal example showing how to visualise ternary prevalences on the simplex.
+The plot combines a cloud of posterior triplets, a few reference prevalences, a
+confidence ellipse induced by the cloud, and a smooth density centred around the
+true prevalence.
+"""
+
+rng = np.random.default_rng(0)
+true_prev = np.array([0.20, 0.35, 0.45])
+train_prev = np.array([0.50, 0.30, 0.20])
+pred_prev = np.array([0.18, 0.39, 0.43])
+posterior_cloud = rng.dirichlet(alpha=45 * true_prev, size=250)
+
+point_layers = [
+    {
+        'points': posterior_cloud,
+        'label': 'posterior cloud',
+        'style': {'s': 12, 'alpha': 0.25, 'color': 'steelblue', 'edgecolors': 'none'},
+    },
+    {
+        'points': true_prev,
+        'label': 'true prevalence',
+        'style': {'s': 70, 'color': 'black'},
+    },
+    {
+        'points': pred_prev,
+        'label': 'predicted prevalence',
+        'style': {'s': 70, 'color': 'crimson'},
+    },
+    {
+        'points': train_prev,
+        'label': 'training prevalence',
+        'style': {'s': 70, 'color': 'darkorange'},
+    },
+]
+
+confidence_region = ConfidenceIntervals(posterior_cloud, confidence_level=0.95, bonferroni_correction=True)
+
+region_layers = [
+    {
+        'fn': lambda p: float(p in confidence_region),
+        'label': '95% confidence intervals',
+        'color': 'seagreen',
+        'alpha': 0.15,
+    }
+]
+
+density = lambda p: np.exp(-45 * np.sum((p - true_prev) ** 2, axis=1))
+
+qp.plot.plot_simplex(
+    point_layers=point_layers,
+    region_layers=region_layers,
+    density_function=density,
+    density_color='royalblue',
+    class_names=['class A', 'class B', 'class C'],
+    title='Ternary prevalence visualisation',
+    legend_ncol=3,
+    figsize=(7.2, 5.8),
+    class_name_fontsize=9,
+    title_fontsize=10,
+    legend_fontsize=8,
+    savepath='./plots/simplex_visualization.png',
+)
diff --git a/examples/20.cifar10_quantification.py b/examples/20.cifar10_quantification.py
new file mode 100644
index 0000000..55e8f3f
--- /dev/null
+++ b/examples/20.cifar10_quantification.py
@@ -0,0 +1,64 @@
+import quapy as qp
+from quapy.data.datasets import fetch_image_embeddings
+from quapy.method.aggregative import EMQ, RLLS
+from quapy.classification.calibration import TemperatureScalingFromLogits
+from quapy.protocol import UPP
+from sklearn.linear_model import LogisticRegression
+
+
+# This example illustrates how to run experiments with image datasets, in this case with CIFAR10
+# The datasets available in quapy do not consist of raw image files, but are instead pre-generated 
+# embeddings (see the manuals for further information).
+
+if __name__ == '__main__':
+
+    # Let us begin with a typical case in which the embeddings come from the penultimate layer of a neural
+    # model (in this case, a resnet18). We get these representations by specifying embedding='features'
+
+    print('fetching cifar10 embeddings')
+    train, test = fetch_image_embeddings(dataset_name='cifar10', embedding='features').train_test
+
+    print('training:', train)
+    print('test:', test)
+
+    Xtr, ytr = train.Xy
+
+    # let us train an Expectation Maximazion Quantifier (EMQ), aka Maximum Likelihood for Label Shift (MLLS)
+    # using a logistic regressor as the underlying classifier, with Bias Corrected Temperature Scaling (BCTS)
+
+    bcts_emq = EMQ(classifier=LogisticRegression(), calib='bcts', val_split=5)
+
+    print(f'fitting quantifier {bcts_emq}')
+    bcts_emq.fit(Xtr, ytr)
+
+    # we generate many samples exhibiting prior probability shift with the artificial prevalence protocol
+    # (we use the multiclass variant UPP instead of the grid-based APP)
+
+    qp.environ["SAMPLE_SIZE"] = 500  # when the sample size is common to all experiments, it is conveniet to set it once and for all
+    artificial_prev_prot = UPP(test, repeats=200)
+    print('generating 200 test bags of 500 instances each')
+
+    bctsemq_report = qp.evaluation.evaluation_report(bcts_emq, protocol=artificial_prev_prot, error_metrics=['mae', 'mrae'])
+    print(bctsemq_report.mean(numeric_only=True))
+
+    # we could instead use the pre-generated logits of the resnet18
+    print('fetching cifar10 logits')
+    train, test = fetch_image_embeddings(dataset_name='cifar10', embedding='logits').train_test
+    Xtr, ytr = train.Xy
+
+    # in this case, the representations are already classification-related outputs;
+    # we can convert them into (hopefully well-) calibrated outputs via TemperatureScaling
+    
+    print('generating posterior probabilities out of logits via temperature scaling')    
+    emq = EMQ(classifier=TemperatureScalingFromLogits(bias_corrected=True))
+    emq.fit(Xtr, ytr)
+
+    print('generating 200 test bags of 500 instances each')
+    artificial_prev_prot = UPP(test, repeats=200)
+    emq_report = qp.evaluation.evaluation_report(emq, protocol=artificial_prev_prot, error_metrics=['mae', 'mrae'])
+    print(emq_report.mean(numeric_only=True))
+
+
+
+
+
diff --git a/quapy/data/datasets.py b/quapy/data/datasets.py
index 63e9314..fd2383a 100644
--- a/quapy/data/datasets.py
+++ b/quapy/data/datasets.py
@@ -1113,8 +1113,6 @@ def _fetch_image_embedding_splits(dataset_name, embedding, data_home=None) -> tu
     val = LabelledCollection(embedding_dict['val'], labels_dict['val'], classes=train.classes)
     test = LabelledCollection(embedding_dict['test'], labels_dict['test'], classes=train.classes)
 
-    print(f'{len(train)} | {len(val)} | {len(test)} | {train.X.shape[1]} | {train.n_classes} | {train.n_classes}')
-
     return train, val, test
 
 
diff --git a/quapy/method/_helper.py b/quapy/method/_helper.py
new file mode 100644
index 0000000..0e5792a
--- /dev/null
+++ b/quapy/method/_helper.py
@@ -0,0 +1,107 @@
+"""
+Internal helper utilities shared by quantification methods.
+"""
+
+import numpy as np
+from sklearn.metrics import confusion_matrix
+from sklearn.preprocessing import LabelEncoder
+
+
+def _get_abstention_calibrators():
+    try:
+        from abstention.calibration import NoBiasVectorScaling, TempScaling, VectorScaling
+    except ImportError as exc:
+        raise ImportError(
+            "Posterior calibration for EMQ requires the optional 'abstention' package."
+        ) from exc
+    return {
+        'nbvs': NoBiasVectorScaling(),
+        'bcts': TempScaling(bias_positions='all'),
+        'ts': TempScaling(),
+        'vs': VectorScaling(),
+    }
+
+
+def _get_cvxpy():
+    try:
+        import cvxpy as cp
+    except ImportError as exc:
+        raise ImportError(
+            "RLLS requires the optional 'cvxpy' package."
+        ) from exc
+    return cp
+
+
+def _labels_to_indices(labels, classes):
+    encoder = LabelEncoder().fit(classes)
+    return encoder.transform(labels)
+
+
+def _rlls_check_mode(mode):
+    valid = {'soft', 'hard'}
+    if mode not in valid:
+        raise ValueError(f'unknown mode {mode!r}; valid ones are {valid}')
+    return mode
+
+
+def _rlls_joint_distribution(posteriors, labels, classes, mode='soft'):
+    mode = _rlls_check_mode(mode)
+    posteriors = np.asarray(posteriors, dtype=float)
+    labels = np.asarray(labels)
+    n_samples, n_classes = posteriors.shape
+    assert n_classes == len(classes), 'wrong number of posterior columns'
+
+    if mode == 'hard':
+        pred = np.argmax(posteriors, axis=1)
+        encoded_labels = _labels_to_indices(labels, classes)
+        joint = confusion_matrix(encoded_labels, pred, labels=np.arange(n_classes)).T.astype(float)
+        return joint / n_samples
+
+    joint = np.zeros((n_classes, n_classes), dtype=float)
+    for class_index, class_ in enumerate(classes):
+        idx = labels == class_
+        if idx.any():
+            joint[:, class_index] = posteriors[idx].sum(axis=0)
+    return joint / n_samples
+
+
+def _rlls_predicted_marginal(posteriors, mode='soft'):
+    mode = _rlls_check_mode(mode)
+    posteriors = np.asarray(posteriors, dtype=float)
+    if mode == 'soft':
+        return posteriors.mean(axis=0)
+
+    pred = np.argmax(posteriors, axis=1)
+    counts = np.bincount(pred, minlength=posteriors.shape[1]).astype(float)
+    return counts / counts.sum()
+
+
+def _rlls_compute_3deltaC(n_classes, n_train, delta):
+    return 3 * (
+        2 * np.log(2 * n_classes / delta) / (3 * n_train)
+        + np.sqrt(2 * np.log(2 * n_classes / delta) / n_train)
+    )
+
+
+def _rlls_compute_weights(C_zy, qz, pz, rho, clip=False):
+    cp = _get_cvxpy()
+
+    n_classes = C_zy.shape[0]
+    theta = cp.Variable(n_classes)
+    b = qz - pz
+    objective = cp.Minimize(cp.pnorm(C_zy @ theta - b) + rho * cp.pnorm(theta))
+    constraints = [-1 <= theta]
+    problem = cp.Problem(objective, constraints)
+
+    try:
+        problem.solve(verbose=False, solver=cp.SCS)
+    except cp.error.SolverError:
+        problem.solve(verbose=False, solver=cp.SCS, use_indirect=True)
+
+    if theta.value is None:
+        raise RuntimeError('RLLS optimization failed to produce a solution')
+
+    w = 1 + np.asarray(theta.value, dtype=float)
+    if clip and np.any(w < 0):
+        w = np.clip(w, 0, None)
+    return w
diff --git a/quapy/tests/test_functional.py b/quapy/tests/test_functional.py
new file mode 100644
index 0000000..97e2307
--- /dev/null
+++ b/quapy/tests/test_functional.py
@@ -0,0 +1,27 @@
+import unittest
+
+import numpy as np
+
+import quapy.functional as F
+
+
+class TestFunctional(unittest.TestCase):
+
+    def test_ternary_search_binary(self):
+        def loss(prev):
+            return (prev[1] - 0.37) ** 2
+
+        result = F.argmin_prevalence(loss, n_classes=2, method='ternary_search')
+        self.assertTrue(np.allclose(result.sum(), 1.0))
+        self.assertAlmostEqual(result[1], 0.37, places=3)
+
+    def test_ternary_search_multiclass_not_supported(self):
+        def loss(prev):
+            return np.sum((prev - np.array([0.2, 0.3, 0.5])) ** 2)
+
+        with self.assertRaises(AssertionError):
+            F.argmin_prevalence(loss, n_classes=3, method='ternary_search')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/quapy/tests/test_plot.py b/quapy/tests/test_plot.py
new file mode 100644
index 0000000..3f5a5a5
--- /dev/null
+++ b/quapy/tests/test_plot.py
@@ -0,0 +1,41 @@
+import unittest
+
+try:
+    import matplotlib
+    matplotlib.use('Agg')
+    import matplotlib.pyplot as plt
+    HAS_MATPLOTLIB = True
+except ImportError:
+    plt = None
+    HAS_MATPLOTLIB = False
+
+import numpy as np
+
+import quapy as qp
+
+
+@unittest.skipUnless(HAS_MATPLOTLIB and qp.plot is not None, 'matplotlib is not available')
+class TestPlot(unittest.TestCase):
+
+    def test_plot_simplex_smoke(self):
+        rng = np.random.default_rng(0)
+        true_prev = np.array([0.2, 0.3, 0.5])
+        cloud = rng.dirichlet(alpha=30 * true_prev, size=50)
+
+        fig, ax = plt.subplots(figsize=(5, 5))
+        fig, ax = qp.plot.plot_simplex(
+            point_layers=[
+                {'points': cloud, 'label': 'cloud', 'style': {'s': 8, 'alpha': 0.2}},
+                {'points': true_prev, 'label': 'target', 'style': {'s': 50, 'color': 'black'}},
+            ],
+            region_layers=[
+                {'fn': lambda p: p[:, 2] >= 0.4, 'label': 'high class-3', 'color': 'green', 'alpha': 0.2},
+            ],
+            density_function=lambda p: np.exp(-25 * np.sum((p - true_prev) ** 2, axis=1)),
+            class_names=['A', 'B', 'C'],
+            ax=ax,
+        )
+
+        self.assertIs(fig, ax.figure)
+        self.assertGreaterEqual(len(ax.collections), 2)
+        plt.close(fig)