import time
import numpy as np
import matplotlib.pyplot as plt
from itertools import permutations
import math
import openjij as oj

# ==================== 可調整參數區 ====================
# 問題規模設定
CITIES = 30         # 初始核心問題規模 (測試建議先用 12-15)
CITIES_BOUND = 10     # 問題規模的上下範圍 (單一數字控制：以 CITIES 為中心，例如 15±5，即 10~20)
RANDOM = False  
RANDOM_SEED = 42   
COORD_RANGE = (0.0, 10.0)  

# 【核心創新】：環境變異數控制與複雜度壓力測試參數
N_LIST = list(range(CITIES - CITIES_BOUND, CITIES + CITIES_BOUND + 1, 4))  # 自動產生測試的問題規模範圍
COORD_STD = 15.0  # 空間分佈變異數 (越大代表城市分佈越不均勻、越崎嶇)
RISK_STD = 10.0    # 擾動風險變異數 (越大代表某些路段特別危險，H-infinity 衝突極大)

# 演算法執行設定
NUM_RUNS = 10      # 正式比較的執行次數
BETA = 15.0        # 🌟 稍微調高 (原為 10.0)，讓低溫結冰得更紮實

# QUBO 參數
PENALTY = 500.0       # 🌟 從 3000 大幅降回 500 (解開高爾夫球場效應！)
BIG_PENALTY = 2000.0  # 🌟 起點約束降到 2000 即可

# TSP 求解參數
EXACT_LIMIT = 8    

# 魯棒優化參數 (H-infinity Robust QUBO)
USE_ROBUST = True  # 強制開啟，這是本論文的核心
GAMMA = 0.5        
SIGMA = 1.0        
ALPHA = 10.0       
# =====================================================

# ----------------- 基礎輔助函數 -----------------
def generate_distance_matrix(num_cities, random=True, seed=None, coord_range=None):
    if seed is None: seed = RANDOM_SEED
    if coord_range is None: coord_range = COORD_RANGE
    if not random: np.random.seed(seed)
    else: np.random.seed(None)
    low, high = coord_range
    coords = np.random.uniform(low, high, size=(num_cities, 2))
    D = np.zeros((num_cities, num_cities), dtype=float)
    for i in range(num_cities):
        for j in range(i + 1, num_cities):
            dist = np.linalg.norm(coords[i] - coords[j])
            D[i, j] = dist
            D[j, i] = dist
    return np.round(D, 2), coords

def generate_disturbance_matrix(num_cities, seed=None):
    if seed is None: seed = RANDOM_SEED
    np.random.seed(seed)
    F = np.random.rand(num_cities, num_cities)
    F = (F + F.T) / 2
    np.fill_diagonal(F, 0.0)
    return F

def idx_mtsp(k, i, p, N):
    return k * (N * N) + i * N + p

def build_robust_qubo(D, F, gamma=GAMMA, sigma=SIGMA, alpha=ALPHA, penalty=PENALTY, big_penalty=BIG_PENALTY):
    N = D.shape[0]
    Q = {}
    def addQ(u, v, w):
        if w == 0: return
        if u > v: u, v = v, u
        Q[(u, v)] = Q.get((u, v), 0.0) + w

    for k in range(2):
        for p in range(N):
            q = (p + 1) % N
            for i in range(N):
                for j in range(N):
                    dij, fij = D[i, j], F[i, j]
                    if dij == 0: continue
                    # 這裡就是 H-infinity 的核心風險項！完全保留著！
                    risk_term = (sigma / (gamma**2)) * (fij**2)
                    total_weight = dij + (alpha * risk_term)
                    u, v = idx_mtsp(k, i, p, N), idx_mtsp(k, j, q, N)
                    addQ(u, v, total_weight)

    for k in range(2):
        for p in range(N):
            vars_pos = [idx_mtsp(k, i, p, N) for i in range(N)]
            for u in vars_pos: addQ(u, u, -penalty)
            for a in range(N):
                for b in range(a+1, N): addQ(vars_pos[a], vars_pos[b], 2*penalty)

    for i in range(1, N):
        vars_city = [idx_mtsp(k, i, p, N) for k in range(2) for p in range(N)]
        for u in vars_city: addQ(u, u, -penalty)
        L = len(vars_city)
        for a in range(L):
            for b in range(a+1, L): addQ(vars_city[a], vars_city[b], 2*penalty)

    for k in range(2):
        for i in range(1, N): addQ(idx_mtsp(k, i, 0, N), idx_mtsp(k, i, 0, N), big_penalty)
        addQ(idx_mtsp(k, 0, 0, N), idx_mtsp(k, 0, 0, N), -big_penalty)

    return Q, N

# ---------------------------------------------------------
# 🌟 [新增] 產生平滑的 SQA 與 Pause SQA Schedule
# ---------------------------------------------------------
def get_smooth_sqa_schedule(beta, total_sweeps, num_steps=100):
    """產生標準 SQA 的平滑退火排程 (s 從 0.0 遞增到 1.0)"""
    schedule = []
    sweeps_per_step = max(1, int(total_sweeps / num_steps))
    for s in np.linspace(0.0, 1.0, num_steps):
        schedule.append([float(s), float(beta), int(sweeps_per_step)])
    return schedule

def get_pause_sqa_schedule(s_star, beta, total_anneal_sweeps, pause_sweeps, num_steps=100):
    """產生帶有 Pause 策略的平滑退火排程"""
    schedule = []
    
    # 階段 1：從 s=0.0 平滑退火到 s_star
    steps_p1 = max(1, int(num_steps * s_star))
    sweeps_p1 = max(1, int((total_anneal_sweeps / 2) / steps_p1))
    for s in np.linspace(0.0, s_star, steps_p1, endpoint=False):
        schedule.append([float(s), float(beta), int(sweeps_p1)])
        
    # 階段 2：在相變點 s_star 暫停 (Pause)，讓量子系統充分熱化
    schedule.append([float(s_star), float(beta), int(pause_sweeps)])
    
    # 階段 3：從 s_star 平滑退火到 s=1.0
    steps_p3 = max(1, num_steps - steps_p1)
    sweeps_p3 = max(1, int((total_anneal_sweeps / 2) / steps_p3))
    for s in np.linspace(s_star, 1.0, steps_p3):
        schedule.append([float(s), float(beta), int(sweeps_p3)])
        
    return schedule

def decode_slots(sample, N):
    slots = []
    for k in range(2):
        row = []
        for p in range(N):
            chosen = 0
            for i in range(N):
                if sample.get(idx_mtsp(k, i, p, N), 0) == 1:
                    chosen = i
                    break
            row.append(chosen)
        slots.append(row)
    return slots[0], slots[1]

def repair_routes_from_slots(u1_slots, u2_slots, N):
    count1, count2 = [0]*N, [0]*N
    for c in u1_slots: 
        if 0 <= c < N: count1[c] += 1
    for c in u2_slots: 
        if 0 <= c < N: count2[c] += 1
    assign1, assign2 = [], []
    for city in range(1, N):
        if count1[city] > count2[city]: assign1.append(city)
        elif count2[city] > count1[city]: assign2.append(city)
        else:
            if len(assign1) <= len(assign2): assign1.append(city)
            else: assign2.append(city)
    return sorted(assign1), sorted(assign2)

def uav_cost(path, D):
    if not path or len(path) == 0: return 0.0
    cost = D[0, path[0]]
    for i in range(len(path)-1): cost += D[path[i], path[i+1]]
    cost += D[path[-1], 0]
    return float(cost)

def uav_disturbance_energy(path, F, gamma=GAMMA, sigma=SIGMA):
    if not path or len(path) == 0: return 0.0
    risk = (sigma / (gamma**2)) * (F[0, path[0]]**2)
    for i in range(len(path)-1): risk += (sigma / (gamma**2)) * (F[path[i], path[i+1]]**2)
    risk += (sigma / (gamma**2)) * (F[path[-1], 0]**2)
    return float(risk)

def best_order_for_cities(cities, D, exact_limit=EXACT_LIMIT):
    cities = list(cities)
    if len(cities) <= 1: return cities, uav_cost(cities, D)
    if len(cities) <= exact_limit:
        best_perm, best_cost = None, float('inf')
        for perm in permutations(cities):
            c = uav_cost(list(perm), D)
            if c < best_cost: best_cost, best_perm = c, list(perm)
        return best_perm, best_cost
    else: # 簡單的 nearest neighbor heuristic 當作 fallback
        rem = set(cities)
        curr = list(rem)[0]
        rem.remove(curr)
        route = [curr]
        while rem:
            nxt = min(rem, key=lambda x: D[curr, x])
            route.append(nxt)
            rem.remove(nxt)
            curr = nxt
        return route, uav_cost(route, D)

def get_makespan_and_risk(sample, N, D, F):
    u1, u2 = decode_slots(sample, N)
    a1, a2 = repair_routes_from_slots(u1, u2, N)
    p1, c1 = best_order_for_cities(a1, D)
    p2, c2 = best_order_for_cities(a2, D)
    makespan = max(c1, c2)
    dist_energy = uav_disturbance_energy(p1, F) + uav_disturbance_energy(p2, F)
    # 確保回傳路徑，以便印出
    return makespan, dist_energy, p1, p2

# ============================
# 核心創新功能區 (Pause Strategy)
# ============================
def pilot_search_phase_transition(Q, N, D, F):
    print("\n🔍 [階段一] 前導偵查 (Pilot Search): 尋找最佳暫停點 s*")
    test_s_list = np.arange(0.2, 1.0, 0.1)
    variances = []
    makespans = []
    
    sampler = oj.SQASampler()
    
    for s in test_s_list:
        sched_var = [[s, BETA, 30]] 
        res_var = sampler.sample_qubo(Q, schedule=sched_var, num_reads=10)
        energies = [d.energy for d in res_var.data()]
        variances.append(np.var(energies))
            
        sched_pause = [[s, BETA, 20], [s, BETA, 40], [1.0, BETA, 20]]
        res_pause = sampler.sample_qubo(Q, schedule=sched_pause, num_reads=5)
        mk, _, p1, p2 = get_makespan_and_risk(res_pause.first.sample, N, D, F)
        makespans.append(mk)
        print(f"   測試 s={s:.1f} | 變異數: {variances[-1]:.0f} | Makespan: {mk:.2f}")

    best_s = test_s_list[np.argmin(makespans)]
    print(f"⭐ 鎖定最佳暫停點 s* = {best_s:.1f}")
    return test_s_list, variances, makespans, best_s

def test_pause_durations(Q, N, D, F, s_star):
    print(f"\n⏳ [階段二] 暫停時間測試 (在 s*={s_star:.1f} 進行熱化)")
    durations = [0, 20, 50, 100, 150]
    duration_makespans = []
    
    sampler = oj.SQASampler()
    for d in durations:
        if d == 0:
            sched = get_smooth_sqa_schedule(BETA, total_sweeps=1000)
        else:
            sched = get_pause_sqa_schedule(s_star, BETA, total_anneal_sweeps=1000, pause_sweeps=d)
            
        res = sampler.sample_qubo(Q, schedule=sched, num_reads=10)
        mk, _, p1, p2 = get_makespan_and_risk(res.first.sample, N, D, F)
        duration_makespans.append(mk)
        print(f"   暫停 {d} 步 -> Makespan: {mk:.2f}")
        
    return durations, duration_makespans

def run_comparative_evaluations(Q, N, D, F, s_star):
    print(f"\n⚔️ [階段三] 正式對決 (SA vs SQA vs Pause-SQA, {NUM_RUNS} runs)")
    
    results = {
        'sa': {'mk': [], 'risk': [], 'energy': []},
        'sqa': {'mk': [], 'risk': [], 'energy': []},
        'pause': {'mk': [], 'risk': [], 'energy': []}
    }
    
    sqa_sampler = oj.SQASampler()
    sa_sampler = oj.SASampler()
    
    # 🌟 核心修改：建立 1000 Sweeps 的退火排程
    TOTAL_SWEEPS = 1000
    PAUSE_SWEEPS = 500  # 讓 Pause-SQA 在相變點多停留 500 步熱化
    
    smooth_sched = get_smooth_sqa_schedule(beta=BETA, total_sweeps=TOTAL_SWEEPS, num_steps=100)
    pause_sched = get_pause_sqa_schedule(s_star=s_star, beta=BETA, total_anneal_sweeps=TOTAL_SWEEPS, pause_sweeps=PAUSE_SWEEPS, num_steps=100)

    # (選配) 在 pilot_search 和 test_pause 中也等比例放大 sweeps 
    # 例如 pilot_search 改成 [[s, BETA, 300]]
    
    for r in range(NUM_RUNS):
        # SA
        res_sa = sa_sampler.sample_qubo(Q, num_reads=10)
        mk, risk, p1, p2 = get_makespan_and_risk(res_sa.first.sample, N, D, F)
        results['sa']['mk'].append(mk)
        results['sa']['risk'].append(risk)
        results['sa']['energy'].append(res_sa.first.energy)
        print(f"      [SA 執行 {r+1}/{NUM_RUNS}] 路徑 1: [0, {', '.join(map(str, p1))}, 0], 路徑 2: [0, {', '.join(map(str, p2))}, 0] | 成本: {mk:.2f}")
        
        # Standard SQA
        t0 = time.time()
        res_std = sqa_sampler.sample_qubo(Q, schedule=smooth_sched, num_reads=10)
        mk, risk, p1, p2 = get_makespan_and_risk(res_std.first.sample, N, D, F)
        results['sqa']['mk'].append(mk)
        results['sqa']['risk'].append(risk)
        results['sqa']['energy'].append(res_std.first.energy)
        print(f"      [Std-SQA 執行 {r+1}/{NUM_RUNS}] 耗時: {time.time()-t0:.2f}s | 路徑 1: [0, {', '.join(map(str, p1))}, 0], 路徑 2: [0, {', '.join(map(str, p2))}, 0] | 成本: {mk:.2f}")
        
        # Pause SQA
        t0 = time.time()
        res_pause = sqa_sampler.sample_qubo(Q, schedule=pause_sched, num_reads=10)
        mk, risk, p1, p2 = get_makespan_and_risk(res_pause.first.sample, N, D, F)
        results['pause']['mk'].append(mk)
        results['pause']['risk'].append(risk)
        results['pause']['energy'].append(res_pause.first.energy)
        print(f"      [Pause-SQA 執行 {r+1}/{NUM_RUNS}] 耗時: {time.time()-t0:.2f}s | 路徑 1: [0, {', '.join(map(str, p1))}, 0], 路徑 2: [0, {', '.join(map(str, p2))}, 0] | 成本: {mk:.2f}")
        
    return results

# ============================
# 學術圖表繪製 (Academic Visualization)
# ============================
def create_basic_distribution_charts(results):
    """基本分佈分析圖"""
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    fig.suptitle(f"Algorithm Performance Distributions (N={CITIES})", fontsize=16, fontweight='bold')
    
    labels = ['SA', 'Standard-SQA', 'Pause-SQA']
    colors = ['gray', 'steelblue', 'firebrick']
    
    # 1. Makespan 分布
    parts1 = axes[0].violinplot([results['sa']['mk'], results['sqa']['mk'], results['pause']['mk']], showmeans=True)
    axes[0].set_title("Makespan (Distance Cost) Distribution")
    axes[0].set_ylabel("Distance Cost")
    
    # 2. Risk 分布
    parts2 = axes[1].violinplot([results['sa']['risk'], results['sqa']['risk'], results['pause']['risk']], showmeans=True)
    axes[1].set_title(r"$H_\infty$ Disturbance Risk Distribution")
    axes[1].set_ylabel("Risk Penalty")
    
    # 3. Energy 分布
    parts3 = axes[2].violinplot([results['sa']['energy'], results['sqa']['energy'], results['pause']['energy']], showmeans=True)
    axes[2].set_title("QUBO System Raw Energy")
    axes[2].set_ylabel("System Energy")
    
    # 設定提琴圖外觀
    for i, ax in enumerate(axes):
        ax.set_xticks([1, 2, 3])
        ax.set_xticklabels(labels, fontsize=11)
        ax.grid(alpha=0.3, axis='y')
        parts = [parts1, parts2, parts3][i]
        for pc, color in zip(parts['bodies'], colors):
            pc.set_facecolor(color)
            pc.set_alpha(0.7)

    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.savefig("distribution_comparison.png", dpi=300)
    print("\n🎨 圖表已繪製: distribution_comparison.png")
    # plt.show() (Moved to the end)

def create_academic_charts(s_list, variances, makespans, best_s, durations, dur_makespans, results):
    """原有的圖表二：Pause Strategy 深度分析圖"""
    fig = plt.figure(figsize=(16, 10))
    fig.suptitle(f"Quantum Pause Strategy in H-infinity Robust mTSP (N={CITIES})", fontsize=16, fontweight='bold')
    
    # Chart 1: Phase Transition & U-Valley
    ax1 = plt.subplot(2, 2, 1)
    color1, color2 = 'tab:blue', 'tab:red'
    ax1.set_title("Chart 1: Phase Transition & U-Valley", fontsize=12)
    ax1.set_xlabel("Annealing Parameter $s$")
    ax1.set_ylabel("Trotter Energy Variance", color=color1)
    ax1.plot(s_list, variances, marker='o', color=color1, linewidth=2)
    ax1.tick_params(axis='y', labelcolor=color1)
    
    ax1_twin = ax1.twinx()
    ax1_twin.set_ylabel("Makespan", color=color2)
    ax1_twin.plot(s_list, makespans, marker='s', color=color2, linestyle='--', linewidth=2)
    ax1_twin.tick_params(axis='y', labelcolor=color2)
    ax1.axvline(x=best_s, color='green', linestyle=':', linewidth=2)
    ax1.grid(alpha=0.3)
    
    # Subplot 2 (was Chart 4): Effect of Thermalization Time
    ax2 = plt.subplot(2, 2, 2)
    ax2.set_title(f"Chart 2: Effect of Thermalization Time at $s^*={best_s:.1f}$", fontsize=12)
    ax2.plot(durations, dur_makespans, marker='D', color='darkorange', linewidth=2.5, markersize=8)
    ax2.set_xlabel("Pause Duration (Sweeps)")
    ax2.set_ylabel("Optimized Makespan")
    ax2.grid(alpha=0.3, linestyle='--')

    # Subplot 3 (was Chart 2): Pareto Front
    ax3 = plt.subplot(2, 2, 3)
    ax3.set_title("Chart 3: Robustness Pareto Front", fontsize=12)
    ax3.scatter(results['sqa']['mk'], results['sqa']['risk'], c='steelblue', alpha=0.7, s=80, label='Standard-SQA')
    ax3.scatter(results['pause']['mk'], results['pause']['risk'], c='firebrick', alpha=0.9, s=120, marker='*', label='Pause-SQA')
    ax3.set_xlabel("Makespan")
    ax3.set_ylabel("Disturbance Risk")
    ax3.legend()
    ax3.grid(alpha=0.3, linestyle='--')

    # Subplot 4 (was Chart 3): Reliability Violin Plot
    ax4 = plt.subplot(2, 2, 4)
    ax4.set_title("Chart 4: Final Makespan Comparison", fontsize=12)
    parts = ax4.violinplot([results['sqa']['mk'], results['pause']['mk']], showmeans=True)
    parts['bodies'][0].set_facecolor('steelblue')
    parts['bodies'][1].set_facecolor('firebrick')
    for body in parts['bodies']: body.set_alpha(0.7)
    ax4.set_xticks([1, 2])
    ax4.set_xticklabels(['Standard-SQA', 'Pause-SQA'])
    
    sa_mean = np.mean(results['sa']['mk'])
    ax4.axhline(y=sa_mean, color='gray', linestyle=':', linewidth=2, label=f'SA Mean: {sa_mean:.1f}')
    ax4.legend()
    ax4.grid(alpha=0.3, axis='y')

    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.savefig("academic_pause_strategy.png", dpi=300)
    print("🎨 原有圖表已繪製: academic_pause_strategy.png")
    # plt.show() (Moved to the end)

# ============================
# 複雜度與環境控制功能 (新增)
# ============================
def inject_hidden_shortcut(D, F):
    """
    注入隱藏的黃金捷徑 (Hidden Shortcut)：
    創造一條極度狹窄且成本極低的路線 (例如 4 -> 5 -> 6 -> 7)。
    如果沒照順序走，成本會極高；一旦走對，總成本會大幅下降。
    """
    N = D.shape[0]
    if N < 8:
        return D, F  # 城市太少無法做長捷徑

    # 定義黃金路線的節點 (假設是 4, 5, 6, 7)
    golden_nodes = [4, 5, 6, 7]

    # 1. 築起高牆：先把這幾個城市之間的所有連線，都變成極高成本 (阻止 SA 亂湊)
    for i in golden_nodes:
        for j in golden_nodes:
            if i != j:
                D[i, j] = 8.0
                F[i, j] = 0.8  # 高風險

    # 2. 挖出深谷：只開通 4->5, 5->6, 6->7 這條唯一且完美的捷徑
    for i in range(len(golden_nodes) - 1):
        n1 = golden_nodes[i]
        n2 = golden_nodes[i + 1]
        
        # 距離極短，風險為 0
        D[n1, n2] = 0.05
        D[n2, n1] = 0.05
        F[n1, n2] = 0.0
        F[n2, n1] = 0.0

    print("\n✨ [隱藏捷徑已佈署] 演算法將挑戰尋找極狹窄的黃金路線 (4->5->6->7)！")
    return D, F

def inject_deceptive_trap(D, F, alpha=10.0, gamma=0.5, sigma=1.0):
    """
    注入欺騙性陷阱 (Deceptive Trap)：
    改造矩陣，創造一條「看似完美的捷徑」，測試演算法是否會陷入局部陷阱。
    """
    N = D.shape[0]
    if N < 4:
        return D, F  # 城市太少無法做陷阱

    # 定義陷阱節點
    trap_start = 1
    trap_end = 2
    safe_detour = 3

    # ==========================================
    # 陷阱 1：致命捷徑 (距離極度誘人，但風險爆表)
    # ==========================================
    D[trap_start, trap_end] = 0.1
    D[trap_end, trap_start] = 0.1
    F[trap_start, trap_end] = 0.99
    F[trap_end, trap_start] = 0.99

    # ==========================================
    # 陷阱 2：安全繞路 (距離較遠，但完全無風險)
    # ==========================================
    # 路線: trap_start -> safe_detour -> trap_end
    D[trap_start, safe_detour] = 4.0
    D[safe_detour, trap_start] = 4.0
    D[safe_detour, trap_end] = 4.0
    D[trap_end, safe_detour] = 4.0

    F[trap_start, safe_detour] = 0.01
    F[safe_detour, trap_start] = 0.01
    F[safe_detour, trap_end] = 0.01
    F[trap_end, safe_detour] = 0.01

    
    # 計算並印出真實成本，讓您確認陷阱是否成立
    # 真實成本 = 距離 + alpha * (sigma / gamma^2) * F^2
    risk_multiplier = alpha * (sigma / (gamma**2))
    
    shortcut_cost = 0.1 + risk_multiplier * (0.99**2)
    detour_cost = (4.0 + risk_multiplier * (0.01**2)) + (4.0 + risk_multiplier * (0.01**2))
    
    print(f"   ➤ 致命捷徑 (1->2) 表觀距離: 0.1 | 隱藏真實成本: {shortcut_cost:.2f}")
    print(f"   ➤ 安全繞路 (1->3->2) 表觀距離: 8.0 | 隱藏真實成本: {detour_cost:.2f}")
    
    return D, F

def generate_controlled_matrices(num_cities, coord_std, risk_std, seed=42):
    """根據指定的標準差(變異數)生成 D 與 F 矩陣，精準控制問題複雜度"""
    np.random.seed(seed)
    
    # 1. 生成距離矩陣 D (使用常態分佈控制空間聚集度)
    coords = np.random.normal(loc=0.0, scale=coord_std, size=(num_cities, 2))
    D = np.zeros((num_cities, num_cities), dtype=float)
    for i in range(num_cities):
        for j in range(i + 1, num_cities):
            dist = np.linalg.norm(coords[i] - coords[j])
            D[i, j] = D[j, i] = dist
            
    # 2. 生成擾動矩陣 F (使用常態分佈控制風險極端值)
    F = np.abs(np.random.normal(loc=1.0, scale=risk_std, size=(num_cities, num_cities)))
    F = (F + F.T) / 2  # 確保對稱
    np.fill_diagonal(F, 0.0)
    
    return np.round(D, 2), np.round(F, 2)

def decode_and_eval(sample, N, D):
    # 快速解碼與貪婪/暴力求成本
    u1, u2 = [], []
    for p in range(N):
        for i in range(N):
            if sample.get(idx_mtsp(0, i, p, N), 0) == 1: u1.append(i)
            if sample.get(idx_mtsp(1, i, p, N), 0) == 1: u2.append(i)
    
    # 簡單修復 (去除重複，補齊缺失)
    c1, c2 = set(u1), set(u2)
    a1, a2 = [], []
    for city in range(1, N):
        if city in c1 and city not in c2: a1.append(city)
        elif city in c2 and city not in c1: a2.append(city)
        else:
            if len(a1) <= len(a2): a1.append(city)
            else: a2.append(city)
            
    # 評估成本
    def eval_cost(cities):
        if not cities: return 0
        if len(cities) <= EXACT_LIMIT:
            bc = float('inf')
            for p in permutations(cities):
                c = D[0, p[0]] + sum(D[p[i], p[i+1]] for i in range(len(p)-1)) + D[p[-1], 0]
                bc = min(bc, c)
            return bc
        # 貪婪
        rem, path, curr = set(cities), [], 0
        while rem:
            nxt = min(rem, key=lambda x: D[curr, x])
            path.append(nxt)
            rem.remove(nxt)
            curr = nxt
        return D[0, path[0]] + sum(D[path[i], path[i+1]] for i in range(len(path)-1)) + D[path[-1], 0]

    return max(eval_cost(a1), eval_cost(a2))

def run_complexity_scaling():
    print(f"\n🚀 開始複雜度壓力測試...")
    print(f"⚙️ 空間變異數(COORD_STD)={COORD_STD}, 風險變異數(RISK_STD)={RISK_STD}")
    
    std_means, pause_means = [], []
    std_errs, pause_errs = [], []
    
    sampler = oj.SQASampler()
    
    for N in N_LIST:
        print(f"\n📊 測試規模 N={N} ...")
        D, F = generate_controlled_matrices(N, COORD_STD, RISK_STD)
        
        # 動態計算懲罰值
        max_d = np.max(D)
        max_f = np.max(F) if F is not None else 0
        max_possible_cost = max_d + max_f
        dyn_penalty = max_possible_cost * 12.0
        dyn_big_penalty = dyn_penalty * 5.0
        
        Q, _ = build_robust_qubo(D, F, penalty=dyn_penalty, big_penalty=dyn_big_penalty)
        
        std_results, pause_results = [], []
        
        # 為了節省時間，在壓力測試中我們使用固定的經驗暫停點 s=0.7
        # 總步數統一為 200 步
        std_sched = get_smooth_sqa_schedule(BETA, total_sweeps=200, num_steps=20)
        pause_sched = get_pause_sqa_schedule(0.7, BETA, total_anneal_sweeps=200, pause_sweeps=100, num_steps=30)
        
        # 這裡設定每次 N 跑的次數，與原本 NUM_RUNS 獨立，建議 10 次
        runs_for_scale = 10 
        for r in range(runs_for_scale):
            # Standard SQA
            res_std = sampler.sample_qubo(Q, schedule=std_sched, num_reads=5)
            std_results.append(decode_and_eval(res_std.first.sample, N, D))
            
            # Pause SQA
            res_pause = sampler.sample_qubo(Q, schedule=pause_sched, num_reads=5)
            pause_results.append(decode_and_eval(res_pause.first.sample, N, D))
            
        mean_std, err_std = np.mean(std_results), np.std(std_results)
        mean_pause, err_pause = np.mean(pause_results), np.std(pause_results)
        
        std_means.append(mean_std)
        std_errs.append(err_std)
        pause_means.append(mean_pause)
        pause_errs.append(err_pause)
        
        print(f"   Standard-SQA 平均: {mean_std:.2f} ± {err_std:.2f}")
        print(f"   Pause-SQA 平均: {mean_pause:.2f} ± {err_pause:.2f}")
        
    return std_means, std_errs, pause_means, pause_errs

def plot_crossover(std_means, std_errs, pause_means, pause_errs):
    plt.figure(figsize=(10, 6))
    plt.title(f"Complexity Scaling: Algorithm Performance vs Problem Size (N)\n(Coord Std={COORD_STD}, Risk Std={RISK_STD})", fontsize=14, fontweight='bold')
    
    # 畫線與誤差棒
    plt.errorbar(N_LIST, std_means, yerr=std_errs, fmt='-o', color='steelblue', 
                 linewidth=2.5, capsize=5, markersize=8, label='Standard-SQA')
    plt.errorbar(N_LIST, pause_means, yerr=pause_errs, fmt='-s', color='firebrick', 
                 linewidth=2.5, capsize=5, markersize=8, label='Pause-SQA')
    
    plt.xlabel("Problem Scale (Number of Cities $N$)", fontsize=12)
    plt.ylabel("Optimized Makespan (Distance Cost)", fontsize=12)
    plt.xticks(N_LIST)
    plt.grid(alpha=0.4, linestyle='--')
    plt.legend(fontsize=11)
    
    plt.tight_layout()
    plt.savefig("complexity_crossover.png", dpi=300)
    print("\n🎨 交叉點分析圖已儲存為: complexity_crossover.png")
    # plt.show() (Moved to the end)

# ============================
# 主程式執行入口
# ============================
if __name__ == "__main__":
    D, city_coords = generate_distance_matrix(CITIES, random=RANDOM)
    F = generate_disturbance_matrix(CITIES)
    
    # 1. 佈署死亡陷阱 (測試演算法是否會被騙)
    D, F = inject_deceptive_trap(D, F, alpha=ALPHA, gamma=GAMMA, sigma=SIGMA)
    
    # 2. 🌟 佈署黃金捷徑 (測試演算法能否找到狹窄深谷)
    D, F = inject_hidden_shortcut(D, F)
    
    # ======== 🌟 新增：動態自適應懲罰值 ========
    # 找出矩陣中最遙遠的距離與最大的擾動
    max_d = np.max(D)
    max_f = np.max(F) if F is not None else 0
    max_possible_cost = max_d + max_f
    
    # 根據最大可能成本，動態設定懲罰值 (通常設為 10 倍 ~ 15 倍最穩定)
    PENALTY = max_possible_cost * 12.0
    BIG_PENALTY = PENALTY * 5.0
    
    print(f"🔧 [自動調校] 偵測到單趟最大成本: {max_possible_cost:.2f}")
    print(f"🔧 [自動調校] 動態 PENALTY 設為: {PENALTY:.2f}, BIG_PENALTY 設為: {BIG_PENALTY:.2f}")
    # ==========================================

    print("==================================================")
    print(f"🚀 Quantum Pause Strategy Optimization (mTSP N={CITIES})")
    print("==================================================")
    
    print("\n[初始化] 距離矩陣 D:")
    print(D)
    print("\n[初始化] 擾動矩陣 F:")
    print(F)
    
    Q, _ = build_robust_qubo(D, F, penalty=PENALTY, big_penalty=BIG_PENALTY)
    
    # 階段 1：找相變點
    s_list, variances, pilot_mks, best_s = pilot_search_phase_transition(Q, CITIES, D, F)
    
    # 階段 2：測試暫停時間
    durations, dur_mks = test_pause_durations(Q, CITIES, D, F, best_s)
    
    # 階段 3：正式評估 (一次抓齊所有數據)
    results = run_comparative_evaluations(Q, CITIES, D, F, best_s)
    
    # 階段 4：複雜度壓力測試與交叉點分析
    s_m, s_e, p_m, p_e = run_complexity_scaling()
    
    print("\n🏁 所有運算完成！開始繪製圖表...")
    
    # 繪製圖表一：三大基礎分布對比
    create_basic_distribution_charts(results)
    
    # 繪製圖表二：Pause SQA 深度分析
    create_academic_charts(s_list, variances, pilot_mks, best_s, durations, dur_mks, results)
    
    # 繪製圖表三：複雜度壓力測試與交叉點分析
    plot_crossover(s_m, s_e, p_m, p_e)
    
    print("🏁 所有實驗與圖表繪製完成！")
    plt.show()