Visualization of 5D Operating Window with 2D Slices.

This script loads results from the parameter sweep and creates static matplotlib visualizations showing 2D slices through the 5D parameter space with constraint boundaries.

Usage:

python visualize_operating_window.py

Requires:

  • sweep_output/operating_window_results.npz (generated by 2_optim_parameter_sweep_5d.py)

import matplotlib.pyplot as plt
import numpy as np


def load_results(results_path="sweep_output/operating_window_results.npz"):
    """Load results from NPZ file.

    Parameters
    ----------
    results_path : str
        Path to NPZ results file

    Returns
    -------
    results : dict
        Dictionary containing sweep results
    """
    print(f"Loading results from: {results_path}")

    # Load NPZ file
    data = np.load(results_path)

    # Extract results (handle both old and new variable names)
    results = {
        "V": data["V"],
        "e": data["e"],
        "d_t": data["d_t"],
        "f": data["f"],
        "E_p": data["E_p"],
        "breakdown_margin": data["breakdown_margin"],
        "coverage_margin": data["coverage_margin"],
        "density_margin": data["density_margin"],
        "power_error": data["power_error"],
    }

    print(f"Loaded {len(results['V'])} cases")
    return results


def create_2d_slices(results):
    """Create 2D slice visualizations through 5D parameter space.

    Parameters
    ----------
    results : dict
        Results from parameter sweep
    """
    print("\n" + "=" * 70)
    print("  Creating 2D Slice Visualizations")
    print("=" * 70)

    # Define feasibility
    breakdown_ok = results["breakdown_margin"] > 0
    coverage_ok = results["coverage_margin"] > 0
    density_ok = results["density_margin"] > 0
    power_ok = results["power_error"] < 0.1

    feasible = breakdown_ok & coverage_ok & density_ok & power_ok & ~np.isnan(results["V"])

    n_feasible = np.sum(feasible)
    print(f"\nFeasible points: {n_feasible} / {len(results['V'])}")

    if n_feasible == 0:
        print("WARNING: No feasible points found!")

    # Create figure with multiple 2D slices
    fig = plt.figure(figsize=(16, 12))
    fig.suptitle(
        "VedEf Operating Window - 2D Slices Through 5D Space", fontsize=16, fontweight="bold"
    )

    # Slice 1: V vs e (frequency and energy at median values)
    ax1 = fig.add_subplot(2, 3, 1)
    plot_2d_slice(ax1, results, "V", "e", feasible, "V (kV)", "e (mm)", scale_x=1e-3, scale_y=1e3)

    # Slice 2: V vs f
    ax2 = fig.add_subplot(2, 3, 2)
    plot_2d_slice(ax2, results, "V", "f", feasible, "V (kV)", "f (kHz)", scale_x=1e-3, scale_y=1e-3)

    # Slice 3: f vs E_p
    ax3 = fig.add_subplot(2, 3, 3)
    plot_2d_slice(
        ax3, results, "f", "E_p", feasible, "f (kHz)", "E_p (mJ)", scale_x=1e-3, scale_y=1e3
    )

    # Slice 4: e vs d_t
    ax4 = fig.add_subplot(2, 3, 4)
    plot_2d_slice(
        ax4, results, "e", "d_t", feasible, "e (mm)", "d_t (mm)", scale_x=1e3, scale_y=1e3
    )

    # Slice 5: V vs E_p
    ax5 = fig.add_subplot(2, 3, 5)
    plot_2d_slice(
        ax5, results, "V", "E_p", feasible, "V (kV)", "E_p (mJ)", scale_x=1e-3, scale_y=1e3
    )

    # Slice 6: d_t vs f
    ax6 = fig.add_subplot(2, 3, 6)
    plot_2d_slice(
        ax6, results, "d_t", "f", feasible, "d_t (mm)", "f (kHz)", scale_x=1e3, scale_y=1e-3
    )

    plt.tight_layout()

    # Save figure
    output_path = "operating_window_2d_slices.png"
    plt.savefig(output_path, dpi=150, bbox_inches="tight")
    print(f"\n2D slices saved to: {output_path}")

    plt.show()


def plot_2d_slice(
    ax, results, param_x, param_y, feasible, label_x, label_y, scale_x=1.0, scale_y=1.0
):
    """Plot a 2D slice with feasibility coloring.

    Parameters
    ----------
    ax : matplotlib.axes.Axes
        Axes to plot on
    results : dict
        Results dictionary
    param_x : str
        Parameter name for x-axis
    param_y : str
        Parameter name for y-axis
    feasible : ndarray
        Boolean array indicating feasible points
    label_x : str
        Label for x-axis
    label_y : str
        Label for y-axis
    scale_x : float
        Scaling factor for x values
    scale_y : float
        Scaling factor for y values
    """
    x = results[param_x] * scale_x
    y = results[param_y] * scale_y

    # Plot infeasible points
    infeasible = ~feasible & ~np.isnan(x) & ~np.isnan(y)
    if np.sum(infeasible) > 0:
        ax.scatter(
            x[infeasible], y[infeasible], c="red", marker="x", s=30, alpha=0.3, label="Infeasible"
        )

    # Plot feasible points
    if np.sum(feasible) > 0:
        ax.scatter(
            x[feasible],
            y[feasible],
            c="green",
            marker="o",
            s=40,
            alpha=0.7,
            edgecolors="darkgreen",
            linewidth=1,
            label="Feasible",
        )

    ax.set_xlabel(label_x, fontsize=11)
    ax.set_ylabel(label_y, fontsize=11)
    ax.set_title(f"{label_x} vs {label_y}", fontweight="bold")
    ax.grid(True, alpha=0.3)
    ax.legend(loc="best", fontsize=8)


def create_constraint_analysis(results):
    """Create constraint violation analysis plots.

    Parameters
    ----------
    results : dict
        Results from parameter sweep
    """
    print("\nCreating constraint analysis...")

    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))
    fig.suptitle("Constraint Violation Analysis", fontsize=16, fontweight="bold")

    # Filter out NaN values
    valid = ~np.isnan(results["V"])

    # Plot 1: Breakdown margin distribution
    breakdown_margin = results["breakdown_margin"][valid]
    ax1.hist(
        breakdown_margin[breakdown_margin > -50e3],
        bins=50,
        color="blue",
        alpha=0.7,
        edgecolor="black",
    )
    ax1.axvline(x=0, color="r", linestyle="--", linewidth=2, label="Constraint boundary")
    ax1.set_xlabel("Breakdown Margin (V)", fontsize=11)
    ax1.set_ylabel("Count", fontsize=11)
    ax1.set_title("Breakdown Constraint", fontweight="bold")
    ax1.grid(True, alpha=0.3)
    ax1.legend()

    # Plot 2: Coverage margin distribution
    coverage_margin = results["coverage_margin"][valid]
    ax2.hist(coverage_margin, bins=50, color="green", alpha=0.7, edgecolor="black")
    ax2.axvline(x=0, color="r", linestyle="--", linewidth=2, label="Constraint boundary")
    ax2.set_xlabel("Coverage Margin", fontsize=11)
    ax2.set_ylabel("Count", fontsize=11)
    ax2.set_title("Coverage Constraint", fontweight="bold")
    ax2.grid(True, alpha=0.3)
    ax2.legend()

    # Plot 3: Density margin distribution
    density_margin = results["density_margin"][valid]
    ax3.hist(
        density_margin[density_margin > -1e10],
        bins=50,
        color="orange",
        alpha=0.7,
        edgecolor="black",
    )
    ax3.axvline(x=0, color="r", linestyle="--", linewidth=2, label="Constraint boundary")
    ax3.set_xlabel("Energy Density Margin (J/m³)", fontsize=11)
    ax3.set_ylabel("Count", fontsize=11)
    ax3.set_title("Arc Density Constraint", fontweight="bold")
    ax3.grid(True, alpha=0.3)
    ax3.legend()

    # Plot 4: Power error distribution
    power_error = results["power_error"][valid]
    ax4.hist(power_error, bins=50, color="purple", alpha=0.7, edgecolor="black")
    ax4.axvline(x=0.1, color="r", linestyle="--", linewidth=2, label="10% tolerance")
    ax4.axvline(x=-0.1, color="r", linestyle="--", linewidth=2)
    ax4.set_xlabel("Power Error (fraction)", fontsize=11)
    ax4.set_ylabel("Count", fontsize=11)
    ax4.set_title("Power Constraint (25 kW target)", fontweight="bold")
    ax4.grid(True, alpha=0.3)
    ax4.legend()

    plt.tight_layout()

    output_path = "constraint_analysis.png"
    plt.savefig(output_path, dpi=150, bbox_inches="tight")
    print(f"Constraint analysis saved to: {output_path}")

    plt.show()


def main():
    """Load results and create visualizations."""
    try:
        results = load_results()
    except Exception as e:
        print(f"\nError loading results: {e}")
        print("\nPlease run example_operating_window_5d.py first to generate results.")
        return

    # Create visualizations
    create_2d_slices(results)
    create_constraint_analysis(results)

    print("\n" + "=" * 70)
    print("  Visualization complete!")
    print("=" * 70)


if __name__ == "__main__":
    main()