canvod.audit API Reference¶

def compare_datasets(
    ds_a: Any,
    ds_b: Any,
    *,
    variables: list[str] | None = None,
    tier: ToleranceTier = ToleranceTier.NUMERICAL,
    tolerance_overrides: dict[str, Tolerance] | None = None,
    label: str = "",
    align: bool = True,
    metadata: dict[str, Any] | None = None,
    report_coverage: bool = False,
) -> ComparisonResult:
    """Compare two xarray Datasets variable-by-variable.

    Parameters
    ----------
    ds_a, ds_b : xarray.Dataset
        Datasets to compare. ``ds_a`` is typically the candidate (canvodpy),
        ``ds_b`` the reference.
    variables : list[str], optional
        Variables to compare. If None, uses the intersection of both datasets'
        data variables.
    tier : ToleranceTier
        Comparison strictness level.
    tolerance_overrides : dict[str, Tolerance], optional
        Per-variable tolerance overrides.
    label : str
        Human-readable label for this comparison.
    align : bool
        If True, align datasets on shared (epoch, sid) coordinates.
    metadata : dict, optional
        Free-form metadata to attach to the result.

    Returns
    -------
    ComparisonResult
    """
    alignment = None
    if align:
        ds_a, ds_b, alignment = _align_datasets(ds_a, ds_b)

        # Guard: fail if alignment produced zero overlap
        if alignment.n_shared_epochs == 0 or alignment.n_shared_sids == 0:
            return ComparisonResult(
                label=label or f"{tier.value} comparison",
                variable_stats={},
                tier=tier,
                passed=False,
                failures={"_alignment": "No shared coordinates after alignment"},
                metadata=metadata or {},
                alignment=alignment,
            )
        # Warn when >5 % of the union is dropped — significant data loss
        total_epochs = (
            alignment.n_shared_epochs
            + alignment.n_dropped_epochs_a
            + alignment.n_dropped_epochs_b
        )
        total_sids = (
            alignment.n_shared_sids
            + alignment.n_dropped_sids_a
            + alignment.n_dropped_sids_b
        )
        drop_epoch_pct = (
            (alignment.n_dropped_epochs_a + alignment.n_dropped_epochs_b) / total_epochs
            if total_epochs > 0
            else 0.0
        )
        drop_sid_pct = (
            (alignment.n_dropped_sids_a + alignment.n_dropped_sids_b) / total_sids
            if total_sids > 0
            else 0.0
        )
        if drop_epoch_pct > 0.05 or drop_sid_pct > 0.05:
            warnings.warn(
                f"[{label or 'compare_datasets'}] Alignment dropped "
                f"{drop_epoch_pct:.1%} of epochs "
                f"({alignment.n_dropped_epochs_a} from A, "
                f"{alignment.n_dropped_epochs_b} from B) and "
                f"{drop_sid_pct:.1%} of SIDs "
                f"({alignment.n_dropped_sids_a} from A, "
                f"{alignment.n_dropped_sids_b} from B). "
                "Results cover only the intersection.",
                stacklevel=2,
            )

    # Determine variables to compare
    vars_a = set(ds_a.data_vars)
    vars_b = set(ds_b.data_vars)
    if variables is None:
        variables = sorted(vars_a & vars_b)

    # Compute per-variable stats
    variable_stats: dict[str, VariableStats] = {}
    failures: dict[str, str] = {}
    cov_valid_both: dict[str, int] = {}
    cov_valid_a_only: dict[str, int] = {}
    cov_valid_b_only: dict[str, int] = {}
    cov_neither: dict[str, int] = {}

    for var in variables:
        if var not in ds_a.data_vars or var not in ds_b.data_vars:
            continue

        try:
            a_vals = ds_a[var].values.astype(np.float64)
            b_vals = ds_b[var].values.astype(np.float64)
        except ValueError, TypeError:
            failures[var] = "non-numeric dtype, skipped"
            continue

        if a_vals.shape != b_vals.shape:
            failures[var] = f"shape mismatch: {a_vals.shape} vs {b_vals.shape}"
            continue

        # Strip attrs from the DataArray and all its coordinates before
        # comparing — da.equals() already ignores attrs per xarray spec, but
        # this makes the contract explicit.  Coordinate VALUES and dtypes are
        # still checked; if they differ between stores that is meaningful.
        exact_match = bool(_strip_attrs(ds_a[var]).equals(_strip_attrs(ds_b[var])))

        vs = compute_variable_stats(var, a_vals, b_vals, exact_match=exact_match)
        variable_stats[var] = vs

        # Tolerance check is annotation-only — never overrides exact_match for `passed`
        tol = get_tolerance(var, tier, tolerance_overrides)
        reason = _check_tolerance(vs, tol)
        if reason:
            failures[var] = reason

        # Structural diagnostic: values are identical but exact_match=False —
        # use xr.testing.assert_equal to pinpoint what structural element
        # differs (coord dtype, missing coord, dimension order, etc.).
        # Only runs when no numeric or NaN-rate failure was already recorded,
        # so it doesn't shadow a real science difference.
        if not exact_match and vs.max_abs_diff == 0.0 and var not in failures:
            msg = _structural_diff(ds_a[var], ds_b[var])
            failures[var] = f"structural (values identical): {msg}"

        if report_coverage:
            mask_a = ~np.isnan(a_vals)
            mask_b = ~np.isnan(b_vals)
            cov_valid_both[var] = int(np.sum(mask_a & mask_b))
            cov_valid_a_only[var] = int(np.sum(mask_a & ~mask_b))
            cov_valid_b_only[var] = int(np.sum(~mask_a & mask_b))
            cov_neither[var] = int(np.sum(~mask_a & ~mask_b))

    coverage = (
        CoverageReport(
            vars_a_only=sorted(vars_a - vars_b),
            vars_b_only=sorted(vars_b - vars_a),
            valid_both=cov_valid_both,
            valid_a_only=cov_valid_a_only,
            valid_b_only=cov_valid_b_only,
            neither_valid=cov_neither,
        )
        if report_coverage
        else None
    )

    # passed iff every compared variable is within tolerance for the given tier.
    # For EXACT tier (atol=0): a variable with any diff is in failures, so
    # len(failures)==0 implies bit-identical.  For SCIENTIFIC tier: implies
    # within stated tolerances.  exact_match is still recorded per-variable for
    # informational purposes (shown in the Typst ✓/✗ column).
    passed = bool(variable_stats) and len(failures) == 0

    return ComparisonResult(
        label=label or f"{tier.value} comparison",
        variable_stats=variable_stats,
        tier=tier,
        passed=passed,
        failures=failures,
        metadata=metadata or {},
        alignment=alignment,
        coverage=coverage,
    )

def compute_variable_stats(
    name: str,
    a: np.ndarray,
    b: np.ndarray,
    exact_match: bool | None = None,
) -> VariableStats:
    """Compute all comparison statistics for a single variable.

    Parameters
    ----------
    name : str
        Variable name.
    a, b : np.ndarray
        Arrays of equal shape to compare.
    exact_match : bool, optional
        Pre-computed exact equality (e.g. from ``da.equals()``).
        If None, inferred as ``max_abs_diff == 0``.
    """
    a = a.ravel().astype(np.float64)
    b = b.ravel().astype(np.float64)

    n_total = len(a)
    mask_a = np.isfinite(a)
    mask_b = np.isfinite(b)
    mask = mask_a & mask_b

    n_nan_a = int(np.isnan(a).sum())
    n_nan_b = int(np.isnan(b).sum())
    n_compared = int(mask.sum())
    n_a_only = int(np.sum(mask_a & ~mask_b))
    n_b_only = int(np.sum(~mask_a & mask_b))
    n_neither = int(np.sum(~mask_a & ~mask_b))

    if n_compared == 0:
        _nan = float("nan")
        return VariableStats(
            name=name,
            exact_match=exact_match if exact_match is not None else True,
            n_total=n_total,
            n_compared=0,
            n_a_only=n_a_only,
            n_b_only=n_b_only,
            n_neither=n_neither,
            n_nonzero_diff=0,
            max_abs_diff=_nan,
            rmse=_nan,
            mae=_nan,
            bias=_nan,
            correlation=_nan,
            p50=_nan,
            p90=_nan,
            p99=_nan,
            n_nan_a=n_nan_a,
            n_nan_b=n_nan_b,
            pct_nan_a=n_nan_a / n_total if n_total > 0 else 0.0,
            pct_nan_b=n_nan_b / n_total if n_total > 0 else 0.0,
            nan_agreement_rate=nan_agreement(a, b),
        )

    # Compute diff once, reuse for all stats
    diff = a[mask] - b[mask]
    absdiff = np.abs(diff)
    n_nonzero_diff = int(np.sum(absdiff > 0))

    _max_abs = float(np.max(absdiff))
    if exact_match is None:
        exact_match = _max_abs == 0.0

    return VariableStats(
        name=name,
        exact_match=exact_match,
        n_total=n_total,
        n_compared=n_compared,
        n_a_only=n_a_only,
        n_b_only=n_b_only,
        n_neither=n_neither,
        n_nonzero_diff=n_nonzero_diff,
        max_abs_diff=_max_abs,
        rmse=float(np.sqrt(np.mean(diff**2))),
        mae=float(np.mean(absdiff)),
        bias=float(np.mean(diff)),
        correlation=correlation(a, b),
        p50=float(np.percentile(absdiff, 50)),
        p90=float(np.percentile(absdiff, 90)),
        p99=float(np.percentile(absdiff, 99)),
        n_nan_a=n_nan_a,
        n_nan_b=n_nan_b,
        pct_nan_a=n_nan_a / n_total if n_total > 0 else 0.0,
        pct_nan_b=n_nan_b / n_total if n_total > 0 else 0.0,
        nan_agreement_rate=nan_agreement(a, b),
    )

def rmse(a: np.ndarray, b: np.ndarray) -> float:
    """Root Mean Square Error over mutually valid elements."""
    mask = _valid_mask(a, b)
    if not mask.any():
        return float("nan")
    diff = a[mask] - b[mask]
    return float(np.sqrt(np.mean(diff**2)))

def bias(a: np.ndarray, b: np.ndarray) -> float:
    """Mean difference (a - b) over mutually valid elements."""
    mask = _valid_mask(a, b)
    if not mask.any():
        return float("nan")
    return float(np.mean(a[mask] - b[mask]))

def mae(a: np.ndarray, b: np.ndarray) -> float:
    """Mean Absolute Error over mutually valid elements."""
    mask = _valid_mask(a, b)
    if not mask.any():
        return float("nan")
    return float(np.mean(np.abs(a[mask] - b[mask])))

def max_abs_diff(a: np.ndarray, b: np.ndarray) -> float:
    """Maximum absolute difference over mutually valid elements."""
    mask = _valid_mask(a, b)
    if not mask.any():
        return float("nan")
    return float(np.max(np.abs(a[mask] - b[mask])))

def correlation(a: np.ndarray, b: np.ndarray) -> float:
    """Pearson correlation over mutually valid elements."""
    mask = _valid_mask(a, b)
    if mask.sum() < 2:
        return float("nan")
    r = np.corrcoef(a[mask], b[mask])
    return float(r[0, 1])

def nan_agreement(a: np.ndarray, b: np.ndarray) -> float:
    """Fraction of elements where NaN status agrees (both NaN or both finite)."""
    nan_a = np.isnan(a)
    nan_b = np.isnan(b)
    agree = nan_a == nan_b
    return float(agree.mean())

def get_tolerance(
    variable: str,
    tier: ToleranceTier,
    overrides: dict[str, Tolerance] | None = None,
) -> Tolerance:
    """Look up the tolerance for a variable at a given tier.

    Resolution order:
    1. ``overrides[variable]`` if provided
    2. ``SCIENTIFIC_DEFAULTS[variable]`` if tier is SCIENTIFIC
    3. ``TIER_DEFAULTS[tier]`` (catch-all)
    """
    if overrides and variable in overrides:
        return overrides[variable]
    if tier == ToleranceTier.SCIENTIFIC and variable in SCIENTIFIC_DEFAULTS:
        return SCIENTIFIC_DEFAULTS[variable]
    return TIER_DEFAULTS[tier]

def compare_sbf_vs_rinex(
    store_sbf: Any,
    store_rinex: Any,
    *,
    group: str,
    variables: list[str] | None = None,
) -> ComparisonResult:
    """Compare SBF and RINEX reader outputs from the same observation period.

    Parameters
    ----------
    store_sbf : MyIcechunkStore
        Store populated from SBF files.
    store_rinex : MyIcechunkStore
        Store populated from RINEX files.
    group : str
        Date group key (e.g. "2025001").
    variables : list[str], optional
        Variables to compare. Defaults to intersection.

    Returns
    -------
    ComparisonResult
        SCIENTIFIC tier — expects SNR quantization differences (0.25 dB)
        and potentially different satellite coverage.
    """
    ds_sbf = store_sbf.read_group(group)
    ds_rinex = store_rinex.read_group(group)

    return compare_datasets(
        ds_sbf,
        ds_rinex,
        variables=variables,
        tier=ToleranceTier.SCIENTIFIC,
        label=f"SBF vs RINEX — {group}",
        metadata={
            "comparison_type": "internal",
            "reader_a": "sbf",
            "reader_b": "rinex",
            "group": group,
        },
    )

def compare_ephemeris_sources(
    store_broadcast: Any,
    store_agency: Any,
    *,
    group: str,
    variables: list[str] | None = None,
) -> ComparisonResult:
    """Compare broadcast vs agency (SP3/CLK) ephemeris augmentation.

    Parameters
    ----------
    store_broadcast : MyIcechunkStore
        Store augmented with broadcast ephemeris (from SBF SatVisibility).
    store_agency : MyIcechunkStore
        Store augmented with agency products (SP3 + CLK).
    group : str
        Date group key.

    Returns
    -------
    ComparisonResult
        SCIENTIFIC tier — broadcast ephemeris has lower accuracy than
        precise products (~1-2m vs ~2cm orbit accuracy).
    """
    ds_broadcast = store_broadcast.read_group(group)
    ds_agency = store_agency.read_group(group)

    # Default to satellite coordinate variables
    if variables is None:
        coord_vars = ["sat_x", "sat_y", "sat_z", "phi", "theta"]
        variables = [
            v
            for v in coord_vars
            if v in ds_broadcast.data_vars and v in ds_agency.data_vars
        ]

    return compare_datasets(
        ds_broadcast,
        ds_agency,
        variables=variables,
        tier=ToleranceTier.SCIENTIFIC,
        label=f"Broadcast vs Agency ephemeris — {group}",
        metadata={
            "comparison_type": "internal",
            "ephemeris_a": "broadcast",
            "ephemeris_b": "agency",
            "group": group,
        },
    )

def freeze_checkpoint(
    ds: Any,
    output_path: Path,
    *,
    metadata: dict[str, Any] | None = None,
) -> Path:
    """Save a dataset as a NetCDF checkpoint for future regression testing.

    Parameters
    ----------
    ds : xarray.Dataset
        The known-good reference output.
    output_path : Path
        Where to write the checkpoint file (.nc).
    metadata : dict, optional
        Metadata to store as dataset attributes (git hash, date, config).

    Returns
    -------
    Path
        The written checkpoint path.
    """
    output_path = Path(output_path)
    output_path.parent.mkdir(parents=True, exist_ok=True)

    if metadata:
        ds = ds.assign_attrs(**{f"checkpoint_{k}": str(v) for k, v in metadata.items()})

    ds.to_netcdf(output_path)
    return output_path

def compare_against_checkpoint(
    ds: Any,
    checkpoint_path: Path,
    *,
    variables: list[str] | None = None,
    tier: ToleranceTier = ToleranceTier.EXACT,
    label: str = "",
) -> ComparisonResult:
    """Compare a dataset against a frozen checkpoint.

    Parameters
    ----------
    ds : xarray.Dataset
        Current output to verify.
    checkpoint_path : Path
        Path to the reference checkpoint (.nc file).
    variables : list[str], optional
        Variables to compare. Defaults to intersection.
    tier : ToleranceTier
        Default EXACT — regressions should produce identical output.
    label : str
        Comparison label.

    Returns
    -------
    ComparisonResult
    """
    checkpoint_path = Path(checkpoint_path)
    ds_ref = xr.open_dataset(checkpoint_path)

    return compare_datasets(
        ds,
        ds_ref,
        variables=variables,
        tier=tier,
        label=label or f"Regression check vs {checkpoint_path.name}",
        metadata={
            "comparison_type": "regression",
            "checkpoint_path": str(checkpoint_path),
        },
    )

def compare_vs_gnssvod(
    ds_canvod: Any,
    ds_or_df_gnssvod: Any,
    *,
    band_filter: str = "L1|C",
    snr_col: str = "S1C",
    vod_col: str | None = "VOD1",
    variables: list[str] | None = None,
    label: str = "canvodpy vs gnssvod",
    gnssvod_epoch_col: str = "Epoch",
    gnssvod_sid_col: str = "SV",
) -> ComparisonResult:
    """Compare canvodpy output against gnssvod output for one band.

    Uses ``GnssvodAdapter`` to project canvodpy into gnssvod variable
    space before comparison, ensuring identical variable names, units,
    and conventions.

    Parameters
    ----------
    ds_canvod : xarray.Dataset
        canvodpy output (epoch, sid) dataset.
    ds_or_df_gnssvod : xarray.Dataset or pandas.DataFrame
        gnssvod output. If a pandas DataFrame, it is automatically
        converted to xarray.
    band_filter : str
        SID band suffix, e.g. ``"L1|C"`` or ``"L2|W"``.
    snr_col : str
        gnssvod SNR column name for this band.
    vod_col : str or None
        gnssvod VOD column name for this band.
    variables : list[str], optional
        Variables to compare. Defaults to all shared variables.
    label : str
        Comparison label.
    gnssvod_epoch_col, gnssvod_sid_col : str
        Column names in gnssvod DataFrame.

    Returns
    -------
    ComparisonResult
        SCIENTIFIC tier — two independent implementations.
    """
    import pandas as pd

    if isinstance(ds_or_df_gnssvod, pd.DataFrame):
        ds_gnssvod = gnssvod_df_to_xarray(
            ds_or_df_gnssvod,
            epoch_col=gnssvod_epoch_col,
            sid_col=gnssvod_sid_col,
        )
    else:
        ds_gnssvod = ds_or_df_gnssvod

    # Project canvodpy into gnssvod space
    adapter = GnssvodAdapter(
        ds_canvod,
        band_filter=band_filter,
        snr_col=snr_col,
        vod_col=vod_col,
    )
    ds_adapted = adapter.to_gnssvod_dataset()

    # Handle azimuth wrap-around
    ds_adapted = _wrap_aware_azimuth_diff(ds_adapted, ds_gnssvod)

    # Build tolerance overrides
    tol_overrides = {}
    if snr_col in ds_adapted.data_vars:
        tol_overrides[snr_col] = GNSSVOD_TOLERANCES["SNR"]
    if "Azimuth" in ds_adapted.data_vars:
        tol_overrides["Azimuth"] = GNSSVOD_TOLERANCES["Azimuth"]
    if "Elevation" in ds_adapted.data_vars:
        tol_overrides["Elevation"] = GNSSVOD_TOLERANCES["Elevation"]
    if vod_col and vod_col in ds_adapted.data_vars:
        tol_overrides[vod_col] = GNSSVOD_TOLERANCES["VOD"]

    return compare_datasets(
        ds_adapted,
        ds_gnssvod,
        variables=variables,
        tier=ToleranceTier.SCIENTIFIC,
        tolerance_overrides=tol_overrides,
        label=label,
        metadata={
            "comparison_type": "external",
            "implementation_a": "canvodpy",
            "implementation_b": "gnssvod",
            "band": band_filter,
            "snr_col": snr_col,
            "vod_col": vod_col,
        },
    )

def to_polars(result: ComparisonResult) -> object:
    """Return per-variable stats as a polars DataFrame."""
    return result.to_polars()

def to_markdown(result: ComparisonResult) -> str:
    """Render per-variable stats as a Markdown table."""
    df = result.to_polars()
    if df.is_empty():
        return "_No variables compared._"

    cols = df.columns
    lines = []
    lines.append("| " + " | ".join(cols) + " |")
    lines.append("| " + " | ".join(["---"] * len(cols)) + " |")

    for row in df.iter_rows(named=True):
        cells = []
        for c in cols:
            v = row[c]
            if isinstance(v, float):
                cells.append(f"{v:.6g}")
            else:
                cells.append(str(v))
        lines.append("| " + " | ".join(cells) + " |")

    return "\n".join(lines)

def to_latex(
    result: ComparisonResult,
    *,
    caption: str = "",
    label: str = "",
) -> str:
    """Render per-variable stats as a LaTeX table.

    Produces a ``tabular`` environment inside a ``table`` float, ready
    to paste into a paper.
    """
    df = result.to_polars()
    if df.is_empty():
        return "% No variables compared."

    cols = df.columns
    col_spec = "l" + "r" * (len(cols) - 1)

    lines = [
        r"\begin{table}[htbp]",
        r"\centering",
        f"\\caption{{{caption or result.label}}}",
    ]
    if label:
        lines.append(f"\\label{{{label}}}")

    lines.append(f"\\begin{{tabular}}{{{col_spec}}}")
    lines.append(r"\toprule")

    # Header
    header = " & ".join(f"\\textbf{{{c}}}" for c in cols) + r" \\"
    lines.append(header)
    lines.append(r"\midrule")

    # Rows
    for row in df.iter_rows(named=True):
        cells = []
        for c in cols:
            v = row[c]
            if isinstance(v, float):
                if abs(v) < 1e-3 or abs(v) > 1e4:
                    cells.append(f"${v:.2e}$")
                else:
                    cells.append(f"{v:.4f}")
            else:
                cells.append(str(v).replace("_", r"\_"))
        lines.append(" & ".join(cells) + r" \\")

    lines.append(r"\bottomrule")
    lines.append(r"\end{tabular}")
    lines.append(r"\end{table}")

    return "\n".join(lines)

def plot_diff_histogram(
    ds_a: Any,
    ds_b: Any,
    variable: str,
    *,
    ax: Any = None,
    bins: int = 100,
    title: str | None = None,
) -> Any:
    """Histogram of element-wise differences (a - b) for a single variable.

    Returns a matplotlib Figure.
    """
    import matplotlib.pyplot as plt

    a = ds_a[variable].values.ravel().astype(np.float64)
    b = ds_b[variable].values.ravel().astype(np.float64)
    mask = np.isfinite(a) & np.isfinite(b)
    diff = a[mask] - b[mask]

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

    ax.hist(diff, bins=bins, color="#5D7D5B", edgecolor="#375D3B", alpha=0.8)
    ax.axvline(0, color="#C75050", linewidth=1, linestyle="--")
    ax.set_xlabel(f"Difference ({variable})")
    ax.set_ylabel("Count")
    ax.set_title(title or f"Distribution of differences — {variable}")

    # Annotate with stats
    _text = f"mean={np.mean(diff):.2e}\nstd={np.std(diff):.2e}\nmax|diff|={np.max(np.abs(diff)):.2e}"
    ax.text(
        0.98,
        0.95,
        _text,
        transform=ax.transAxes,
        ha="right",
        va="top",
        fontsize=8,
        family="monospace",
        bbox={"boxstyle": "round", "facecolor": "#E1E6B9", "alpha": 0.8},
    )

    fig.tight_layout()
    return fig

def plot_scatter(
    ds_a: Any,
    ds_b: Any,
    variable: str,
    *,
    ax: Any = None,
    title: str | None = None,
    label_a: str = "Dataset A",
    label_b: str = "Dataset B",
    max_points: int = 50000,
) -> Any:
    """Scatter plot of variable values: a vs b with 1:1 line.

    Returns a matplotlib Figure.
    """
    import matplotlib.pyplot as plt

    a = ds_a[variable].values.ravel().astype(np.float64)
    b = ds_b[variable].values.ravel().astype(np.float64)
    mask = np.isfinite(a) & np.isfinite(b)
    a, b = a[mask], b[mask]

    # Subsample for plotting performance
    if len(a) > max_points:
        idx = np.random.default_rng(42).choice(len(a), max_points, replace=False)
        a, b = a[idx], b[idx]

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

    ax.scatter(b, a, s=1, alpha=0.3, color="#5D7D5B")

    # 1:1 line
    lims = [min(a.min(), b.min()), max(a.max(), b.max())]
    ax.plot(lims, lims, "k--", linewidth=0.8, label="1:1")

    ax.set_xlabel(f"{label_b} — {variable}")
    ax.set_ylabel(f"{label_a} — {variable}")
    ax.set_title(title or f"{variable}: {label_a} vs {label_b}")
    ax.set_aspect("equal")
    ax.legend(loc="upper left")

    fig.tight_layout()
    return fig