Add empirical attainment surface

package/visualization/plot_empirical_attainment_surface/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2025 Shuhei Watanabe
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

package/visualization/plot_empirical_attainment_surface/README.md

@@ -0,0 +1,100 @@
+---
+author: Shuhei Watanabe
+title: Visualizing Variability of Pareto Fronts over Multiple Runs (Empirical Attainment Surface)
+description: This module enables visualizing the uncertainty bands for bi-objective problems.
+tags: [visualization, multi-objective optimization, empirical attainment surface]
+optuna_versions: [4.5.0]
+license: MIT License
+---
+
+## Abstract
+
+Hyperparameter optimization is crucial to achieving high performance in deep learning.
+On top of the performance, other criteria such as inference time or memory requirement often need to be optimized due to some practical reasons.
+This motivates research on multi-objective optimization (MOO).
+However, Pareto fronts of MOO methods are often shown without considering the variability caused by random seeds, making the performance stability evaluation difficult.
+This package provides empirical attainment surface implementation based on [the original implementation](https://github.com/nabenabe0928/empirical-attainment-func).
+
+The details of empirical attainment surface are available in [`Python Tool for Visualizing Variability of Pareto Fronts over Multiple Runs`](https://arxiv.org/abs/2305.08852).
+
+## APIs
+
+- `plot_empirical_attainment_surface(study_list: list[optuna.Study], attainment_ratios: list[int], ax: plt.Axes | None = None, color: str | None = None, label: str | None = None, linestyle: str | None = None, marker: str | None = None, log_scale_inds: list[int] | None = None, **ax_plot_kwargs)`
+  - `study_list`: A list of studies used for the uncertainty visualization of the attainment surface.
+  - `attainment_ratios`: The ratios of studies that need to dominate the surface. The length must be either 1 or 3. If the length is 1, the corresponding attainment surface will be visualized. If the length is 3, the area between the attainment surfaces with the lowest and the highest ratios will be visualized as the uncertainty.
+  - `ax`: The axes of Matplotlib.
+  - `color`: The color of the empirical attainment surface.
+  - `label`: The label of the empirical attainment surface.
+  - `linestyle`: The line style of the empirical attainment surface.
+  - `marker`: The marker of the empirical attainment surface.
+  - `log_scale_inds`: The indices of values that should be in the log scale. Either 0 or 1 can be specified.
+  - `**ax_plot_kwargs`: The `kwargs` input to `ax.plot`.
+- `plot_multiple_empirical_attainment_surfaces(multiple_study_list: list[list[optuna.Study]], attainment_ratios: list[int], ax: plt.Axes | None = None, colors: list[str] | None = None, labels: list[str] | None = None, linestyles: list[str] | None = None, markers: list[str] | None = None, log_scale_inds: list[int] | None = None, **ax_plot_kwargs)`
+  - `multiple_study_list`: A list of study list used for the uncertainty visualization of each attainment surface.
+  - `attainment_ratios`: The ratios of studies that need to dominate the surface. The length must be either 1 or 3. If the length is 1, the corresponding attainment surface will be visualized. If the length is 3, the area between the attainment surfaces with the lowest and the highest ratios will be visualized as the uncertainty.
+  - `ax`: The axes of Matplotlib.
+  - `colors`: A list of the color of each empirical attainment surface.
+  - `labels`: A list of the label of each empirical attainment surface.
+  - `linestyles`: A list of the line style of each empirical attainment surface.
+  - `markers`: A list of the marker of each empirical attainment surface.
+  - `log_scale_inds`: The indices of values that should be in the log scale. Either 0 or 1 can be specified.
+  - `**ax_plot_kwargs`: The `kwargs` input to `ax.plot`.
+
+A typical `attainment_ratios` is `[0.25, 0.5, 0.75]` for the uncertainty visualization and `[0.5]` for the representative Pareto front.
+
+## Installation
+
+```shell
+$ pip install optuna matplotlib
+```
+
+## Example
+
+**Example**
+
+```python
+from __future__ import annotations
+
+import matplotlib.pyplot as plt
+import optuna
+import optunahub
+
+
+def objective(trial: optuna.Trial) -> tuple[float, float]:
+    x = trial.suggest_float("x", -5, 5)
+    y = trial.suggest_float("y", -5, 5)
+    return x**2 + y**2, (x - 2)**2 + (y - 2)**2
+
+multiple_study_list = []
+for sampler_cls in [optuna.samplers.RandomSampler, optuna.samplers.TPESampler]:
+    study_list = []
+    for seed in range(10):
+        sampler = sampler_cls(seed=seed)
+        study = optuna.create_study(directions=["minimize", "minimize"], sampler=sampler)
+        study.optimize(objective, n_trials=50)
+        study_list.append(study)
+    multiple_study_list.append(study_list)
+
+plot_multiple_empirical_attainment_surfaces = optunahub.load_module(
+    "visualization/plot_empirical_attainment_surface"
+).plot_multiple_empirical_attainment_surfaces
+ax = plot_multiple_empirical_attainment_surfaces(
+    multiple_study_list, attainment_ratios=[0.25, 0.5, 0.75], labels=["Random", "TPE"]
+)
+plt.show()
+```
+
+## Others
+
+### Bibtex
+
+Please cite [our paper](https://arxiv.org/abs/2305.08852), when using this module.
+
+```bibtex
+@article{watanabe2023python,
+  title={{P}ython tool for visualizing variability of {P}areto fronts over multiple runs},
+  author={Watanabe, Shuhei},
+  journal={arXiv preprint arXiv:2305.08852},
+  year={2023}
+}
+```

package/visualization/plot_empirical_attainment_surface/__init__.py

@@ -0,0 +1,5 @@
+from ._wrapper import plot_empirical_attainment_surface
+from ._wrapper import plot_multiple_empirical_attainment_surfaces
+
+
+__all__ = ["plot_empirical_attainment_surface", "plot_multiple_empirical_attainment_surfaces"]

package/visualization/plot_empirical_attainment_surface/_eaf.py

@@ -0,0 +1,144 @@
+from __future__ import annotations
+
+import numpy as np
+import optuna
+
+
+EPS = 1e-300
+
+
+def _compute_emp_att_surf(
+    X: np.ndarray, pareto_list: list[np.ndarray], levels: np.ndarray
+) -> np.ndarray:
+    """
+    Compute the empirical attainment surface of the given Pareto front sets.
+
+    Args:
+        x:
+            The first objective values appeared in pareto_list.
+            This array is sorted in the ascending order.
+            The shape is (number of possible values, ).
+        levels:
+            A list of `level` described below:
+                Control the k in the k-% attainment surface.
+                    k = level / n_independent_runs
+                must hold.
+                level must be in [1, n_independent_runs].
+                level=1 leads to the best attainment surface,
+                level=n_independent_runs leads to the worst attainment surface,
+                level=n_independent_runs//2 leads to the median attainment surface.
+        pareto_list:
+            The list of the Pareto front sets.
+            The shape is (trial number, Pareto solution index, objective index).
+            Note that each pareto front set is sorted based on the ascending order of
+            the first objective.
+
+    Returns:
+        emp_att_surfs (np.ndarray):
+            The vertices of the empirical attainment surfaces for each level.
+            If emp_att_surf[i, j, 1] takes np.inf, this is not actually on the surface.
+            The shape is (levels.size, X.size, 2).
+
+    Reference:
+        Title: On the Computation of the Empirical Attainment Function
+        Authors: Carlos M. Fonseca et al.
+        https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.705.1929&rep=rep1&type=pdf
+
+    NOTE:
+        Our algorithm is slightly different from the original one, but the result will be same.
+        More details below:
+            When we define N = n_independent_runs, K = X.size, and S = n_samples,
+            the original algorithm requires O(NK + K log K)
+            and this algorithm requires O(NK log K).
+            Although our algorithm is a bit worse than the original algorithm,
+            since the enumerating Pareto solutions requires O(NS log S),
+            which might be smaller complexity but will take more time in Python,
+            the time complexity will not dominate the whole process.
+    """
+    n_levels = len(levels)
+    emp_att_surfs = np.zeros((n_levels, X.size, 2))
+    emp_att_surfs[..., 0] = X
+    n_independent_runs = len(pareto_list)
+    y_candidates = np.zeros((X.size, n_independent_runs))
+    for i, pf_set in enumerate(pareto_list):
+        ub = np.searchsorted(pf_set[:, 0], X, side="right")
+        y_min = np.minimum.accumulate(np.hstack([np.inf, pf_set[:, 1]]))
+        y_candidates[:, i] = y_min[ub]
+    else:
+        y_candidates = np.sort(y_candidates, axis=-1)
+
+    y_sol = y_candidates[:, levels - 1].T
+    emp_att_surfs[..., 1] = y_sol
+
+    for emp_att_surf in emp_att_surfs:
+        idx = np.sum(emp_att_surf[:, 1] == np.inf)
+        emp_att_surf[:idx, 0] = emp_att_surf[idx, 0]
+
+    return emp_att_surfs
+
+
+def _get_empirical_attainment_surfaces(
+    study_list: list[optuna.Study], levels: list[int], log_scale_inds: list[int] | None = None
+) -> np.ndarray:
+    """
+    Get the empirical attainment surfaces given a list of studies.
+
+    Args:
+        study_list:
+            A list of studies to visualize the empirical attainment function.
+        levels:
+            A list of `level` described below:
+                Control the k in the k-% attainment surface.
+                    k = level / n_independent_runs
+                must hold.
+                level must be in [1, n_independent_runs].
+                level=1 leads to the best attainment surface,
+                level=n_independent_runs leads to the worst attainment surface,
+                level=n_independent_runs//2 leads to the median attainment surface.
+        log_scale_inds:
+            The indices of the log scale.
+            For example, if you would like to plot the first objective in the log scale,
+            you need to feed log_scale_inds=[0].
+            In principle, log_scale_inds changes the minimum value of the axes
+            from -np.inf to a small positive value.
+
+    Returns:
+        emp_att_surfs (np.ndarray):
+            The surfaces attained by (level / n_independent_runs) * 100% (for each level in levels)
+            of the trials. In other words, (level / n_independent_runs) * 100% of runs dominate or
+            at least include those solutions in their Pareto front. Note that we only return the
+            Pareto front of attained solutions.
+    """
+    if any(len(study.directions) != 2 for study in study_list):
+        raise NotImplementedError("Three or more objectives are not supported.")
+    if not all(1 <= level <= len(study_list) for level in levels):
+        raise ValueError(
+            f"All elements in levels must be in [1, {len(study_list)=}], but got {levels=}."
+        )
+    if not np.all(np.maximum.accumulate(levels) == levels):
+        raise ValueError(f"levels must be an increasing sequence, but got {levels}.")
+    larger_is_better_objectives = [
+        i
+        for i, d in enumerate(study_list[0].directions)
+        if d == optuna.study.StudyDirection.MAXIMIZE
+    ]
+    pareto_list = []
+    for study in study_list:
+        trials = study._get_trials(deepcopy=False, states=(optuna.trial.TrialState.COMPLETE,))
+        sign_ = np.ones(2, dtype=float)
+        sign_[larger_is_better_objectives] = -1
+        loss_vals = sign_ * np.array([t.values for t in trials])
+        on_front = optuna.study._multi_objective._is_pareto_front(
+            loss_vals, assume_unique_lexsorted=False
+        )
+        pareto_list.append(loss_vals[on_front])
+
+    log_scale_inds = log_scale_inds or []
+    pareto_sols = np.vstack(pareto_list)
+    X = np.unique(np.hstack([EPS if 0 in log_scale_inds else -np.inf, pareto_sols[:, 0], np.inf]))
+    emp_att_surfs = _compute_emp_att_surf(X=X, pareto_list=pareto_list, levels=np.asarray(levels))
+    emp_att_surfs[..., larger_is_better_objectives] *= -1
+    if larger_is_better_objectives is not None and 0 in larger_is_better_objectives:
+        emp_att_surfs = np.flip(emp_att_surfs, axis=1)
+
+    return emp_att_surfs