diagnostic_validation.py

import json
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import os

# --------------------
# Paths
# --------------------
INPUT_DIR = "results"
OUTPUT_DIR = os.path.join(INPUT_DIR, "diagnostic")
os.makedirs(OUTPUT_DIR, exist_ok=True)

RESULT_FILES = {
    "T1": "results_T1.json",
    "T2": "results_T2.json",
    "T3": "results_T3.json",
}

THRESHOLDS_FILE = "thresholds_t0_nadd3_pctl99.json"
THRESHOLD_OUTPUT_FILE = os.path.join(OUTPUT_DIR, "diagnostic_thresholds.json")

# --------------------
# Load rupture thresholds (per-dimension)
# --------------------
with open(os.path.join(INPUT_DIR, THRESHOLDS_FILE), "r") as f:
    THRESHOLDS = json.load(f)

metrics = ['l2'] + [f'proj{i}' for i in range(10)]
distance_types = ['energy', 'w1']

# --------------------
# Structural metrics
# --------------------
def compute_R(run_results):
    R = []
    for metric in metrics:
        for dist in distance_types:
            observed = run_results[metric][dist]
            tau = THRESHOLDS[metric][f"tau_{dist}"]
            R.append(1 if observed > tau else 0)
    return np.array(R)

def density(R):
    return np.mean(R)

def stability(R_now, R_prev):
    inter = np.sum(np.logical_and(R_now, R_prev))
    union = np.sum(np.logical_or(R_now, R_prev))
    return inter / union if union > 0 else 0

# --------------------
# Load one scenario
# --------------------
def load_case(case_file):
    with open(os.path.join(INPUT_DIR, case_file), "r") as f:
        data = json.load(f)

    runs = []
    for t in range(1, 6):
        R_t = compute_R(data["runs"][str(t)]["results"])
        runs.append({
            "R": R_t,
            "delta": density(R_t)
        })
    return runs

# --------------------
# Build diagnostic summary
# --------------------
def build_summary(runs, max_T=5):
    summary = {"density": [], "persistence": [], "stability": []}

    for T in range(1, max_T + 1):
        deltas = [runs[i]["delta"] for i in range(T)]
        summary["density"].append(deltas[-1])
        summary["persistence"].append(np.mean(deltas))

        if T > 1:
            summary["stability"].append(
                stability(runs[T-1]["R"], runs[T-2]["R"])
            )
        else:
            summary["stability"].append(None)  # estabilidad no definida en T=1

    return summary

# --------------------
# Plotting 2D with thresholds
# --------------------
def plot_2d(metric_data, ylabel, filename, threshold=None, threshold_label=None):
    plt.figure(figsize=(6, 4))

    CASE_COLORS = {
    "T1": "green",   # benigno
    "T2": "red",     # ataque fuerte
    "T3": "orange"   # ataque sigiloso
    }
	
    for case, values in metric_data.items():
        plt.plot(
    		range(1, 6),
    		values,
    		label=case,
    		color=CASE_COLORS.get(case)
		)

    if threshold is not None:
        plt.axhline(
            y=threshold,
            color="gray",
            linestyle="dashed",
            linewidth=1.2,
            label=threshold_label
        )

    plt.xticks(range(1, 6), fontsize=12)
    plt.yticks(fontsize=12)
    plt.xlabel("T", fontsize=14)
    plt.ylabel(ylabel, fontsize=14)
    plt.legend(fontsize=12)
    plt.grid(True)
    plt.tight_layout()
    plt.savefig(filename)
    plt.close()

# --------------------
# Plotting 3D (no thresholds)
# --------------------
def plot_3d(summary_data, filename):
    fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(111, projection='3d')

    colors = {"T1": "green", "T2": "red", "T3": "orange"}

    for case, summary in summary_data.items():
        xs, ys, zs = [], [], []
        for d, p, s in zip(
            summary["density"],
            summary["persistence"],
            summary["stability"]
        ):
            if s is not None:
                xs.append(d)
                ys.append(p)
                zs.append(s)

        ax.plot(xs, ys, zs, marker='o', label=case, color=colors[case])

    ax.set_xlabel("δ(T)", fontsize=12)
    ax.set_ylabel("π(T)", fontsize=12)
    ax.set_zlabel("S(T)", fontsize=12)
    ax.legend()
    plt.tight_layout()
    plt.savefig(filename)
    plt.close()

# --------------------
# Main
# --------------------
if __name__ == "__main__":
    summaries = {}

    # Build summaries
    for case, file in RESULT_FILES.items():
        runs = load_case(file)
        summaries[case] = build_summary(runs)

    # --------------------
    # Save diagnostic summary
    # --------------------
    with open(os.path.join(OUTPUT_DIR, "diagnostic_summary.json"), "w") as f:
        json.dump(summaries, f, indent=2)

    # --------------------
    # Compute diagnostic thresholds from T1 (percentile 95)
    # --------------------
    T1 = summaries["T1"]
    stability_vals = [s for s in T1["stability"] if s is not None]

    diagnostic_thresholds = {
        "density_threshold": float(np.percentile(T1["density"], 95)),
        "persistence_threshold": float(np.percentile(T1["persistence"], 95)),
        "stability_threshold": float(np.percentile(stability_vals, 95))
    }

    with open(THRESHOLD_OUTPUT_FILE, "w") as f_out:
        json.dump(diagnostic_thresholds, f_out, indent=2)

    # --------------------
    # Plot 2D metrics with thresholds
    # --------------------
    labels = {
        "density": r"$\tau_\delta$",
        "persistence": r"$\tau_\pi$",
        "stability": r"$\tau_S$"
    }

    for metric in ["density", "persistence", "stability"]:
        plot_2d(
            metric_data={c: summaries[c][metric] for c in summaries},
            ylabel=f"${metric[0]}(T)$",
            filename=os.path.join(OUTPUT_DIR, f"{metric}_plot.pdf"),
            threshold=diagnostic_thresholds[f"{metric}_threshold"],
            threshold_label=labels[metric]
        )

    # --------------------
    # Plot 3D structural space
    # --------------------
    plot_3d(
        summaries,
        os.path.join(OUTPUT_DIR, "structural_space_3d.pdf")
    )