def prevalence_plot(
  method_names: str | Sequence[str],
  true_prevalences: np.ndarray | Sequence[np.ndarray],
  estim_prevalences: np.ndarray | Sequence[np.ndarray],
  class_name: str = "positive class",
  train_prevalence: float | Sequence[float] = None,
  n_bins: int = 21,
  show_std: bool = True,
  colors: Sequence = None,
  markers: Sequence[str] = None,
  fig_size: tuple[float, float] = (6, 6),
  font_size: int = 11,
  legend_font_size: int = None,
  marker_size: float = 5.0,
  line_width: float = 1.5,
  band_alpha: float = 0.2,
  title: str = "Prevalence Plot",
  grid: bool = True,
  ax: matplotlib.axes.Axes = None,
) -> matplotlib.figure.Figure:
  """Plot predicted vs. true prevalence for a target class (binary diagonal plot).

  For each unique entry in `method_names`, all `(true, estim)` pairs across
  every experiment sharing that name are pooled, binned by true prevalence
  into `n_bins` equal-width intervals, and summarized by the mean predicted
  prevalence per bin (connected by a line) with an optional +/- 1 standard
  deviation shaded band. A dashed diagonal (`y = x`) marks the ideal,
  unbiased quantifier.

  Parameters
  ----------
  method_names: str | Sequence[str]
    Name of the method for each experiment. A name can repeat (e.g. one entry per dataset/fold); all
    matching experiments are merged before binning.
  true_prevalences: np.ndarray | Sequence[np.ndarray]
    True prevalence of `class_name` for each test bag, one 1D array per experiment
    (aligned with `method_names`).
  estim_prevalences: np.ndarray | Sequence[np.ndarray]
    Predicted prevalence of `class_name`, same shape as `true_prevalences`.
  class_name: str, default = "positive class"
    Label used on the axes/legend for the target class. Defaults to "positive class".
  train_prevalence: float | Sequence[float], default = None
    One or more training prevalences to mark on the diagonal. Defaults to None.
  n_bins: int, default = 21
    Number of equal-width bins over `[0, 1]` used to aggregate repeated experiments. Defaults to 21.
  show_st: bool, default = True
    Whether to draw +/- 1 std shaded bands around each method's line. Defaults to True.
  colors: Sequence, default = colorblind-sage palette
    Custom colors, one per unique method. Defaults to a colorblind-safe palette, auto-extended as needed.
  markers: Sequence[str], default = None
    Custom marker styles, one per unique method. Defaults to a built-in marker cycle.
  fig_size: tuple[float, float], default = (6, 6)
    Figure size in inches. Defaults to (6, 6).
  font_size: int, default = 11
    Base font size for axis labels/title. Defaults to 11.
  legend_font_size: int, default = None
    Font size for the legend. Defaults to `font_size - 1` when None.
  marker_size: float, default = 5.0
    Marker size. Defaults to 5.0.
  line_width: float, default = 1.5
    Line width. Defaults to 1.5.
  band_alpha: float, default = 0.2
    Opacity of the +/- 1 std band. Defaults to 0.2.
  title: str, default = "Prevalence Plot"
    Plot title. Defaults to "Prevalence Plot".
  grid: bool, default = True
    Whether to draw a background grid. Defaults to True.
  ax: matplotlib.axes.Axes, default = None
    Existing axes to draw on. A new figure/axes pair is created when None. Defaults to None.

  Returns
  -------
  fig: matplotlib.figure.Figure
    The generated figure. Call `fig.savefig(path)` (png, pdf, svg, ... — any Matplotlib-supported format)
    to persist it.

  Examples
  --------
  >>> import numpy as np
  >>> from quack.visualization import prevalence_plot
  >>> rng = np.random.default_rng(0)
  >>> true_prev = rng.uniform(0, 1, 200)
  >>> estim_prev = np.clip(true_prev + rng.normal(0, 0.05, 200), 0, 1)
  >>> fig = prevalence_plot("CC", true_prev, estim_prev, train_prevalence=0.5)
  >>> fig.savefig("prevalence.png", dpi=300)
  """
  method_names = [method_names] if isinstance(method_names, str) else list(method_names)
  true_prevalences = _normalize_experiments(true_prevalences)
  estim_prevalences = _normalize_experiments(estim_prevalences)

  if not (len(method_names) == len(true_prevalences) == len(estim_prevalences)):
    raise ValueError(
      "method_names, true_prevalences and estim_prevalences must have the "
      "same length (one entry per experiment)."
    )

  unique_methods = list(dict.fromkeys(method_names))  # preserves first-seen order
  colors = get_color_palette(len(unique_methods), palette=colors)
  markers = get_marker_cycle(len(unique_methods), markers=markers)
  legend_font_size = legend_font_size if legend_font_size is not None else max(font_size - 1, 6)

  bin_edges = np.linspace(0.0, 1.0, n_bins + 1)
  bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

  own_axes = ax is None
  if own_axes:
    fig, ax = plt.subplots(figsize=fig_size)
  else:
    fig = ax.get_figure()

  ax.set_aspect('equal')
  ax.plot([0, 1], [0, 1], color=REFERENCE_COLOR, linestyle='--', linewidth=1,
          label='Ideal quantifier (y=x)', zorder=1)

  for idx, method in enumerate(unique_methods):
    true_pool = np.concatenate([
      np.atleast_1d(true_prevalences[i]) for i, name in enumerate(method_names) if name == method
    ])
    estim_pool = np.concatenate([
      np.atleast_1d(estim_prevalences[i]) for i, name in enumerate(method_names) if name == method
    ])

    bin_idx = np.clip(np.digitize(true_pool, bin_edges[1:-1]), 0, n_bins - 1)
    means = np.full(n_bins, np.nan)
    stds = np.full(n_bins, np.nan)
    for b in range(n_bins):
      values = estim_pool[bin_idx == b]
      if values.size > 0:
        means[b] = values.mean()
        stds[b] = values.std() if values.size > 1 else 0.0

    valid = ~np.isnan(means)
    color = colors[idx]
    ax.plot(bin_centers[valid], means[valid], label=method, color=color,
            marker=markers[idx], markersize=marker_size, linewidth=line_width,
            zorder=3)

    if show_std:
      ax.fill_between(bin_centers[valid], np.clip(means[valid] - stds[valid], 0, 1),
                       np.clip(means[valid] + stds[valid], 0, 1), color=color,
                       alpha=band_alpha, zorder=2, linewidth=0)

  if train_prevalence is not None:
    train_prevalence = (train_prevalence if isinstance(train_prevalence, (list, tuple, np.ndarray))
                        else [train_prevalence])
    for p in train_prevalence:
      ax.scatter(p, p, s=(marker_size * 9), color='red', edgecolors=REFERENCE_COLOR,
                 marker='*', linewidth=1.2, zorder=5,
                 label=f'Training prevalence (p={p:.2f})')

  ax.set_xlim(0, 1)
  ax.set_ylim(0, 1)
  ax.set_xlabel(f"True prevalence ({class_name})", fontsize=font_size)
  ax.set_ylabel(f"Estimated prevalence ({class_name})", fontsize=font_size)
  if title:
    ax.set_title(title, fontsize=font_size + 2)
  if grid:
    ax.grid(alpha=0.3)

  handles, labels = ax.get_legend_handles_labels()
  by_label = dict(zip(labels, handles))  # de-duplicate repeated train-prevalence labels
  ax.legend(by_label.values(), by_label.keys(), fontsize=legend_font_size,
            loc='upper center', bbox_to_anchor=(0.5, -0.15),
            ncol=min(3, len(by_label)), frameon=False)

  if own_axes:
    fig.tight_layout()

  return fig

def bias_plot(
  method_names: str | Sequence[str],
  true_prevalences: np.ndarray | Sequence[np.ndarray],
  estim_prevalences: np.ndarray | Sequence[np.ndarray],
  class_name: str = "positive class",
  n_bins: int = 1,
  colors: Sequence = None,
  fig_size: tuple[float, float] = (8, 6),
  font_size: int = 11,
  legend_font_size: int = None,
  title: str = "Bias Plot",
  box_width: float = 0.6,
  grid: bool = True,
  ax: matplotlib.axes.Axes = None,
) -> matplotlib.figure.Figure:
  """Plot the distribution of signed prevalence errors per method as box plots.

  The bias for a single test bag is defined as
  `bias = estimated_prevalence - true_prevalence` for `class_name`. A value
  of 0 indicates a perfectly unbiased estimate; positive values indicate a
  tendency to overestimate the class, negative values a tendency to
  underestimate it.

  When `n_bins > 1`, the true test prevalence range `[0, 1]` is split into
  `n_bins` equal-width intervals and one group of boxes (one box per
  method) is drawn per interval, mirroring QuaPy's `binary_bias_bins`, so
  that bias-vs-prevalence patterns invisible in the global view can be
  detected.

  Parameters
  ----------
  method_names: str | Sequence[str]
    Name of the method for each experiment (can repeat across datasets/folds;
    matching experiments are pooled together).
  true_prevalences: np.ndarray | Sequence[np.ndarray]
    True prevalence of `class_name` per test bag, one 1D array per experiment.
  estim_prevalences: np.ndarray | Sequence[np.ndarray]
    Predicted prevalence of `class_name`, same shape as `true_prevalences`.
  class_name: str, default = "positive class"
    Target class label used in axis/legend text. Defaults to "positive class".
  n_bins: int, default = 1
    Number of equal-width true-prevalence bins. Use `1` for a single global box
    per method; use `>1` to break it down by true-prevalence range.
    Defaults to 1.
  colors: Sequence, default = colorblind-safe palette
    Custom colors, one per unique method. Defaults to a colorblind-safe palette,
    auto-extended as needed.
  fig_size: tuple[float, float], default = (8, 6)
    Figure size in inches. Defaults to (8, 6).
  font_size: int, default = 11
    Base font size for axis labels/title. Defaults to 11.
  legend_font_size: int, default = None
    Legend font size. Defaults to `font_size - 1` when None.
  title: str, default = "Bias Plot"
    Plot title. Defaults to "Bias Plot".
  box_width: float, default = 0.6
    Width of each individual box. Defaults to 0.6.
  grid: bool, default = True
    Whether to draw a background grid. Defaults to True.
  ax: matplotlib.axes.Axes, default = None
    Existing axes to draw on. A new figure/axes pair is created
    when None. Defaults to None.

  Returns
  -------
  fig: matplotlib.figure.Figure
    The generated figure. Call `fig.savefig(path)` to persist it in any
    Matplotlib-supported format.

  Examples
  --------
  >>> import numpy as np
  >>> from quack.visualization import quantification_bias_plot
  >>> rng = np.random.default_rng(0)
  >>> true_prev = rng.uniform(0, 1, 200)
  >>> estim_prev = np.clip(true_prev + 0.1 + rng.normal(0, 0.05, 200), 0, 1)
  >>> fig = quantification_bias_plot("CC", true_prev, estim_prev, n_bins=3)
  """
  method_names = [method_names] if isinstance(method_names, str) else list(method_names)
  true_prevalences = _normalize_experiments(true_prevalences)
  estim_prevalences = _normalize_experiments(estim_prevalences)

  if not (len(method_names) == len(true_prevalences) == len(estim_prevalences)):
    raise ValueError(
      "method_names, true_prevalences and estim_prevalences must have the "
      "same length (one entry per experiment)."
    )

  unique_methods = list(dict.fromkeys(method_names))
  n_methods = len(unique_methods)
  colors = get_color_palette(n_methods, palette=colors)
  legend_font_size = legend_font_size if legend_font_size is not None else max(font_size - 1, 6)

  pooled_true, pooled_bias = {}, {}
  for method in unique_methods:
    true_pool = np.concatenate([
      np.atleast_1d(true_prevalences[i]) for i, name in enumerate(method_names) if name == method
    ])
    estim_pool = np.concatenate([
      np.atleast_1d(estim_prevalences[i]) for i, name in enumerate(method_names) if name == method
    ])
    pooled_true[method] = true_pool
    pooled_bias[method] = estim_pool - true_pool

  own_axes = ax is None
  if own_axes:
    fig, ax = plt.subplots(figsize=fig_size)
  else:
    fig = ax.get_figure()

  bin_edges = np.linspace(0.0, 1.0, n_bins + 1)
  intra_gap = box_width * 1.15  # spacing between method boxes within a prevalence-bin group

  legend_handles = []
  for m_idx, method in enumerate(unique_methods):
    offset = (m_idx - (n_methods - 1) / 2) * intra_gap
    box_data, positions = [], []

    for b in range(n_bins):
      lo, hi = bin_edges[b], bin_edges[b + 1]
      in_bin = ((pooled_true[method] >= lo) &
                (pooled_true[method] <= hi if b == n_bins - 1 else pooled_true[method] < hi))
      values = pooled_bias[method][in_bin]
      if values.size > 0:
        box_data.append(values)
        positions.append(b + offset)

    if box_data:
      bp = ax.boxplot(box_data, positions=positions, widths=box_width,
                       patch_artist=True, showfliers=False, manage_ticks=False)
      for patch in bp['boxes']:
        patch.set_facecolor(colors[m_idx])
        patch.set_alpha(0.75)
        patch.set_edgecolor(REFERENCE_COLOR)
      for element in ('whiskers', 'caps', 'medians'):
        for artist in bp[element]:
          artist.set_color(REFERENCE_COLOR)
      legend_handles.append(plt.Line2D([0], [0], marker='s', linestyle='',
                                        markerfacecolor=colors[m_idx],
                                        markeredgecolor=REFERENCE_COLOR,
                                        markersize=10, label=method))

  ax.axhline(0.0, color=REFERENCE_COLOR, linestyle='--', linewidth=1, zorder=0)

  if n_bins == 1:
    ax.set_xticks([0])
    ax.set_xticklabels([""])
    ax.set_xlabel(f"Method (target: {class_name})", fontsize=font_size)
  else:
    tick_positions = list(range(n_bins))
    tick_labels = [f"[{bin_edges[b]:.2f}, {bin_edges[b+1]:.2f}]" for b in range(n_bins)]
    ax.set_xticks(tick_positions)
    ax.set_xticklabels(tick_labels, rotation=30, ha='right')
    ax.set_xlabel(f"True prevalence range ({class_name})", fontsize=font_size)

  ax.set_ylabel(f"Bias (estimated - true) for {class_name}", fontsize=font_size)
  if title:
    ax.set_title(title, fontsize=font_size + 2)
  if grid:
    ax.grid(axis='y', alpha=0.3)

  ref_handle = plt.Line2D([0], [0], color=REFERENCE_COLOR, linestyle='--', label='Unbiased (bias = 0)')
  ax.legend(handles=legend_handles + [ref_handle], fontsize=legend_font_size,
            loc='upper center', bbox_to_anchor=(0.5, -0.2),
            ncol=min(4, n_methods + 1), frameon=False)

  if own_axes:
    fig.tight_layout()

  return fig

def class_distribution_plot(
  y: np.ndarray,
  normalize: bool = False,
  title: str = "Class Distribution",
  fig_size: tuple[float, float] = (8, 5),
  colors: Sequence = None,
  font_size: int = 11,
  bar_label_font_size: int = None,
  show_bar_labels: bool = True,
  horizontal: bool = False,
  grid: bool = True,
  ax: matplotlib.axes.Axes = None,
) -> matplotlib.figure.Figure:
  """Plot the (optionally normalized) class distribution of a label array.

  Parameters
  ----------
  y: np.ndarray
    Array-like with all labels.
  normalize: bool, default = False
    If True, plot relative frequencies (prevalences summing to 1.0)
    instead of raw counts. Defaults to False.
  title: str, default = "Class Distribution"
    Plot title. Defaults to "Class Distribution".
  fig_size: tuple[float, float], default = (8, 5)
    Figure size in inches. Defaults to (8, 5).
  colors: Sequence, default = colorblind-safe palette
    Custom colors, one per class. Defaults to a colorblind-safe palette,
    auto-extended for many classes.
  font_size: int, default = 11
    Base font size for axis labels/title. Defaults to 11.
  bar_label_font_size: int, default = None
    Font size for the value labels drawn on top of each bar.
    Defaults to `font_size - 1` when None.
  show_bar_labels: bool, default = True
    Whether to annotate each bar with its value. Defaults to True.
  horizontal: bool, default = False
    If True, draws horizontal bars
    (useful for many classes / long class names). Defaults to False.
  grid: bool, default =True
    Whether to draw a background grid. Defaults to True.
  ax: matplotlib.axes.Axes, default = None
    Existing axes to draw on. A new figure/axes pair is
    created when None. Defaults to None.

  Returns
  -------
  fig: matplotlib.figure.Figure
    The generated figure. Call `fig.savefig(path)` to persist it in any
    Matplotlib-supported format.

  Examples
  --------
  >>> import numpy as np
  >>> from quack.visualization import class_distribution_plot
  >>> y = np.random.choice([0, 1, 2], size=500, p=[0.6, 0.3, 0.1])
  >>> fig = class_distribution_plot(y, normalize=True)
  >>> fig.savefig("class_distribution.pdf")
  """
  bar_label_font_size = bar_label_font_size if bar_label_font_size is not None else max(font_size - 1, 6)

  series = pd.Series(np.asarray(y), name="class")
  counts = series.value_counts().sort_index()
  values = (counts / counts.sum()) if normalize else counts

  classes = [str(c) for c in values.index]
  colors = get_color_palette(len(classes), palette=colors)

  own_axes = ax is None
  if own_axes:
    fig, ax = plt.subplots(figsize=fig_size)
  else:
    fig = ax.get_figure()

  if horizontal:
    bars = ax.barh(classes, values.values, color=colors, edgecolor=REFERENCE_COLOR, linewidth=0.5)
    ax.set_xlabel("Proportion" if normalize else "# instances", fontsize=font_size)
    ax.set_ylabel("Class", fontsize=font_size)
  else:
    bars = ax.bar(classes, values.values, color=colors, edgecolor=REFERENCE_COLOR, linewidth=0.5)
    ax.set_ylabel("Proportion" if normalize else "# instances", fontsize=font_size)
    ax.set_xlabel("Class", fontsize=font_size)

  if show_bar_labels:
    fmt = "{:.3f}" if normalize else "{:.0f}"
    ax.bar_label(bars, labels=[fmt.format(v) for v in values.values],
                 fontsize=bar_label_font_size, padding=2)

  if title:
    ax.set_title(title, fontsize=font_size + 2)
  if grid:
    ax.grid(axis='x' if horizontal else 'y', alpha=0.3)

  if own_axes:
    fig.tight_layout()

  return fig

def prevalence_coverage_plot(
  prevalences: np.ndarray | Sequence[np.ndarray],
  labels: str | Sequence[str] = None,
  class_name: str = "positive class",
  train_prevalence: float | Sequence[float] = None,
  n_bins: int = 20,
  show_rug: bool = True,
  show_stats: bool = True,
  density: bool = False,
  colors: Sequence = None,
  fig_size: tuple[float, float] = (8, 5),
  font_size: int = 11,
  legend_font_size: int = None,
  bar_alpha: float = 0.65,
  rug_height: float = 0.04,
  title: str = "Prevalence Coverage",
  grid: bool = True,
  ax: matplotlib.axes.Axes = None,
) -> matplotlib.figure.Figure:
  """Plot how well a set of bags covers the `[0, 1]` prevalence range for a class.

  Each entry in `prevalences` is a 1D array with the realized prevalence
  of `class_name` for every generated bag in one experiment (typically
  `generator.sampled_prevalences_[:, class_index]` from a
  `quack.bag_generator.BaseBagGenerator` subclass). A histogram over
  `n_bins` equal-width bins shows how many bags fall in each prevalence
  range, while an optional rug plot marks every individual bag along the
  x-axis, so isolated or empty regions of the simplex are easy to spot
  even when the histogram bin is technically non-empty.

  Multiple series can be overlaid (e.g. to compare `PriorShiftBagGenerator`
  vs. `CovariateShiftBagGenerator`, or different `sampling_strategy`/
  `dirichlet_alpha` configurations) using semi-transparent, colorblind-safe
  colors.

  Parameters
  ----------
  prevalences: np.ndarray | Sequence[np.ndarray]
    One 1D array of per-bag prevalences for `class_name` per experiment/series. A
    single array plots one series.
  labels: str | Sequence[str], default = None
    Name of each series, used in the legend and coverage statistics. Defaults to
    `"Bags"` for a single series, or `"Series 1"`, `"Series 2"`, ... for multiple.
  class_name: str, default = "positive class"
    Label used on the x-axis for the target class. Defaults to "positive class".
  train_prevalence: float | Sequence[float], default = None
    One or more training prevalences to mark as vertical reference lines.
    Defaults to None.
  n_bins: int, default = 20
    Number of equal-width histogram bins over `[0, 1]`. Defaults to 20.
  show_rug: bool, default = True
    Whether to draw a rug plot (one tick per bag) below the histogram.
    Defaults to True.
  show_stats: bool, default = True
    Whether to annotate the plot with per-series min/max/mean prevalence
    and simplex-bin coverage (the fraction of the `n_bins` bins that contain
    at least one bag). Defaults to True.
  density: bool, default = False
    If True, normalize histograms to a density (area sums to 1) instead of raw
    bag counts — useful when comparing series generated with a different number
    of bags. Defaults to False.
  colors: Sequence, default = None 
    Custom colors, one per series. Defaults to a colorblind-safe palette,
    auto-extended as needed.
  fig_size: tuple[float, float], default = (8, 6)
    Figure size in inches. Defaults to (8, 5).
  font_size: int, default = 11
    Base font size for axis labels/title. Defaults to 11.
  legend_font_size: int, default = None
    Legend font size. Defaults to `font_size - 1` when None.
  bar_alpha: float, default = 0.65
    Opacity of the histogram bars, low enough for overlapping
    series to remain distinguishable. Defaults to 0.65.
  rug_height: float, default = 0.04
    Height of each rug row, as a fraction of the histogram's y-range.
    Defaults to 0.04.
  title: str, default = "Prevalence Coverage"
    Plot title. Defaults to "Prevalence Coverage".
  grid: bool, default = True
    Whether to draw a background grid. Defaults to True.
  ax: matplotlib.axes.Axes, default = None
    Existing axes to draw on. A new figure/axes pair is created
    when None. Defaults to None.

  Returns
  -------
  fig: matplotlib.figure.Figure
    The generated figure. Call `fig.savefig(path)` to persist it
    in any Matplotlib-supported format.

  Raises
  ------
  ValueError: If `labels` length does not match the number of series,
    any series is not 1D, or any prevalence value falls outside
    `[0, 1]`.

  Examples
  --------
  >>> import numpy as np
  >>> from quack.bag_generator import PriorShiftBagGenerator, CovariateShiftBagGenerator
  >>> from quack.visualization import prevalence_coverage_plot
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=500, n_classes=2, random_state=0)
  >>> prior_gen = PriorShiftBagGenerator(n_bags=200, bag_size=100, random_state=0)
  >>> cov_gen = CovariateShiftBagGenerator(n_bags=200, bag_size=100, random_state=0)
  >>> prior_gen.to_list(X, y)
  >>> cov_gen.to_list(X, y)
  >>> fig = prevalence_coverage_plot(
  ...     [prior_gen.sampled_prevalences_[:, 1], cov_gen.sampled_prevalences_[:, 1]],
  ...     labels=["Prior Shift", "Covariate Shift"],
  ...     class_name="positive class",
  ... )
  >>> fig.savefig("prevalence_coverage.png", dpi=300)
  """
  prevalences = _normalize_experiments(prevalences)
  n_series = len(prevalences)

  if labels is None:
    labels = [f"Series {i + 1}" for i in range(n_series)] if n_series > 1 else ["Bags"]
  else:
    labels = [labels] if isinstance(labels, str) else list(labels)

  if len(labels) != n_series:
    raise ValueError(
      f"labels must have the same length as prevalences ({n_series}), got {len(labels)}."
    )

  clean_prevalences = []
  for label, p in zip(labels, prevalences):
    p = np.asarray(p, dtype=float)
    if p.ndim != 1:
      raise ValueError(f"Each prevalence array must be 1D, got shape {p.shape} for series '{label}'.")
    if p.size and np.any((p < 0) | (p > 1)):
      raise ValueError(f"Prevalence values must lie within [0, 1]; series '{label}' violates this.")
    clean_prevalences.append(p)
  prevalences = clean_prevalences

  colors = get_color_palette(n_series, palette=colors)
  legend_font_size = legend_font_size if legend_font_size is not None else max(font_size - 1, 6)

  own_axes = ax is None
  if own_axes:
    fig, ax = plt.subplots(figsize=fig_size)
  else:
    fig = ax.get_figure()

  bin_edges = np.linspace(0.0, 1.0, n_bins + 1)

  for i, p in enumerate(prevalences):
    ax.hist(p, bins=bin_edges, density=density, color=colors[i], alpha=bar_alpha,
            edgecolor=REFERENCE_COLOR, linewidth=0.5, label=labels[i], zorder=2)

  if show_rug:
    y_min, y_max = ax.get_ylim()
    rug_span = rug_height * (y_max - y_min) if y_max > y_min else rug_height
    for i, p in enumerate(prevalences):
      row_offset = -rug_span * (i + 1) * 1.3
      ax.plot(p, np.full_like(p, y_min + row_offset), '|', color=colors[i],
              markersize=8, markeredgewidth=1.2, alpha=0.85, zorder=3)
    ax.set_ylim(y_min - rug_span * 1.3 * (n_series + 1), y_max)

  if train_prevalence is not None:
    train_prevalence = (train_prevalence if isinstance(train_prevalence, (list, tuple, np.ndarray))
                        else [train_prevalence])
    for p in train_prevalence:
      ax.axvline(p, color=REFERENCE_COLOR, linestyle='--', linewidth=1.3, zorder=4,
                 label=f'Training prevalence (p={p:.2f})')

  ax.set_xlim(0, 1)
  ax.set_xlabel(f"Prevalence ({class_name})", fontsize=font_size)
  ax.set_ylabel("Density" if density else "# bags", fontsize=font_size)
  if title:
    ax.set_title(title, fontsize=font_size + 2)
  if grid:
    ax.grid(axis='y', alpha=0.3)

  if show_stats:
    stats_lines = []
    for label, p in zip(labels, prevalences):
      if p.size == 0:
        stats_lines.append(f"{label}: no bags")
        continue
      bin_idx = np.clip(np.digitize(p, bin_edges[1:-1]), 0, n_bins - 1)
      coverage = np.unique(bin_idx).size / n_bins
      stats_lines.append(
        f"{label}: min={p.min():.2f} max={p.max():.2f} mean={p.mean():.2f} coverage={coverage:.0%}"
      )
    ax.text(0.01, 0.98, "\n".join(stats_lines), transform=ax.transAxes,
            fontsize=max(font_size - 2, 6), va='top', ha='left', color=REFERENCE_COLOR,
            bbox=dict(boxstyle='round', facecolor='white', edgecolor=REFERENCE_COLOR, alpha=0.8))

  handles, hlabels = ax.get_legend_handles_labels()
  by_label = dict(zip(hlabels, handles))  # de-duplicate repeated train-prevalence labels
  ax.legend(by_label.values(), by_label.keys(), fontsize=legend_font_size,
            loc='upper center', bbox_to_anchor=(0.5, -0.15),
            ncol=min(3, len(by_label)), frameon=False)

  if own_axes:
    fig.tight_layout()

  return fig

def get_color_palette(n_colors: int, palette: Sequence = None) -> list:
  """Build a list of `n_colors` visually distinct colors.

  Falls back to the colorblind-safe base palette while there are enough
  colors available. When more colors than the base palette are requested
  (e.g. many quantifiers being compared at once), it extends the palette
  by uniformly sampling a perceptually-uniform colormap so that all
  colors remain distinguishable from each other.

  Parameters
  ----------
  n_colors: int
    Number of distinct colors needed.
  palette: Sequence, default = None
    User-provided palette (hex strings or RGBA tuples) to use instead of
    the default colorblind-safe one.

  Returns
  -------
  colors_list: list
    List of length `n_colors` containing hex strings or RGBA tuples.
  """
  base = list(palette) if palette is not None else list(COLORBLIND_PALETTE)

  if n_colors <= len(base):
    return base[:n_colors]

  # extend deterministically using a perceptually-uniform colormap so
  # additional colors stay maximally separated from one another
  cmap = plt.get_cmap('turbo')
  n_extra = n_colors - len(base)
  extra = [cmap(x) for x in np.linspace(0.05, 0.95, n_extra)]
  return base + extra

def get_marker_cycle(n_markers: int, markers: Sequence[str] = None) -> list:
  """Cycle through a fixed list of distinguishable marker shapes.

  Parameters
  ----------
  n_markers: int
    Number of markers needed.
  markers: Sequence[str], default = None
    Custom marker list.

  Returns
  -------
  markers_shapes: list
    List of length `n_markers` with matplotlib marker style strings.
  """
  base = list(markers) if markers is not None else list(MARKERS)
  return [base[i % len(base)] for i in range(n_markers)]