# ============================================================
# Install dependency (Google Colaboratory)
# ============================================================
!pip install scikit-optimize -q

# ============================================================
# Anomaly Detection Parameter Optimization
# Maximizing F1 Score via Bayesian Optimization
# ============================================================

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from sklearn.ensemble import IsolationForest
from sklearn.metrics import f1_score, precision_score, recall_score, confusion_matrix
from sklearn.model_selection import train_test_split
from scipy.interpolate import griddata
from skopt import gp_minimize
from skopt.space import Integer, Real
from skopt.utils import use_named_args
import warnings
warnings.filterwarnings("ignore")

# ------------------------------------------------------------
# 0. Reproducibility
# ------------------------------------------------------------
SEED = 42

# ------------------------------------------------------------
# 1. Synthetic Dataset Generation
# ------------------------------------------------------------
def generate_dataset(n_normal=1000, n_anomaly=88, seed=42):
    rng_local = np.random.default_rng(seed)
    n1 = n_normal // 2
    n2 = n_normal - n1
    X_norm1 = rng_local.multivariate_normal(
        mean=[-2.0, -2.0],
        cov=[[0.6, 0.2], [0.2, 0.6]],
        size=n1
    )
    X_norm2 = rng_local.multivariate_normal(
        mean=[2.0, 2.0],
        cov=[[0.8, -0.3], [-0.3, 0.8]],
        size=n2
    )
    X_normal = np.vstack([X_norm1, X_norm2])
    X_anom   = rng_local.uniform(low=-6.0, high=6.0, size=(n_anomaly, 2))
    X = np.vstack([X_normal, X_anom])
    y = np.array([1] * n_normal + [-1] * n_anomaly)
    return X, y

X, y = generate_dataset()
X_train, X_val, y_train, y_val = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=SEED
)
print(f"Train: {X_train.shape[0]} samples | Val: {X_val.shape[0]} samples")
print(f"Anomaly rate (val): {(y_val == -1).mean():.3f}")

# ------------------------------------------------------------
# 2. Bayesian Optimization
# ------------------------------------------------------------
search_space = [
    Integer(50,   400,  name="n_estimators"),
    Real(   0.1,  1.0,  name="max_samples"),
    Real(   0.01, 0.25, name="contamination"),
]

call_history = []

@use_named_args(search_space)
def objective(n_estimators, max_samples, contamination):
    model = IsolationForest(
        n_estimators=int(n_estimators),
        max_samples=float(max_samples),
        contamination=float(contamination),
        random_state=SEED,
        n_jobs=-1
    )
    model.fit(X_train)
    y_pred = model.predict(X_val)
    f1 = f1_score(y_val, y_pred, pos_label=-1)
    call_history.append({
        "n_estimators":  int(n_estimators),
        "max_samples":   float(max_samples),
        "contamination": float(contamination),
        "f1": f1
    })
    return -f1   # minimize negative F1

print("\nRunning Bayesian Optimization (50 calls)...")
result = gp_minimize(
    func=objective,
    dimensions=search_space,
    n_calls=50,
    n_initial_points=10,
    acq_func="EI",
    random_state=SEED,
    verbose=False
)

best_params = {
    "n_estimators":  int(result.x[0]),
    "max_samples":   float(result.x[1]),
    "contamination": float(result.x[2]),
}
best_f1 = -result.fun
print(f"\nBest Parameters : {best_params}")
print(f"Best F1 Score   : {best_f1:.4f}")

# ------------------------------------------------------------
# 3. Final Model Evaluation
# ------------------------------------------------------------
best_model = IsolationForest(
    n_estimators=best_params["n_estimators"],
    max_samples=best_params["max_samples"],
    contamination=best_params["contamination"],
    random_state=SEED,
    n_jobs=-1
)
best_model.fit(X_train)
y_pred_best = best_model.predict(X_val)

prec    = precision_score(y_val, y_pred_best, pos_label=-1)
rec     = recall_score(y_val,    y_pred_best, pos_label=-1)
f1_opt  = f1_score(y_val,        y_pred_best, pos_label=-1)
cm_best = confusion_matrix(y_val, y_pred_best)

print(f"\nPrecision: {prec:.4f} | Recall: {rec:.4f} | F1: {f1_opt:.4f}")
print(f"Confusion Matrix:\n{cm_best}")

baseline = IsolationForest(contamination=0.08, random_state=SEED, n_jobs=-1)
baseline.fit(X_train)
y_pred_base = baseline.predict(X_val)
f1_base = f1_score(y_val, y_pred_base, pos_label=-1)
print(f"Baseline F1     : {f1_base:.4f}")

# ------------------------------------------------------------
# 4. Visualization
# ------------------------------------------------------------
f1_history  = [-v for v in result.func_vals]
calls_idx   = np.arange(1, len(f1_history) + 1)
best_so_far = np.maximum.accumulate(f1_history)

# Decision boundary meshgrid
x0_min, x0_max = X[:, 0].min() - 0.5, X[:, 0].max() + 0.5
x1_min, x1_max = X[:, 1].min() - 0.5, X[:, 1].max() + 0.5
xx, yy = np.meshgrid(np.linspace(x0_min, x0_max, 300),
                     np.linspace(x1_min, x1_max, 300))
Z = best_model.decision_function(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)

# Build arrays from call history
hist_arr   = np.array([[h["contamination"], h["max_samples"],
                         h["n_estimators"], h["f1"]]
                        for h in call_history])
cont_vals  = hist_arr[:, 0]
ms_vals    = hist_arr[:, 1]
ne_vals    = hist_arr[:, 2]
f1_vals_3d = hist_arr[:, 3]

# 3D surface 1: contamination × max_samples → F1
cont_grid = np.linspace(cont_vals.min(), cont_vals.max(), 60)
ms_grid   = np.linspace(ms_vals.min(),   ms_vals.max(),  60)
CG, MG    = np.meshgrid(cont_grid, ms_grid)
F1G       = griddata((cont_vals, ms_vals), f1_vals_3d,
                     (CG, MG), method="cubic")

# 3D surface 2: n_estimators × contamination → F1
ne_grid  = np.linspace(ne_vals.min(),   ne_vals.max(),   40)
co_grid  = np.linspace(cont_vals.min(), cont_vals.max(), 40)
NEG, COG = np.meshgrid(ne_grid, co_grid)
F1G2     = griddata((ne_vals, cont_vals), f1_vals_3d,
                    (NEG, COG), method="cubic")

# ---- Style constants ----
BG_FIG       = "#0d0d0d"
BG_AX        = "#141414"
TITLE_COLOR  = "#e8e0d0"
LABEL_COLOR  = "#b0a898"
GRID_COLOR   = "#2a2a2a"
ACCENT1      = "#ff6b35"
ACCENT2      = "#00d4ff"
ACCENT_GREEN = "#39ff14"

def style_ax(ax, title="", xlabel="", ylabel=""):
    ax.set_facecolor(BG_AX)
    for sp in ax.spines.values():
        sp.set_edgecolor("#333333")
    ax.tick_params(colors=LABEL_COLOR, labelsize=9)
    ax.xaxis.label.set_color(LABEL_COLOR)
    ax.yaxis.label.set_color(LABEL_COLOR)
    ax.grid(color=GRID_COLOR, linewidth=0.5, linestyle="--")
    if title:
        ax.set_title(title, color=TITLE_COLOR, fontsize=11,
                     fontweight="bold", pad=8)
    if xlabel:
        ax.set_xlabel(xlabel, fontsize=9)
    if ylabel:
        ax.set_ylabel(ylabel, fontsize=9)

fig = plt.figure(figsize=(20, 18), facecolor=BG_FIG)
gs  = gridspec.GridSpec(3, 3, figure=fig,
                        hspace=0.45, wspace=0.35,
                        left=0.06, right=0.97,
                        top=0.93,  bottom=0.06)

# ---- Panel 1: Convergence ----
ax1 = fig.add_subplot(gs[0, 0])
style_ax(ax1, "Bayesian Optimization Convergence",
         "Iteration", "F1 Score (anomaly class)")
ax1.scatter(calls_idx[:10], f1_history[:10],
            color=LABEL_COLOR, s=30, zorder=3,
            label="Random init", alpha=0.7)
ax1.scatter(calls_idx[10:], f1_history[10:],
            color=ACCENT2, s=30, zorder=3,
            label="BO calls", alpha=0.8)
ax1.plot(calls_idx, best_so_far, color=ACCENT1,
         linewidth=2.0, label="Best so far", zorder=4)
ax1.axhline(f1_base, color=ACCENT_GREEN, linewidth=1.2,
            linestyle="--", label=f"Baseline={f1_base:.3f}")
ax1.legend(fontsize=8, facecolor="#1a1a1a",
           labelcolor=TITLE_COLOR, edgecolor="#444444")

# ---- Panel 2: Decision boundary ----
ax2 = fig.add_subplot(gs[0, 1:])
style_ax(ax2, "Isolation Forest Decision Boundary (Optimal Parameters)",
         "Feature 1", "Feature 2")
cf = ax2.contourf(xx, yy, Z, levels=30, cmap="RdYlGn", alpha=0.6)
ax2.contour(xx, yy, Z, levels=[0], colors=ACCENT1, linewidths=1.8)
cb = plt.colorbar(cf, ax=ax2)
cb.ax.tick_params(colors=LABEL_COLOR)
cb.set_label("Anomaly Score", color=LABEL_COLOR, fontsize=9)
mask_n  = y_val == 1
mask_a  = y_val == -1
fp_mask = (y_val == 1)  & (y_pred_best == -1)
fn_mask = (y_val == -1) & (y_pred_best == 1)
ax2.scatter(X_val[mask_n, 0], X_val[mask_n, 1],
            c=ACCENT2, s=12, alpha=0.5, label="Normal", zorder=3)
ax2.scatter(X_val[mask_a, 0], X_val[mask_a, 1],
            c=ACCENT1, s=40, marker="x", linewidths=1.5,
            label="Anomaly (true)", zorder=4)
ax2.scatter(X_val[fp_mask, 0], X_val[fp_mask, 1],
            c="yellow", s=60, marker="^", alpha=0.9,
            label="FP (normal→anomaly)", zorder=5)
ax2.scatter(X_val[fn_mask, 0], X_val[fn_mask, 1],
            c="magenta", s=60, marker="v", alpha=0.9,
            label="FN (anomaly→normal)", zorder=5)
ax2.legend(fontsize=8, facecolor="#1a1a1a",
           labelcolor=TITLE_COLOR, edgecolor="#444444")

# ---- Panel 3: 3D surface contamination × max_samples ----
ax3 = fig.add_subplot(gs[1, :2], projection="3d")
ax3.set_facecolor(BG_AX)
for axis in [ax3.xaxis, ax3.yaxis, ax3.zaxis]:
    axis.pane.fill = False
    axis.pane.set_edgecolor("#333333")
surf3 = ax3.plot_surface(CG, MG, np.nan_to_num(F1G, nan=0.0),
                         cmap="plasma", alpha=0.85, edgecolor="none")
ax3.scatter(best_params["contamination"],
            best_params["max_samples"],
            best_f1,
            color=ACCENT1, s=200, marker="*", zorder=10, label="Optimum")
ax3.set_xlabel("Contamination", color=LABEL_COLOR, fontsize=9, labelpad=8)
ax3.set_ylabel("Max Samples",   color=LABEL_COLOR, fontsize=9, labelpad=8)
ax3.set_zlabel("F1 Score",      color=LABEL_COLOR, fontsize=9, labelpad=8)
ax3.set_title("F1 Surface: contamination × max_samples",
              color=TITLE_COLOR, fontsize=11, fontweight="bold", pad=10)
ax3.tick_params(colors=LABEL_COLOR, labelsize=7)
plt.colorbar(surf3, ax=ax3, shrink=0.4, pad=0.12).ax.tick_params(colors=LABEL_COLOR)
ax3.legend(fontsize=8, facecolor="#1a1a1a",
           labelcolor=TITLE_COLOR, edgecolor="#444444")

# ---- Panel 4: 3D surface contamination × n_estimators ----
ax4 = fig.add_subplot(gs[1, 2], projection="3d")
ax4.set_facecolor(BG_AX)
for axis in [ax4.xaxis, ax4.yaxis, ax4.zaxis]:
    axis.pane.fill = False
    axis.pane.set_edgecolor("#333333")
surf4 = ax4.plot_surface(COG, NEG, np.nan_to_num(F1G2, nan=0.0),
                         cmap="viridis", alpha=0.85, edgecolor="none")
ax4.set_xlabel("Contamination", color=LABEL_COLOR, fontsize=7, labelpad=6)
ax4.set_ylabel("n_estimators",  color=LABEL_COLOR, fontsize=7, labelpad=6)
ax4.set_zlabel("F1",            color=LABEL_COLOR, fontsize=7, labelpad=6)
ax4.set_title("F1 Surface: contamination × n_estimators",
              color=TITLE_COLOR, fontsize=9, fontweight="bold", pad=8)
ax4.tick_params(colors=LABEL_COLOR, labelsize=6)

# ---- Panel 5: Sampled parameter space ----
ax5 = fig.add_subplot(gs[2, 0])
style_ax(ax5, "Sampled Parameters Colored by F1",
         "Contamination", "Max Samples")
sc = ax5.scatter(
    [h["contamination"] for h in call_history],
    [h["max_samples"]   for h in call_history],
    c=[h["f1"]          for h in call_history],
    cmap="plasma", s=50, alpha=0.85, edgecolors="none"
)
ax5.scatter(best_params["contamination"],
            best_params["max_samples"],
            color=ACCENT1, s=200, marker="*", zorder=5, label="Optimum")
cb5 = plt.colorbar(sc, ax=ax5)
cb5.ax.tick_params(colors=LABEL_COLOR)
cb5.set_label("F1", color=LABEL_COLOR, fontsize=8)
ax5.legend(fontsize=8, facecolor="#1a1a1a",
           labelcolor=TITLE_COLOR, edgecolor="#444444")

# ---- Panel 6: F1 per call bar chart ----
ax6 = fig.add_subplot(gs[2, 1])
style_ax(ax6, "F1 Score per BO Call", "Call Index", "F1 Score")
colors_bar = [ACCENT2 if i >= 10 else LABEL_COLOR
              for i in range(len(f1_history))]
ax6.bar(calls_idx, f1_history, color=colors_bar, alpha=0.75, width=0.7)
ax6.axhline(best_f1, color=ACCENT1, linewidth=1.5, linestyle="--",
            label=f"Best={best_f1:.3f}")
ax6.axhline(f1_base, color=ACCENT_GREEN, linewidth=1.2, linestyle="--",
            label=f"Baseline={f1_base:.3f}")
ax6.legend(fontsize=8, facecolor="#1a1a1a",
           labelcolor=TITLE_COLOR, edgecolor="#444444")

# ---- Panel 7: Confusion matrix ----
ax7 = fig.add_subplot(gs[2, 2])
style_ax(ax7, "Confusion Matrix (Optimal Model)")
ax7.set_facecolor(BG_AX)
im = ax7.imshow(cm_best, cmap="YlOrRd", aspect="auto")
for i in range(2):
    for j in range(2):
        ax7.text(j, i, str(cm_best[i, j]),
                 ha="center", va="center",
                 color="white", fontsize=16, fontweight="bold")
ax7.set_xticks([0, 1])
ax7.set_yticks([0, 1])
ax7.set_xticklabels(["Pred Normal", "Pred Anomaly"],
                    color=LABEL_COLOR, fontsize=8)
ax7.set_yticklabels(["True Normal", "True Anomaly"],
                    color=LABEL_COLOR, fontsize=8)
plt.colorbar(im, ax=ax7).ax.tick_params(colors=LABEL_COLOR)

fig.suptitle(
    "Anomaly Detection Parameter Optimization"
    " — Isolation Forest + Bayesian Optimization",
    color=TITLE_COLOR, fontsize=14, fontweight="bold", y=0.97
)
plt.savefig("anomaly_detection_optimization.png", dpi=150,
            bbox_inches="tight", facecolor=fig.get_facecolor())
plt.show()
print("Figure saved.")

# ============================================================
# AML Detection Rule Optimization
# Multi-Objective Combinatorial Optimization
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from mpl_toolkits.mplot3d import Axes3D
from itertools import combinations
import random
import warnings
warnings.filterwarnings('ignore')

# ── Reproducibility ──────────────────────────────────────────
np.random.seed(42)
random.seed(42)

# ============================================================
# 1. PROBLEM DEFINITION
# ============================================================

rules = {
    'names': [
        'Large Cash (>$10K)',
        'Rapid Fund Movement',
        'Structuring Pattern',
        'Shell Company Link',
        'High-Risk Country',
        'Unusual Velocity',
        'Dormant Activation',
        'Round Number Txn',
        'PEP Involvement',
        'Crypto Conversion'
    ],
    'cost':            np.array([3, 5, 4, 6, 2, 3, 4, 2, 7, 8]),
    'detection_rate':  np.array([0.72, 0.85, 0.78, 0.91, 0.65,
                                  0.70, 0.68, 0.55, 0.88, 0.93]),
    'false_pos_rate':  np.array([0.30, 0.20, 0.25, 0.15, 0.40,
                                  0.35, 0.38, 0.50, 0.12, 0.18]),
}

n_rules    = len(rules['names'])
BUDGET     = 20          # max total cost
FP_LIMIT   = 1.50        # max sum of false positive rates
ALPHA      = 0.6         # weight: detection vs false-positive

cost  = rules['cost']
det   = rules['detection_rate']
fp    = rules['false_pos_rate']

def score(subset):
    """Composite score for a given subset (list of indices)."""
    if len(subset) == 0:
        return 0.0
    return ALPHA * det[list(subset)].sum() - (1 - ALPHA) * fp[list(subset)].sum()

def is_feasible(subset):
    """Check budget and false-positive constraints."""
    if len(subset) == 0:
        return False
    idx = list(subset)
    return cost[idx].sum() <= BUDGET and fp[idx].sum() <= FP_LIMIT

# ============================================================
# 2. BRUTE FORCE (Exact Solution)
# ============================================================

def brute_force():
    best_score  = -np.inf
    best_subset = []
    all_results = []

    for r in range(1, n_rules + 1):
        for subset in combinations(range(n_rules), r):
            if is_feasible(subset):
                s = score(subset)
                all_results.append({
                    'subset': subset,
                    'score':  s,
                    'total_det': det[list(subset)].sum(),
                    'total_fp':  fp[list(subset)].sum(),
                    'total_cost': cost[list(subset)].sum(),
                    'n_rules': len(subset)
                })
                if s > best_score:
                    best_score  = s
                    best_subset = subset

    return best_subset, best_score, pd.DataFrame(all_results)

print("Running Brute Force...")
bf_subset, bf_score, bf_results = brute_force()
print(f"  Best subset : {[rules['names'][i] for i in bf_subset]}")
print(f"  Score       : {bf_score:.4f}")
print(f"  Total cost  : {cost[list(bf_subset)].sum()}")
print(f"  Detection Σ : {det[list(bf_subset)].sum():.3f}")
print(f"  False Pos Σ : {fp[list(bf_subset)].sum():.3f}")

# ============================================================
# 3. GREEDY HEURISTIC
# ============================================================

def greedy():
    selected   = []
    remaining  = list(range(n_rules))
    used_cost  = 0
    used_fp    = 0

    while remaining:
        best_idx   = None
        best_ratio = -np.inf
        for i in remaining:
            if used_cost + cost[i] <= BUDGET and used_fp + fp[i] <= FP_LIMIT:
                ratio = (ALPHA * det[i] - (1 - ALPHA) * fp[i]) / cost[i]
                if ratio > best_ratio:
                    best_ratio = ratio
                    best_idx   = i
        if best_idx is None:
            break
        selected.append(best_idx)
        used_cost += cost[best_idx]
        used_fp   += fp[best_idx]
        remaining.remove(best_idx)

    return tuple(selected), score(selected)

print("\nRunning Greedy...")
gr_subset, gr_score = greedy()
print(f"  Best subset : {[rules['names'][i] for i in gr_subset]}")
print(f"  Score       : {gr_score:.4f}")

# ============================================================
# 4. GENETIC ALGORITHM
# ============================================================

POP_SIZE    = 200
N_GEN       = 300
MUT_RATE    = 0.05
ELITE_RATIO = 0.1

def ga_score(chrom):
    subset = [i for i, g in enumerate(chrom) if g == 1]
    if not is_feasible(subset):
        penalty = (max(0, cost[subset].sum() - BUDGET) +
                   max(0, fp[subset].sum()  - FP_LIMIT)) * 2
        return score(subset) - penalty if subset else -999
    return score(subset) if subset else 0

def init_population():
    pop = []
    for _ in range(POP_SIZE):
        chrom = np.zeros(n_rules, dtype=int)
        idx   = random.sample(range(n_rules), random.randint(1, n_rules))
        chrom[idx] = 1
        pop.append(chrom)
    return pop

def tournament(pop, fits, k=5):
    contestants = random.sample(range(len(pop)), k)
    return pop[max(contestants, key=lambda x: fits[x])].copy()

def crossover(p1, p2):
    pt = random.randint(1, n_rules - 1)
    c1 = np.concatenate([p1[:pt], p2[pt:]])
    c2 = np.concatenate([p2[:pt], p1[pt:]])
    return c1, c2

def mutate(chrom):
    for i in range(n_rules):
        if random.random() < MUT_RATE:
            chrom[i] ^= 1
    return chrom

def genetic_algorithm():
    pop  = init_population()
    best_chrom = None
    best_fit   = -np.inf
    history    = []

    n_elite = max(1, int(POP_SIZE * ELITE_RATIO))

    for gen in range(N_GEN):
        fits  = [ga_score(c) for c in pop]
        order = np.argsort(fits)[::-1]

        if fits[order[0]] > best_fit:
            best_fit   = fits[order[0]]
            best_chrom = pop[order[0]].copy()

        history.append(best_fit)
        elites  = [pop[i].copy() for i in order[:n_elite]]
        new_pop = elites[:]

        while len(new_pop) < POP_SIZE:
            p1 = tournament(pop, fits)
            p2 = tournament(pop, fits)
            c1, c2 = crossover(p1, p2)
            new_pop.extend([mutate(c1), mutate(c2)])

        pop = new_pop[:POP_SIZE]

    best_subset = tuple(i for i, g in enumerate(best_chrom) if g == 1)
    return best_subset, score(best_subset), history

print("\nRunning Genetic Algorithm...")
ga_subset, ga_score_val, ga_history = genetic_algorithm()
print(f"  Best subset : {[rules['names'][i] for i in ga_subset]}")
print(f"  Score       : {ga_score_val:.4f}")

# ============================================================
# 5. VISUALIZATION
# ============================================================

fig = plt.figure(figsize=(22, 18))
fig.patch.set_facecolor('#0f0f1a')
plt.rcParams.update({
    'text.color': 'white',
    'axes.labelcolor': 'white',
    'xtick.color': 'white',
    'ytick.color': 'white',
})

CYAN   = '#00e5ff'
GREEN  = '#00ff88'
ORANGE = '#ff6b35'
PURPLE = '#b388ff'
YELLOW = '#ffd740'
PINK   = '#ff4081'

# ── Plot 1: Rule Properties Bar Chart ────────────────────────
ax1 = fig.add_subplot(3, 3, 1)
ax1.set_facecolor('#1a1a2e')
x    = np.arange(n_rules)
w    = 0.3
bars1 = ax1.bar(x - w, det,  w, label='Detection Rate', color=GREEN,  alpha=0.85)
bars2 = ax1.bar(x,      fp,  w, label='False Pos Rate', color=ORANGE, alpha=0.85)
bars3 = ax1.bar(x + w, cost / cost.max(), w, label='Cost (norm)', color=CYAN, alpha=0.85)
ax1.set_xticks(x)
ax1.set_xticklabels([r.split()[0] for r in rules['names']], rotation=45, ha='right', fontsize=7)
ax1.set_title('Rule Properties Overview', color='white', fontweight='bold')
ax1.legend(fontsize=7, facecolor='#1a1a2e', labelcolor='white')
ax1.spines[:].set_color('#333355')
ax1.set_ylim(0, 1.1)

# ── Plot 2: Score Comparison ──────────────────────────────────
ax2 = fig.add_subplot(3, 3, 2)
ax2.set_facecolor('#1a1a2e')
methods = ['Brute Force\n(Exact)', 'Greedy\n(Heuristic)', 'Genetic Alg.\n(Metaheuristic)']
scores  = [bf_score, gr_score, ga_score_val]
colors  = [GREEN, ORANGE, PURPLE]
bars = ax2.bar(methods, scores, color=colors, alpha=0.85, width=0.5)
for bar, s in zip(bars, scores):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005,
             f'{s:.4f}', ha='center', va='bottom', color='white', fontsize=9, fontweight='bold')
ax2.set_title('Score Comparison by Method', color='white', fontweight='bold')
ax2.set_ylabel('Composite Score', color='white')
ax2.spines[:].set_color('#333355')
ax2.set_ylim(0, max(scores) * 1.15)

# ── Plot 3: GA Convergence ────────────────────────────────────
ax3 = fig.add_subplot(3, 3, 3)
ax3.set_facecolor('#1a1a2e')
ax3.plot(ga_history, color=PURPLE, linewidth=1.5, alpha=0.9)
ax3.fill_between(range(len(ga_history)), ga_history, alpha=0.2, color=PURPLE)
ax3.set_title('Genetic Algorithm Convergence', color='white', fontweight='bold')
ax3.set_xlabel('Generation', color='white')
ax3.set_ylabel('Best Score', color='white')
ax3.spines[:].set_color('#333355')
ax3.axhline(bf_score, color=GREEN, linestyle='--', linewidth=1, label=f'Optimal={bf_score:.3f}')
ax3.legend(fontsize=8, facecolor='#1a1a2e', labelcolor='white')

# ── Plot 4: Pareto Frontier (2D) ──────────────────────────────
ax4 = fig.add_subplot(3, 3, 4)
ax4.set_facecolor('#1a1a2e')

feasible = bf_results[bf_results.apply(
    lambda r: cost[list(r['subset'])].sum() <= BUDGET and fp[list(r['subset'])].sum() <= FP_LIMIT,
    axis=1
)].copy()

sc = ax4.scatter(
    feasible['total_fp'], feasible['total_det'],
    c=feasible['score'], cmap='plasma',
    alpha=0.6, s=20, zorder=2
)
plt.colorbar(sc, ax=ax4, label='Score').ax.yaxis.label.set_color('white')

# Mark best solutions
for subset, label, color, marker in [
    (bf_subset, 'Brute Force', GREEN,  '*'),
    (gr_subset, 'Greedy',      ORANGE, 'D'),
    (ga_subset, 'GA',          PURPLE, '^'),
]:
    ax4.scatter(fp[list(subset)].sum(), det[list(subset)].sum(),
                color=color, s=180, marker=marker, zorder=5,
                label=label, edgecolors='white', linewidth=0.8)

ax4.set_xlabel('Total False Positive Rate', color='white')
ax4.set_ylabel('Total Detection Rate', color='white')
ax4.set_title('Pareto Space: Detection vs False Positives', color='white', fontweight='bold')
ax4.legend(fontsize=8, facecolor='#1a1a2e', labelcolor='white')
ax4.spines[:].set_color('#333355')

# ── Plot 5: Rule Selection Heatmap ───────────────────────────
ax5 = fig.add_subplot(3, 3, 5)
ax5.set_facecolor('#1a1a2e')
selection_matrix = np.zeros((3, n_rules))
for j, meth_subset in enumerate([bf_subset, gr_subset, ga_subset]):
    for i in meth_subset:
        selection_matrix[j, i] = 1

im = ax5.imshow(selection_matrix, cmap='YlOrRd', aspect='auto', vmin=0, vmax=1)
ax5.set_yticks([0, 1, 2])
ax5.set_yticklabels(['Brute\nForce', 'Greedy', 'GA'], fontsize=9)
ax5.set_xticks(range(n_rules))
ax5.set_xticklabels([r.split()[0] for r in rules['names']], rotation=45, ha='right', fontsize=7)
ax5.set_title('Rule Selection by Method', color='white', fontweight='bold')
for i in range(3):
    for j in range(n_rules):
        val = int(selection_matrix[i, j])
        ax5.text(j, i, '✓' if val else '', ha='center', va='center',
                 color='white', fontsize=12, fontweight='bold')
ax5.spines[:].set_color('#333355')

# ── Plot 6: Cost-Detection Efficiency Bubble ─────────────────
ax6 = fig.add_subplot(3, 3, 6)
ax6.set_facecolor('#1a1a2e')
efficiency = (ALPHA * det - (1 - ALPHA) * fp) / cost
bubble_colors = plt.cm.cool(efficiency / efficiency.max())
bubbles = ax6.scatter(cost, det, s=fp * 1000, c=efficiency, cmap='cool',
                       alpha=0.8, edgecolors='white', linewidth=0.5)
for i, name in enumerate(rules['names']):
    ax6.annotate(name.split()[0], (cost[i], det[i]),
                 textcoords='offset points', xytext=(5, 3),
                 fontsize=7, color='white', alpha=0.9)
plt.colorbar(bubbles, ax=ax6, label='Efficiency').ax.yaxis.label.set_color('white')
ax6.set_xlabel('Cost', color='white')
ax6.set_ylabel('Detection Rate', color='white')
ax6.set_title('Efficiency Bubble Chart\n(bubble size = false positive rate)', color='white', fontweight='bold')
ax6.spines[:].set_color('#333355')

# ── Plot 7: 3D Pareto Surface ─────────────────────────────────
ax7 = fig.add_subplot(3, 3, 7, projection='3d')
ax7.set_facecolor('#0f0f1a')

sampled = feasible.sample(min(300, len(feasible)), random_state=42)
sc3d = ax7.scatter(
    sampled['total_det'],
    sampled['total_fp'],
    sampled['total_cost'],
    c=sampled['score'], cmap='plasma',
    alpha=0.6, s=20
)
for subset, label, color in [
    (bf_subset, 'Brute Force', GREEN),
    (gr_subset, 'Greedy',      ORANGE),
    (ga_subset, 'GA',          PURPLE),
]:
    ax7.scatter(det[list(subset)].sum(), fp[list(subset)].sum(), cost[list(subset)].sum(),
                color=color, s=200, marker='*', zorder=10, label=label)

ax7.set_xlabel('Detection Rate', color='white', fontsize=8)
ax7.set_ylabel('False Positive', color='white', fontsize=8)
ax7.set_zlabel('Total Cost', color='white', fontsize=8)
ax7.set_title('3D Feasible Solution Space', color='white', fontweight='bold')
ax7.tick_params(colors='white', labelsize=7)
ax7.legend(fontsize=7, facecolor='#0f0f1a', labelcolor='white')
ax7.xaxis.pane.fill = False
ax7.yaxis.pane.fill = False
ax7.zaxis.pane.fill = False

# ── Plot 8: Alpha Sensitivity ─────────────────────────────────
ax8 = fig.add_subplot(3, 3, 8)
ax8.set_facecolor('#1a1a2e')
alphas        = np.linspace(0, 1, 21)
bf_scores_a   = []
n_rules_sel_a = []

for a in alphas:
    best_s = -np.inf
    best_n = 0
    for _, row in feasible.iterrows():
        sub = list(row['subset'])
        s   = a * det[sub].sum() - (1 - a) * fp[sub].sum()
        if s > best_s:
            best_s = s
            best_n = len(sub)
    bf_scores_a.append(best_s)
    n_rules_sel_a.append(best_n)

ax8_twin = ax8.twinx()
ax8.plot(alphas, bf_scores_a,   color=CYAN,   linewidth=2, label='Score')
ax8_twin.plot(alphas, n_rules_sel_a, color=YELLOW, linewidth=2, linestyle='--', label='# Rules')
ax8.axvline(ALPHA, color=PINK, linestyle=':', linewidth=1.5, label=f'Current α={ALPHA}')
ax8.set_xlabel('Alpha (detection weight)', color='white')
ax8.set_ylabel('Best Score', color='white')
ax8_twin.set_ylabel('# Rules Selected', color=YELLOW)
ax8_twin.tick_params(colors='white')
ax8.set_title('Sensitivity to Alpha (α)', color='white', fontweight='bold')
ax8.spines[:].set_color('#333355')
ax8_twin.spines[:].set_color('#333355')
lines1, labels1 = ax8.get_legend_handles_labels()
lines2, labels2 = ax8_twin.get_legend_handles_labels()
ax8.legend(lines1 + lines2, labels1 + labels2, fontsize=7, facecolor='#1a1a2e', labelcolor='white')

# ── Plot 9: Score Distribution ────────────────────────────────
ax9 = fig.add_subplot(3, 3, 9)
ax9.set_facecolor('#1a1a2e')
ax9.hist(feasible['score'], bins=30, color=CYAN, alpha=0.7, edgecolor='white', linewidth=0.3)
ax9.axvline(bf_score, color=GREEN,  linestyle='--', linewidth=2, label=f'Optimal  {bf_score:.3f}')
ax9.axvline(gr_score, color=ORANGE, linestyle='--', linewidth=2, label=f'Greedy   {gr_score:.3f}')
ax9.axvline(ga_score_val, color=PURPLE, linestyle='--', linewidth=2, label=f'GA       {ga_score_val:.3f}')
ax9.set_xlabel('Score', color='white')
ax9.set_ylabel('Frequency', color='white')
ax9.set_title('Score Distribution of Feasible Solutions', color='white', fontweight='bold')
ax9.legend(fontsize=8, facecolor='#1a1a2e', labelcolor='white')
ax9.spines[:].set_color('#333355')

plt.suptitle('AML Detection Rule Optimization Dashboard',
             fontsize=16, color='white', fontweight='bold', y=1.01)
plt.tight_layout()
plt.savefig('aml_optimization.png', dpi=150, bbox_inches='tight',
            facecolor='#0f0f1a')
plt.show()
print("\n[Dashboard saved as aml_optimization.png]")

# ============================================================
# 6. SUMMARY TABLE
# ============================================================

print("\n" + "="*65)
print("  FINAL COMPARISON SUMMARY")
print("="*65)
header = f"{'Method':<20} {'Score':>8} {'Det Σ':>8} {'FP Σ':>8} {'Cost':>6} {'#Rules':>7}"
print(header)
print("-"*65)
for method, subset in [('Brute Force', bf_subset), ('Greedy', gr_subset), ('GA', ga_subset)]:
    idx = list(subset)
    print(f"{method:<20} {score(subset):>8.4f} {det[idx].sum():>8.3f} "
          f"{fp[idx].sum():>8.3f} {cost[idx].sum():>6} {len(idx):>7}")
print("="*65)
print(f"\nBudget limit : {BUDGET}   |   FP limit : {FP_LIMIT}   |   Alpha : {ALPHA}")

# ============================================================
# Credit Card Fraud Detection: Feature Selection Optimization
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from mpl_toolkits.mplot3d import Axes3D

from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import (
    SelectKBest, mutual_info_classif, RFE, SelectFromModel
)
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.metrics import average_precision_score

import warnings
warnings.filterwarnings("ignore")

np.random.seed(42)

# ============================================================
# 1. Synthetic Dataset Generation
# ============================================================

n_samples = 5000
fraud_rate = 0.05

n_fraud = int(n_samples * fraud_rate)
n_legit = n_samples - n_fraud

def generate_transactions(n, is_fraud):
    rng = np.random.RandomState(42 if not is_fraud else 99)
    data = {
        "amount":             rng.exponential(scale=500 if is_fraud else 100, size=n),
        "hour":               rng.choice(np.arange(0, 6), size=n) if is_fraud
                              else rng.randint(8, 22, size=n),
        "distance_from_home": rng.exponential(scale=80 if is_fraud else 10, size=n),
        "velocity_1h":        rng.poisson(lam=5 if is_fraud else 1, size=n),
        "velocity_24h":       rng.poisson(lam=15 if is_fraud else 4, size=n),
        "merchant_risk":      rng.beta(a=8, b=2, size=n) if is_fraud
                              else rng.beta(a=2, b=8, size=n),
        "foreign_transaction":rng.binomial(1, 0.7 if is_fraud else 0.05, size=n),
        "chip_used":          rng.binomial(1, 0.2 if is_fraud else 0.9, size=n),
        "online_order":       rng.binomial(1, 0.8 if is_fraud else 0.3, size=n),
        "noise_1":            rng.randn(n),
        "noise_2":            rng.randn(n),
        "noise_3":            rng.randn(n),
        "noise_4":            rng.randn(n),
        "noise_5":            rng.randn(n),
    }
    return pd.DataFrame(data)

df_legit  = generate_transactions(n_legit,  is_fraud=False)
df_fraud  = generate_transactions(n_fraud,  is_fraud=True)
df_legit["fraud"]  = 0
df_fraud["fraud"]  = 1

df = pd.concat([df_legit, df_fraud], ignore_index=True).sample(
    frac=1, random_state=42
).reset_index(drop=True)

X = df.drop("fraud", axis=1)
y = df["fraud"]

feature_names = list(X.columns)
print(f"Dataset shape: {X.shape}, Fraud rate: {y.mean():.2%}")

# ============================================================
# 2. Feature Selection Methods
# ============================================================

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
base_clf = LogisticRegression(max_iter=1000, random_state=42)

# ---------- 2-A. Filter: Mutual Information ----------
mi_scores = mutual_info_classif(X_scaled, y, random_state=42)
mi_ranking = np.argsort(mi_scores)[::-1]

# ---------- 2-B. Wrapper: RFE ----------
rfe = RFE(estimator=LogisticRegression(max_iter=1000, random_state=42),
          n_features_to_select=7, step=1)
rfe.fit(X_scaled, y)
rfe_support = rfe.support_
rfe_ranking = rfe.ranking_

# ---------- 2-C. Embedded: LASSO (L1) ----------
lasso = LogisticRegression(
    penalty="l1", solver="liblinear", C=0.1,
    max_iter=1000, random_state=42
)
lasso.fit(X_scaled, y)
lasso_coef = np.abs(lasso.coef_[0])

# ---------- 2-D. Tree-based: Random Forest ----------
rf = RandomForestClassifier(
    n_estimators=200, n_jobs=-1, random_state=42
)
rf.fit(X_scaled, y)
rf_importance = rf.feature_importances_

# ============================================================
# 3. Evaluate AP for Different Feature Subsets (Filter Method)
# ============================================================

k_values = list(range(1, len(feature_names) + 1))
ap_scores_k = []

for k in k_values:
    top_k_idx = mi_ranking[:k]
    X_k = X_scaled[:, top_k_idx]
    scores = cross_val_score(
        base_clf, X_k, y,
        cv=cv, scoring="average_precision", n_jobs=-1
    )
    ap_scores_k.append(scores.mean())

best_k = k_values[np.argmax(ap_scores_k)]
best_ap = max(ap_scores_k)
print(f"Best k (Filter/MI): {best_k}, AP: {best_ap:.4f}")

# ============================================================
# 4. Compare All Methods
# ============================================================

def evaluate_subset(indices):
    X_sub = X_scaled[:, indices]
    scores = cross_val_score(
        base_clf, X_sub, y,
        cv=cv, scoring="average_precision", n_jobs=-1
    )
    return scores.mean(), scores.std()

# Filter top-k
filter_idx  = mi_ranking[:best_k]
filter_ap, filter_std = evaluate_subset(filter_idx)

# RFE selected
rfe_idx     = np.where(rfe_support)[0]
rfe_ap, rfe_std = evaluate_subset(rfe_idx)

# LASSO nonzero
lasso_idx   = np.where(lasso_coef > 0)[0]
lasso_ap, lasso_std = evaluate_subset(lasso_idx)

# RF top-k (same k as filter)
rf_idx      = np.argsort(rf_importance)[::-1][:best_k]
rf_ap, rf_std = evaluate_subset(rf_idx)

# All features (baseline)
all_idx     = np.arange(len(feature_names))
all_ap, all_std = evaluate_subset(all_idx)

method_results = {
    "All Features\n(Baseline)": (all_ap,    all_std,    len(all_idx)),
    "Filter\n(MI SelectKBest)": (filter_ap, filter_std, len(filter_idx)),
    "Wrapper\n(RFE)":           (rfe_ap,    rfe_std,    len(rfe_idx)),
    "Embedded\n(LASSO L1)":     (lasso_ap,  lasso_std,  len(lasso_idx)),
    "Tree-based\n(RF Importance)":(rf_ap,   rf_std,     len(rf_idx)),
}
print("\n--- Method Comparison ---")
for m, (ap, std, nf) in method_results.items():
    print(f"{m.replace(chr(10),' '):<35} AP={ap:.4f} ± {std:.4f}  #features={nf}")

# ============================================================
# 5. Visualization
# ============================================================

plt.rcParams.update({
    "figure.facecolor": "white",
    "axes.facecolor":   "#f8f9fa",
    "axes.grid":        True,
    "grid.alpha":       0.4,
    "font.size":        11,
})

fig = plt.figure(figsize=(22, 20))
gs  = gridspec.GridSpec(3, 3, figure=fig, hspace=0.45, wspace=0.38)

colors_feat = plt.cm.RdYlGn(np.linspace(0.15, 0.85, len(feature_names)))

# ---- Plot 1: MI scores (bar) ----
ax1 = fig.add_subplot(gs[0, 0])
sorted_mi  = mi_scores[mi_ranking]
sorted_names = [feature_names[i] for i in mi_ranking]
bars = ax1.barh(sorted_names[::-1], sorted_mi[::-1],
                color=colors_feat, edgecolor="white")
ax1.set_title("Filter Method\nMutual Information Scores", fontweight="bold")
ax1.set_xlabel("MI Score")
for bar, val in zip(bars, sorted_mi[::-1]):
    ax1.text(val + 0.001, bar.get_y() + bar.get_height()/2,
             f"{val:.3f}", va="center", fontsize=8)

# ---- Plot 2: RFE Ranking ----
ax2 = fig.add_subplot(gs[0, 1])
rfe_rank_sorted_idx = np.argsort(rfe_ranking)
rfe_names = [feature_names[i] for i in rfe_rank_sorted_idx]
rfe_ranks = rfe_ranking[rfe_rank_sorted_idx]
bar_colors = ["#2ecc71" if s else "#e74c3c" for s in rfe_support[rfe_rank_sorted_idx]]
ax2.barh(rfe_names, rfe_ranks, color=bar_colors, edgecolor="white")
ax2.set_title("Wrapper Method\nRFE Ranking (green = selected)", fontweight="bold")
ax2.set_xlabel("RFE Rank (lower = better)")
from matplotlib.patches import Patch
ax2.legend(handles=[Patch(color="#2ecc71", label="Selected"),
                    Patch(color="#e74c3c", label="Eliminated")], fontsize=9)

# ---- Plot 3: LASSO Coefficients ----
ax3 = fig.add_subplot(gs[0, 2])
lasso_sorted_idx = np.argsort(lasso_coef)[::-1]
ax3.bar([feature_names[i] for i in lasso_sorted_idx],
        lasso_coef[lasso_sorted_idx],
        color=["#3498db" if lasso_coef[i] > 0 else "#bdc3c7"
               for i in lasso_sorted_idx],
        edgecolor="white")
ax3.set_title("Embedded Method\nLASSO |Coefficients|", fontweight="bold")
ax3.set_ylabel("|Coefficient|")
ax3.set_xticklabels([feature_names[i] for i in lasso_sorted_idx],
                    rotation=45, ha="right", fontsize=8)

# ---- Plot 4: RF Importance ----
ax4 = fig.add_subplot(gs[1, 0])
rf_sorted_idx = np.argsort(rf_importance)[::-1]
cmap_rf = plt.cm.Blues(np.linspace(0.4, 0.9, len(feature_names)))
ax4.bar([feature_names[i] for i in rf_sorted_idx],
        rf_importance[rf_sorted_idx],
        color=cmap_rf[::-1], edgecolor="white")
ax4.set_title("Tree-based Method\nRandom Forest Importance", fontweight="bold")
ax4.set_ylabel("Importance")
ax4.set_xticklabels([feature_names[i] for i in rf_sorted_idx],
                    rotation=45, ha="right", fontsize=8)

# ---- Plot 5: AP vs k (filter method) ----
ax5 = fig.add_subplot(gs[1, 1])
ax5.plot(k_values, ap_scores_k, "o-", color="#9b59b6", lw=2, ms=6)
ax5.axvline(best_k, color="red", ls="--", lw=1.5, label=f"Best k={best_k}")
ax5.axhline(all_ap, color="gray", ls=":", lw=1.5, label=f"Baseline (all) AP={all_ap:.3f}")
ax5.scatter([best_k], [best_ap], color="red", s=120, zorder=5)
ax5.set_title("Filter Method\nAP vs Number of Features (k)", fontweight="bold")
ax5.set_xlabel("k (number of features)")
ax5.set_ylabel("Average Precision (CV)")
ax5.legend(fontsize=9)

# ---- Plot 6: Method Comparison ----
ax6 = fig.add_subplot(gs[1, 2])
method_names = list(method_results.keys())
ap_vals  = [v[0] for v in method_results.values()]
std_vals = [v[1] for v in method_results.values()]
n_feats  = [v[2] for v in method_results.values()]
bar_col  = ["#95a5a6", "#e67e22", "#e74c3c", "#3498db", "#2ecc71"]
bars6 = ax6.bar(method_names, ap_vals, yerr=std_vals,
                color=bar_col, edgecolor="white",
                capsize=5, alpha=0.9)
ax6.set_title("Method Comparison\nAverage Precision (5-Fold CV)", fontweight="bold")
ax6.set_ylabel("Average Precision")
ax6.set_ylim(min(ap_vals) - 0.05, max(ap_vals) + 0.05)
for bar, ap, nf in zip(bars6, ap_vals, n_feats):
    ax6.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005,
             f"AP={ap:.3f}\nn={nf}", ha="center", va="bottom", fontsize=8)
ax6.tick_params(axis="x", labelsize=8)

# ---- Plot 7: 3D — amount × distance × velocity_1h (fraud vs legit) ----
ax7 = fig.add_subplot(gs[2, 0], projection="3d")
idx_legit_s = df[df.fraud == 0].sample(300, random_state=42).index
idx_fraud_s = df[df.fraud == 1].sample(min(200, n_fraud), random_state=42).index
ax7.scatter(df.loc[idx_legit_s, "amount"],
            df.loc[idx_legit_s, "distance_from_home"],
            df.loc[idx_legit_s, "velocity_1h"],
            c="#3498db", alpha=0.4, s=15, label="Legit")
ax7.scatter(df.loc[idx_fraud_s, "amount"],
            df.loc[idx_fraud_s, "distance_from_home"],
            df.loc[idx_fraud_s, "velocity_1h"],
            c="#e74c3c", alpha=0.7, s=25, label="Fraud")
ax7.set_xlabel("Amount", fontsize=8)
ax7.set_ylabel("Distance\nfrom Home", fontsize=8)
ax7.set_zlabel("Velocity\n(1h)", fontsize=8)
ax7.set_title("3D Feature Space\nAmount × Distance × Velocity", fontweight="bold")
ax7.legend(fontsize=8)

# ---- Plot 8: 3D — MI score × RF importance × LASSO coef ----
ax8 = fig.add_subplot(gs[2, 1], projection="3d")
mi_n  = mi_scores / (mi_scores.max() + 1e-9)
rf_n  = rf_importance / (rf_importance.max() + 1e-9)
ls_n  = lasso_coef / (lasso_coef.max() + 1e-9)
noise_mask = np.array(["noise" in f for f in feature_names])
sc = ax8.scatter(mi_n[~noise_mask], rf_n[~noise_mask], ls_n[~noise_mask],
                 c="#2ecc71", s=80, label="Signal", zorder=5)
ax8.scatter(mi_n[noise_mask],  rf_n[noise_mask],  ls_n[noise_mask],
            c="#e74c3c", s=60, marker="x", label="Noise", zorder=5)
for i, name in enumerate(feature_names):
    ax8.text(mi_n[i], rf_n[i], ls_n[i], name, fontsize=6)
ax8.set_xlabel("MI Score\n(normalized)", fontsize=8)
ax8.set_ylabel("RF Importance\n(normalized)", fontsize=8)
ax8.set_zlabel("LASSO |Coef|\n(normalized)", fontsize=8)
ax8.set_title("3D Method Agreement\nMI × RF × LASSO", fontweight="bold")
ax8.legend(fontsize=8)

# ---- Plot 9: Feature score heatmap ----
ax9 = fig.add_subplot(gs[2, 2])
score_matrix = np.column_stack([mi_n, rf_n, ls_n])
im = ax9.imshow(score_matrix, aspect="auto", cmap="YlOrRd")
ax9.set_xticks([0, 1, 2])
ax9.set_xticklabels(["MI", "RF", "LASSO"], fontsize=10)
ax9.set_yticks(range(len(feature_names)))
ax9.set_yticklabels(feature_names, fontsize=9)
ax9.set_title("Feature Score Heatmap\n(normalized per method)", fontweight="bold")
plt.colorbar(im, ax=ax9, fraction=0.046, pad=0.04)

fig.suptitle(
    "Credit Card Fraud Detection — Feature Selection Optimization",
    fontsize=16, fontweight="bold", y=1.01
)
plt.savefig("fraud_feature_selection.png", dpi=150,
            bbox_inches="tight", facecolor="white")
plt.show()
print("Figure saved.")

# ============================================================
# Real-Time Monitoring Alert Threshold Optimization
# Minimizing Operational Load
# ============================================================

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib import cm
from scipy import stats
from scipy.optimize import minimize_scalar, minimize
from scipy.stats import norm
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')

# ── Seaborn-style without importing seaborn ──────────────────
plt.rcParams.update({
    'axes.facecolor':    '#f8f9fa',
    'figure.facecolor':  'white',
    'axes.grid':         True,
    'grid.color':        'white',
    'grid.linewidth':    1.2,
    'axes.spines.top':   False,
    'axes.spines.right': False,
    'font.size':         11,
})

# ============================================================
# 1. System Parameters
# ============================================================
# Normal operation distribution
MU_NORMAL    = 45.0   # mean CPU % under normal load
SIGMA_NORMAL = 10.0   # std  CPU % under normal load

# Incident distribution
MU_INC    = 80.0   # mean CPU % during incident
SIGMA_INC =  8.0   # std  CPU % during incident

# Cost weights
LAMBDA_FP  =  1.0   # false-positive cost (wasted operator time)
LAMBDA_FN  = 10.0   # false-negative cost (missed incident / business impact)
LAMBDA_VOL =  0.5   # alert-volume cost   (noise fatigue)

# Operational context
INCIDENTS_PER_DAY  = 5
CHECKS_PER_HOUR    = 60
HOURS_PER_DAY      = 24
CHECKS_PER_DAY     = CHECKS_PER_HOUR * HOURS_PER_DAY

# ============================================================
# 2. Core Functions (vectorised for speed)
# ============================================================
def fpr(theta, mu=MU_NORMAL, sigma=SIGMA_NORMAL):
    """False Positive Rate: P(X > theta | normal)"""
    return 1.0 - norm.cdf(theta, loc=mu, scale=sigma)

def fnr(theta, mu=MU_INC, sigma=SIGMA_INC):
    """False Negative Rate: P(X <= theta | incident)"""
    return norm.cdf(theta, loc=mu, scale=sigma)

def alert_volume_per_hour(theta):
    """Expected number of alerts fired per hour."""
    # Normal alerts (false positives)
    normal_checks = CHECKS_PER_HOUR
    vol_normal = normal_checks * fpr(theta)

    # Incident alerts (true positives) — incident lasts ~15 min on average
    incident_checks_per_hour = (INCIDENTS_PER_DAY / HOURS_PER_DAY) * 15
    vol_incident = incident_checks_per_hour * (1.0 - fnr(theta))

    return vol_normal + vol_incident

def operational_cost(theta):
    """
    Total cost C(theta) = lambda_fp * FPR + lambda_fn * FNR + lambda_vol * V_bar
    Vectorised: theta can be a scalar or ndarray.
    """
    theta = np.asarray(theta, dtype=float)
    fp = fpr(theta)
    fn = fnr(theta)
    vol = alert_volume_per_hour(theta)
    # Normalise volume to [0,1] range for comparability
    vol_norm = vol / CHECKS_PER_HOUR
    return LAMBDA_FP * fp + LAMBDA_FN * fn + LAMBDA_VOL * vol_norm

# ============================================================
# 3. Optimisation  (scipy.optimize — fast bounded search)
# ============================================================
result = minimize_scalar(
    operational_cost,
    bounds=(MU_NORMAL, MU_INC),
    method='bounded',
    options={'xatol': 1e-6}
)
THETA_OPT = result.x
COST_OPT  = result.fun

# Sensitivity: how cost changes ±5 pp around optimum
DELTA = 5.0
cost_low  = operational_cost(THETA_OPT - DELTA)
cost_high = operational_cost(THETA_OPT + DELTA)

print("=" * 55)
print("  ALERT THRESHOLD OPTIMISATION RESULTS")
print("=" * 55)
print(f"  Optimal threshold  θ*  : {THETA_OPT:.2f}%")
print(f"  Minimum cost       C*  : {COST_OPT:.4f}")
print(f"  FPR at θ*              : {fpr(THETA_OPT)*100:.2f}%")
print(f"  FNR at θ*              : {fnr(THETA_OPT)*100:.2f}%")
print(f"  Alert volume/hr    V̄   : {alert_volume_per_hour(THETA_OPT):.2f} alerts")
print(f"  Cost at θ*-5pp         : {cost_low:.4f}  (Δ={cost_low-COST_OPT:+.4f})")
print(f"  Cost at θ*+5pp         : {cost_high:.4f}  (Δ={cost_high-COST_OPT:+.4f})")
print("=" * 55)

# ============================================================
# 4. Multi-weight Sensitivity Grid (vectorised, fast)
# ============================================================
lam_fp_vals = np.linspace(0.5, 5.0, 40)
lam_fn_vals = np.linspace(1.0, 20.0, 40)
LFP, LFN   = np.meshgrid(lam_fp_vals, lam_fn_vals)

# For each (lam_fp, lam_fn) combo, find optimal theta via a dense grid search
# (faster than calling minimize_scalar 1600 times)
theta_grid = np.linspace(MU_NORMAL, MU_INC, 500)
fp_grid    = fpr(theta_grid)
fn_grid    = fnr(theta_grid)
vol_grid   = alert_volume_per_hour(theta_grid) / CHECKS_PER_HOUR

# Shape: (len_lam_fp, len_lam_fn, len_theta)
cost_3d = (
    LFP[:, :, np.newaxis] * fp_grid[np.newaxis, np.newaxis, :]
    + LFN[:, :, np.newaxis] * fn_grid[np.newaxis, np.newaxis, :]
    + LAMBDA_VOL * vol_grid[np.newaxis, np.newaxis, :]
)
opt_idx    = np.argmin(cost_3d, axis=2)
THETA_SURF = theta_grid[opt_idx]          # (40, 40) surface

# ============================================================
# 5. Simulate a 24-hour time series
# ============================================================
rng = np.random.default_rng(42)
t   = np.arange(0, CHECKS_PER_DAY)

# Base CPU: sinusoidal daily cycle + noise
cpu_base = (
    MU_NORMAL
    + 8.0 * np.sin(2 * np.pi * t / CHECKS_PER_DAY)
    + rng.normal(0, SIGMA_NORMAL * 0.7, size=len(t))
)

# Inject 5 incidents
incident_times = rng.choice(t[200:-200], size=INCIDENTS_PER_DAY, replace=False)
cpu_sim = cpu_base.copy()
for it in incident_times:
    dur = rng.integers(10, 40)
    end = min(it + dur, CHECKS_PER_DAY)
    cpu_sim[it:end] += rng.normal(MU_INC - MU_NORMAL, SIGMA_INC * 0.5, size=end - it)

# Alert flags for multiple thresholds
THETA_NAIVE = 70.0  # typical naive threshold
alerts_opt   = cpu_sim > THETA_OPT
alerts_naive = cpu_sim > THETA_NAIVE

# ============================================================
# 6. Plots
# ============================================================
fig = plt.figure(figsize=(22, 26))
fig.suptitle(
    "Real-Time Monitoring Alert Threshold Optimisation\n"
    "Minimising Operational Load",
    fontsize=16, fontweight='bold', y=0.98
)

COLOR_OPT   = '#2ecc71'
COLOR_NAIVE = '#e74c3c'
COLOR_FPR   = '#3498db'
COLOR_FNR   = '#e67e22'
COLOR_COST  = '#9b59b6'

theta_range = np.linspace(20, 100, 800)

# ── Plot 1: CPU distributions ─────────────────────────────
ax1 = fig.add_subplot(4, 2, 1)
x   = np.linspace(0, 120, 1000)
ax1.fill_between(x, norm.pdf(x, MU_NORMAL, SIGMA_NORMAL),
                 alpha=0.35, color=COLOR_FPR, label='Normal operation')
ax1.fill_between(x, norm.pdf(x, MU_INC, SIGMA_INC),
                 alpha=0.35, color=COLOR_NAIVE, label='Incident')
ax1.plot(x, norm.pdf(x, MU_NORMAL, SIGMA_NORMAL), color=COLOR_FPR, lw=2)
ax1.plot(x, norm.pdf(x, MU_INC,    SIGMA_INC),    color=COLOR_NAIVE, lw=2)
ax1.axvline(THETA_OPT,   color=COLOR_OPT,   lw=2.5, ls='--', label=f'θ* = {THETA_OPT:.1f}%')
ax1.axvline(THETA_NAIVE, color=COLOR_NAIVE, lw=2,   ls=':',  label=f'θ_naive = {THETA_NAIVE}%')
ax1.set_xlabel('CPU Usage (%)'); ax1.set_ylabel('Density')
ax1.set_title('CPU Distribution: Normal vs Incident')
ax1.legend(fontsize=9)

# ── Plot 2: FPR & FNR curves ──────────────────────────────
ax2 = fig.add_subplot(4, 2, 2)
ax2.plot(theta_range, fpr(theta_range) * 100, color=COLOR_FPR, lw=2.5, label='FPR (%)')
ax2.plot(theta_range, fnr(theta_range) * 100, color=COLOR_FNR, lw=2.5, label='FNR (%)')
ax2.axvline(THETA_OPT,   color=COLOR_OPT,   lw=2.5, ls='--', label=f'θ* = {THETA_OPT:.1f}%')
ax2.axvline(THETA_NAIVE, color=COLOR_NAIVE, lw=2,   ls=':',  label=f'θ_naive = {THETA_NAIVE}%')
ax2.set_xlabel('Threshold θ (%)'); ax2.set_ylabel('Rate (%)')
ax2.set_title('False Positive / Negative Rates vs Threshold')
ax2.legend(fontsize=9)

# ── Plot 3: Cost function ─────────────────────────────────
ax3 = fig.add_subplot(4, 2, 3)
cost_vals = operational_cost(theta_range)
ax3.plot(theta_range, cost_vals, color=COLOR_COST, lw=2.5)
ax3.axvline(THETA_OPT,   color=COLOR_OPT,   lw=2.5, ls='--', label=f'θ* = {THETA_OPT:.1f}%')
ax3.axvline(THETA_NAIVE, color=COLOR_NAIVE, lw=2,   ls=':',  label=f'θ_naive = {THETA_NAIVE}%')
ax3.scatter([THETA_OPT], [COST_OPT], color=COLOR_OPT, s=120, zorder=5)
ax3.set_xlabel('Threshold θ (%)'); ax3.set_ylabel('Cost C(θ)')
ax3.set_title('Operational Cost vs Threshold')
ax3.legend(fontsize=9)

# ── Plot 4: Alert volume ───────────────────────────────────
ax4 = fig.add_subplot(4, 2, 4)
vol_vals = alert_volume_per_hour(theta_range)
ax4.plot(theta_range, vol_vals, color='#1abc9c', lw=2.5)
ax4.axvline(THETA_OPT,   color=COLOR_OPT,   lw=2.5, ls='--', label=f'θ* = {THETA_OPT:.1f}%')
ax4.axvline(THETA_NAIVE, color=COLOR_NAIVE, lw=2,   ls=':',  label=f'θ_naive = {THETA_NAIVE}%')
ax4.set_xlabel('Threshold θ (%)'); ax4.set_ylabel('Alerts / Hour')
ax4.set_title('Expected Alert Volume per Hour')
ax4.legend(fontsize=9)

# ── Plot 5 & 6: 24-hour simulation ────────────────────────
hours = t / CHECKS_PER_HOUR
ax5 = fig.add_subplot(4, 2, (5, 6))
ax5.plot(hours, cpu_sim, color='#95a5a6', lw=0.8, alpha=0.7, label='CPU usage')
ax5.axhline(THETA_OPT,   color=COLOR_OPT,   lw=2, ls='--', label=f'θ* = {THETA_OPT:.1f}%')
ax5.axhline(THETA_NAIVE, color=COLOR_NAIVE, lw=2, ls=':',  label=f'θ_naive = {THETA_NAIVE}%')
# Shade alerts
ax5.fill_between(hours, cpu_sim, THETA_OPT,
                 where=alerts_opt, alpha=0.4, color=COLOR_OPT,   label='Alert (opt)')
ax5.fill_between(hours, cpu_sim, THETA_NAIVE,
                 where=alerts_naive, alpha=0.25, color=COLOR_NAIVE, label='Alert (naive)')
# Mark incidents
for it in incident_times:
    ax5.axvline(it / CHECKS_PER_HOUR, color='black', lw=1, ls=':', alpha=0.5)
ax5.set_xlabel('Time (hours)'); ax5.set_ylabel('CPU Usage (%)')
ax5.set_title('Simulated 24-Hour CPU Timeline with Alerts')
ax5.legend(fontsize=9, ncol=3)

# ── Plot 7: ROC-style trade-off ───────────────────────────
ax6 = fig.add_subplot(4, 2, 7)
fpr_vals = fpr(theta_range) * 100
fnr_vals = fnr(theta_range) * 100
sc = ax6.scatter(fpr_vals, fnr_vals,
                 c=theta_range, cmap='plasma', s=8, zorder=3)
plt.colorbar(sc, ax=ax6, label='Threshold θ (%)')
ax6.scatter([fpr(THETA_OPT) * 100], [fnr(THETA_OPT) * 100],
            color=COLOR_OPT, s=200, zorder=5, marker='*', label='θ*')
ax6.scatter([fpr(THETA_NAIVE) * 100], [fnr(THETA_NAIVE) * 100],
            color=COLOR_NAIVE, s=150, zorder=5, marker='D', label='θ_naive')
ax6.set_xlabel('FPR (%)'); ax6.set_ylabel('FNR (%)')
ax6.set_title('FPR vs FNR Trade-Off Curve')
ax6.legend(fontsize=9)

# ── Plot 8: Cost decomposition bar ───────────────────────
ax7 = fig.add_subplot(4, 2, 8)
thetas_bar  = [THETA_OPT, THETA_NAIVE, 60.0, 55.0]
labels_bar  = [f'θ*={THETA_OPT:.1f}%', 'Naive=70%', 'θ=60%', 'θ=55%']
c_fp  = [LAMBDA_FP  * fpr(th) for th in thetas_bar]
c_fn  = [LAMBDA_FN  * fnr(th) for th in thetas_bar]
c_vol = [LAMBDA_VOL * alert_volume_per_hour(th) / CHECKS_PER_HOUR for th in thetas_bar]
x_pos = np.arange(len(thetas_bar))
w = 0.55
ax7.bar(x_pos, c_fp,  w, label='FP cost',     color=COLOR_FPR,  alpha=0.85)
ax7.bar(x_pos, c_fn,  w, bottom=c_fp,
        label='FN cost',     color=COLOR_FNR,  alpha=0.85)
ax7.bar(x_pos, c_vol, w,
        bottom=np.array(c_fp) + np.array(c_fn),
        label='Volume cost', color='#95a5a6', alpha=0.85)
ax7.set_xticks(x_pos); ax7.set_xticklabels(labels_bar, fontsize=10)
ax7.set_ylabel('Cost'); ax7.set_title('Cost Decomposition by Threshold')
ax7.legend(fontsize=9)

plt.tight_layout(rect=[0, 0, 1, 0.97])
plt.savefig('alert_opt_2d.png', dpi=150, bbox_inches='tight')
plt.show()
print("[Figure 1 — 2D dashboard saved]")

# ============================================================
# 7. 3D Surface: Optimal Threshold over (λ_fp, λ_fn) space
# ============================================================
fig2 = plt.figure(figsize=(20, 14))
fig2.suptitle('3D Sensitivity Analysis: Optimal Threshold Surface',
              fontsize=15, fontweight='bold')

# ── 3D surface ────────────────────────────────────────────
ax3d = fig2.add_subplot(2, 2, (1, 2), projection='3d')
surf = ax3d.plot_surface(
    LFP, LFN, THETA_SURF,
    cmap='viridis', alpha=0.88, edgecolor='none'
)
ax3d.set_xlabel('λ_fp (FP weight)',  labelpad=10)
ax3d.set_ylabel('λ_fn (FN weight)',  labelpad=10)
ax3d.set_zlabel('Optimal θ* (%)',    labelpad=10)
ax3d.set_title('Optimal Threshold θ* as a Function of Cost Weights')
fig2.colorbar(surf, ax=ax3d, shrink=0.5, label='θ* (%)')

# Mark current weights on surface
idx_fp = np.argmin(np.abs(lam_fp_vals - LAMBDA_FP))
idx_fn = np.argmin(np.abs(lam_fn_vals - LAMBDA_FN))
ax3d.scatter(
    [LAMBDA_FP], [LAMBDA_FN], [THETA_SURF[idx_fn, idx_fp]],
    color='red', s=150, zorder=10, label='Current weights'
)
ax3d.legend()

# ── 3D cost landscape for fixed weights ───────────────────
ax3d2 = fig2.add_subplot(2, 2, 3, projection='3d')
theta_g = np.linspace(40, 95, 120)
lam_fn_g = np.linspace(1, 20, 120)
TH, LN   = np.meshgrid(theta_g, lam_fn_g)
COST_LAND = (
    LAMBDA_FP  * fpr(TH)
    + LN * fnr(TH)
    + LAMBDA_VOL * alert_volume_per_hour(TH) / CHECKS_PER_HOUR
)
surf2 = ax3d2.plot_surface(TH, LN, COST_LAND, cmap='inferno', alpha=0.85, edgecolor='none')
ax3d2.set_xlabel('Threshold θ (%)', labelpad=8)
ax3d2.set_ylabel('λ_fn',            labelpad=8)
ax3d2.set_zlabel('Cost C(θ)',        labelpad=8)
ax3d2.set_title('Cost Landscape: θ vs λ_fn')
fig2.colorbar(surf2, ax=ax3d2, shrink=0.5, label='Cost')

# ── Heatmap of optimal threshold ──────────────────────────
ax_hm = fig2.add_subplot(2, 2, 4)
hm = ax_hm.contourf(LFP, LFN, THETA_SURF, levels=25, cmap='viridis')
fig2.colorbar(hm, ax=ax_hm, label='θ* (%)')
ax_hm.contour(LFP, LFN, THETA_SURF, levels=10, colors='white', linewidths=0.6, alpha=0.5)
ax_hm.scatter([LAMBDA_FP], [LAMBDA_FN], color='red', s=180,
              marker='*', label='Current', zorder=5)
ax_hm.set_xlabel('λ_fp'); ax_hm.set_ylabel('λ_fn')
ax_hm.set_title('Optimal Threshold Heatmap (λ_fp × λ_fn)')
ax_hm.legend()

plt.tight_layout()
plt.savefig('alert_opt_3d.png', dpi=150, bbox_inches='tight')
plt.show()
print("[Figure 2 — 3D sensitivity surface saved]")

# ============================================================
# 8. Summary statistics table
# ============================================================
print("\n" + "=" * 60)
print("  COMPARATIVE PERFORMANCE TABLE")
print("=" * 60)
print(f"{'Metric':<30} {'θ*='+str(round(THETA_OPT,1))+'%':>12} {'θ=70%':>10} {'θ=60%':>10}")
print("-" * 60)
for th, label in [(THETA_OPT, f'θ*={THETA_OPT:.1f}%'), (70.0, 'θ=70%'), (60.0, 'θ=60%')]:
    pass  # just for reference — table built below

rows = [
    ('FPR (%)',          lambda t: f'{fpr(t)*100:.2f}%'),
    ('FNR (%)',          lambda t: f'{fnr(t)*100:.2f}%'),
    ('Alerts/hr',        lambda t: f'{alert_volume_per_hour(t):.2f}'),
    ('Total cost C(θ)',  lambda t: f'{operational_cost(t):.4f}'),
]
for name, fn in rows:
    vals = [fn(t) for t in [THETA_OPT, 70.0, 60.0]]
    print(f'{name:<30} {vals[0]:>12} {vals[1]:>10} {vals[2]:>10}')
print("=" * 60)