TPS-RAS-repeat.py

import numpy as np
import time
import random
import copy
import torch
import torch.nn as nn
import torch.optim as optim
import collections
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from scipy import stats

# ==========================================
# 1. System Model (Env)
# ==========================================

class FogNode:
    def __init__(self, id, tier, cpu, ram, bw):
        self.id = id
        self.tier = tier  # 1: Edge, 2: Fog, 3: Cloud
        self.capacity = np.array([cpu, ram, bw], dtype=np.float64)
        self.current_free = self.capacity.copy()

    def reset(self):
        self.current_free = self.capacity.copy()

class Task:
    def __init__(self, id, req_cpu, req_ram, req_bw, arrival_time):
        self.id = id
        self.req = np.array([req_cpu, req_ram, req_bw], dtype=np.float64)
        self.arrival_time = arrival_time

def generate_heterogeneous_env(num_tasks=100):
    nodes = []
    # Tier 1: Edge (30)
    for i in range(30):
        nodes.append(FogNode(len(nodes), 1, 
                             random.uniform(500, 1000), 
                             random.uniform(512, 1024), 100))
    # Tier 2: Fog (8)
    for i in range(8):
        nodes.append(FogNode(len(nodes), 2, 
                             random.uniform(3000, 6000), 
                             random.uniform(8096, 16384), 1000))
    # Tier 3: Cloud (1)
    nodes.append(FogNode(len(nodes), 3, 1e8, 1e8, 1e6)) 

    tasks = []
    current_time = 0
    for i in range(num_tasks):
        inter_arrival = random.expovariate(1.0 / 10.0) 
        current_time += inter_arrival
        if random.random() < 0.7: 
            req = [random.uniform(50, 300), random.uniform(20, 100), random.uniform(1, 5)]
        else: 
            req = [random.uniform(200, 800), random.uniform(500, 2000), random.uniform(20, 50)]
        tasks.append(Task(i, req[0], req[1], req[2], current_time))
        
    return nodes, tasks

def plot_box_results(all_data_list, workloads):
    """
    all_data_list: List of 'data' dicts for each workload [data_100, data_300, data_500]
    workloads: [100, 300, 500]
    """
    # 1. Convert Data to DataFrame for Seaborn
    records = []
    for n_task, data_dict in zip(workloads, all_data_list):
        for algo, metrics in data_dict.items():
            # metrics['Util'] is a list of 30 values
            for i in range(len(metrics['Util'])):
                records.append({
                    'Tasks': n_task,
                    'Algorithm': algo,
                    'Fog Util (%)': metrics['Util'][i],
                    'NetCost (Hops)': metrics['NetCost'][i],
                    'Time (s)': metrics['Time'][i]
                })
    
    df = pd.DataFrame(records)

    # 2. Plot Setup
    sns.set(style="whitegrid", font_scale=1.1)
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    
    # Define colors (TPS-RAS distinct)
    palette = {"Edgewards": "#95a5a6", "KubeDSM": "#34495e", "DRL": "#e74c3c", "TPS-RAS": "#2ecc71"}

    # Subplot 1: Fog Utilization
    sns.boxplot(x="Tasks", y="Fog Util (%)", hue="Algorithm", data=df, ax=axes[0], palette=palette, linewidth=1.2, fliersize=3)
    axes[0].set_title("Resource Utilization Stability")
    axes[0].get_legend().remove()

    # Subplot 2: Network Cost
    sns.boxplot(x="Tasks", y="NetCost (Hops)", hue="Algorithm", data=df, ax=axes[1], palette=palette, linewidth=1.2, fliersize=3)
    axes[1].set_title("Network Cost Distribution")
    axes[1].get_legend().remove()

    # Subplot 3: Decision Time (Log Scale if needed, but linear usually fine for Comparison)
    sns.boxplot(x="Tasks", y="Time (s)", hue="Algorithm", data=df, ax=axes[2], palette=palette, linewidth=1.2, fliersize=3)
    axes[2].set_title("Decision Latency (Speed & Jitter)")
    
    # Adjust Legend
    handles, labels = axes[2].get_legend_handles_labels()
    axes[2].legend(handles=handles, labels=labels, loc='upper left', bbox_to_anchor=(1, 1))

    plt.tight_layout()
    plt.savefig("Experiment_BoxPlots.png", dpi=300)
    print("\n[Info] Box plots saved to 'Experiment_BoxPlots.png'")
    plt.show()

# ==========================================
# 2. DRL Agent (Advantage Actor-Critic [26])
# ==========================================


class A2CNetwork(nn.Module):
    def __init__(self, input_dim, output_dim):
        super(A2CNetwork, self).__init__()
        # 公共特征提取层
        self.fc = nn.Sequential(
            nn.Linear(input_dim, 256),
            nn.ReLU(),
            nn.Linear(256, 128),
            nn.ReLU()
        )
        # Actor Head: 输出动作概率
        self.actor = nn.Sequential(
            nn.Linear(128, output_dim),
            nn.Softmax(dim=-1)
        )
        # Critic Head: 输出状态价值 V(s)
        self.critic = nn.Linear(128, 1)
        
    def forward(self, x):
        feat = self.fc(x)
        probs = self.actor(feat)
        value = self.critic(feat)
        return probs, value

class A2CAgent:
    def __init__(self, num_nodes):
        self.num_nodes = num_nodes
        self.state_dim = 3 + num_nodes * 2 
        self.action_dim = num_nodes
        self.device = torch.device("cpu")
        self.net = A2CNetwork(self.state_dim, self.action_dim).to(self.device)
        self.optimizer = optim.Adam(self.net.parameters(), lr=0.003)
        self.gamma = 0.95
        
    def select_action(self, state, valid_mask):
        state_t = torch.FloatTensor(state).unsqueeze(0).to(self.device)
        probs, _ = self.net(state_t)
        
        # 掩盖无效动作 (Masking)
        probs_np = probs.detach().cpu().numpy()[0]
        probs_np[~valid_mask] = 0.0
        
        # 重新归一化
        if np.sum(probs_np) > 0:
            probs_np /= np.sum(probs_np)
        else:
            # 如果 mask 掉后全为 0 (极罕见)，随机选一个 valid 的
            valid_indices = np.where(valid_mask)[0]
            action = np.random.choice(valid_indices)
            return action, torch.tensor([1.0/len(valid_indices)])
            
        # 根据概率采样 (Stochastic Policy)
        action = np.random.choice(self.action_dim, p=probs_np)
        return action, probs[0, action] # 返回动作和该动作的概率log值(用于求导)

    def update(self, state, action, reward, next_state, done):
        # A2C 是在线学习 (On-policy)，通常单步或小批次更新
        # 这里为了简化嵌入到 Simulator，使用单步更新 (One-step Actor-Critic)
        
        state_t = torch.FloatTensor(state).unsqueeze(0).to(self.device)
        next_state_t = torch.FloatTensor(next_state).unsqueeze(0).to(self.device)
        reward_t = torch.FloatTensor([reward]).to(self.device)
        action_t = torch.LongTensor([action]).view(1, -1).to(self.device)
        
        # 前向传播
        probs, value = self.net(state_t)
        _, next_value = self.net(next_state_t)
        
        # 计算 TD Target 和 Advantage
        target = reward_t + self.gamma * next_value.detach() * (1 - int(done))
        advantage = target - value
        
        # Loss 计算
        # Critic Loss: MSE(V(s), Target)
        critic_loss = advantage.pow(2).mean()
        
        # Actor Loss: -log(prob) * Advantage
        # 重新获取当前动作的 log_prob
        dist = torch.distributions.Categorical(probs)
        log_prob = dist.log_prob(action_t)
        actor_loss = -(log_prob * advantage.detach()).mean()
        
        total_loss = actor_loss + critic_loss
        
        self.optimizer.zero_grad()
        total_loss.backward()
        self.optimizer.step()

# ==========================================
# 3. Simulator & Algorithms
# ==========================================

class Simulator:
    def __init__(self, nodes, tasks):
        self.nodes = copy.deepcopy(nodes)
        self.tasks = tasks
    
    def reset_nodes(self):
        for n in self.nodes: n.reset()

    def evaluate(self, allocation, exec_time):
        success = 0
        total_hops = 0
        fog_nodes = [n for n in self.nodes if n.tier in [1, 2]]
        caps = np.sum([n.capacity for n in fog_nodes], axis=0)
        used = np.zeros(3)
        
        for t_id, n_id in allocation.items():
            if n_id is not None:
                success += 1
                n = self.nodes[n_id]
                total_hops += n.tier * (1 if n.tier < 3 else 5)
                if n.tier in [1, 2]:
                    used += self.tasks[t_id].req
        
        util_avg = np.mean(used / (caps + 1e-9)) if success > 0 else 0
        
        return {
            "Util%": util_avg * 100,
            "NetCost": total_hops / (success + 1e-9),
            "Time(s)": exec_time
        }

    # --- Baseline: Edgewards ---
    def run_edgewards(self):
        start = time.time()
        self.reset_nodes()
        alloc = {}
        sorted_nodes = sorted(range(len(self.nodes)), key=lambda i: self.nodes[i].tier)
        for t in self.tasks:
            assigned = None
            for n_idx in sorted_nodes:
                if np.all(t.req <= self.nodes[n_idx].current_free):
                    self.nodes[n_idx].current_free -= t.req
                    assigned = n_idx
                    break
            alloc[t.id] = assigned
        return self.evaluate(alloc, time.time() - start)

    # --- Baseline: Kubernetes ---
    # --- Advanced Baseline: KubeDSM (Dynamic Scheduling & Migration) ---
    def run_kubedsm(self):
        """
        KubeDSM [21]: 扩展的 K8s 调度器。 
        特性：
        1. 优先使用 LeastAllocated (K8s Default) 进行调度。
        2. [核心区别] 当任务无法调度时，触发 'Migration' 机制：
           - 扫描节点，寻找可以通过'挪走小任务'来腾出空间的节点。
           - 执行模拟迁移（增加额外的迁移延迟惩罚）。
        """
        start = time.time()
        self.reset_nodes()
        alloc = {}
        
        # 定义迁移带来的额外开销 (模拟容器冷迁移: 停止-传输-启动)
        # 假设迁移一个任务平均带来 200ms (0.2s) 的系统阻塞/延迟惩罚
        MIGRATION_PENALTY = 0.2 
        total_migration_count = 0

        # 保存已分配的任务位置，方便迁移时查询 {task_id: node_idx}
        task_placement = {} 

        for t in self.tasks:
            # --- 步骤 1: 标准 K8s 调度 (LeastAllocated) ---
            best_node = -1
            max_score = -1
            
            # 尝试找到能直接放下的节点
            candidates = []
            for n_idx, n in enumerate(self.nodes):
                if np.all(t.req <= n.current_free):
                    # LeastAllocated Score: (Free / Cap) 的平均值
                    # 分数越高，剩余比例越大（负载均衡）
                    score = np.mean(n.current_free / (n.capacity + 1e-9))
                    if score > max_score:
                        max_score = score
                        best_node = n_idx
            
            # --- 步骤 2: 如果直接调度成功 ---
            if best_node != -1:
                self.nodes[best_node].current_free -= t.req
                alloc[t.id] = best_node
                task_placement[t.id] = best_node
            
            # --- 步骤 3: [KubeDSM 特性] 触发迁移机制 ---
            else:
                migration_success = False
                
                # 遍历所有节点，寻找"潜在宿主"
                # 条件：节点的总容量够放这个任务，只是现在被占满了
                potential_nodes = [n for n in self.nodes if np.all(n.capacity >= t.req)]
                
                # 按当前剩余资源从大到小排序，优先尝试容易腾出空间的节点
                # 简单启发式：找那个"稍微努努力(迁移几个小任务)"就能放下的节点
                potential_nodes.sort(key=lambda n: np.mean(n.current_free), reverse=True)
                
                for p_node in potential_nodes:
                    # 找到该节点上当前运行的所有任务
                    running_tasks = [tid for tid, nid in task_placement.items() if nid == p_node.id]
                    
                    # 模拟：尝试把该节点上的任务迁移出去
                    # 简单策略：贪心尝试迁移小任务，直到空间足够
                    
                    # 备份节点状态，以便回滚
                    original_free = p_node.current_free.copy()
                    
                    migrated_tasks = []
                    space_cleared = p_node.current_free.copy()
                    
                    # 寻找可以迁移的候选任务 (简单的 First-Fit 尝试)
                    for rt_id in running_tasks:
                        # 如果空间已经够了，停止寻找
                        if np.all(space_cleared >= t.req):
                            break
                        
                        # 尝试把 rt_id 迁移到 *其他* 节点
                        # 这里的 rt_task 我们需要从 self.tasks 里找
                        rt_task = self.tasks[rt_id]
                        
                        # 找一个"备胎"节点 (不能是当前 p_node)
                        target_node = -1
                        for other_n in self.nodes:
                            if other_n.id != p_node.id and np.all(other_n.current_free >= rt_task.req):
                                target_node = other_n.id
                                break
                        
                        if target_node != -1:
                            # 模拟迁移：
                            # 1. 从 p_node 移除 (逻辑上增加 space_cleared)
                            space_cleared += rt_task.req
                            # 2. 记录这个迁移操作
                            migrated_tasks.append((rt_id, target_node, rt_task))
                    
                    # 检查是否腾出了足够空间
                    if np.all(space_cleared >= t.req):
                        # 执行真正的迁移提交
                        # 1. 把任务从 p_node 移走，放到 target_node
                        for mt_id, target_idx, mt_task in migrated_tasks:
                            # 更新 p_node (释放资源)
                            self.nodes[p_node.id].current_free += mt_task.req
                            # 更新 target_node (占用资源)
                            self.nodes[target_idx].current_free -= mt_task.req
                            # 更新索引
                            task_placement[mt_id] = target_idx
                            alloc[mt_id] = target_idx
                            
                            total_migration_count += 1
                        
                        # 2. 把当前大任务 t 放到 p_node
                        self.nodes[p_node.id].current_free -= t.req
                        alloc[t.id] = p_node.id
                        task_placement[t.id] = p_node.id
                        
                        migration_success = True
                        break # 任务 t 安排好了，跳出节点循环
                
                if not migration_success:
                    alloc[t.id] = None # 实在没办法，迁移也腾不出

        # 计算总耗时：基础调度时间 + 迁移惩罚
        exec_time = (time.time() - start) + (total_migration_count * MIGRATION_PENALTY)
        
        return self.evaluate(alloc, exec_time)

    # --- Rival: A2C (Up-to-Date DRL) ---
    def run_drl(self):
        start = time.time()
        self.reset_nodes()
        alloc = {}
        # 替换为 A2C Agent
        agent = A2CAgent(len(self.nodes))
        
        # 为了 A2C 更新，我们需要上一步的状态
        prev_state = None
        prev_action = None
        
        for i, t in enumerate(self.tasks):
            # 1. State Construction
            current_frees = np.array([n.current_free[:2] for n in self.nodes]).flatten()
            state = np.concatenate([t.req, current_frees])
            
            # 2. Validity Mask
            valid = np.array([np.all(t.req <= n.current_free) for n in self.nodes])
            if not np.any(valid): 
                alloc[t.id] = None
                # 任务失败，给予惩罚并更新上一步
                if prev_state is not None:
                     agent.update(prev_state, prev_action, -10, state, True)
                prev_state = None
                continue
                
            # 3. Action Selection (Probabilistic)
            action, _ = agent.select_action(state, valid)
            
            # 4. Reward for PREVIOUS step (Delayed Reward)
            # 在单步模拟中，我们假设 Reward 是即时的 (Immediate Reward)
            # 如果是上一步的动作导致了现在的状态... 这里为了简化，直接算当前步 Reward
            
            n = self.nodes[action]
            n.current_free -= t.req
            alloc[t.id] = action
            
            tier_penalty = {1: 0, 2: 2, 3: 10}[n.tier]
            reward = 10 - tier_penalty
            
            # 5. Online Update (A2C learns step-by-step)
            # 注意：A2C 需要 (s, a, r, s')。我们用当前 s 和 模拟的 s' (扣除资源后)
            # 模拟 Next State
            next_current_frees = np.array([n.current_free[:2] for n in self.nodes]).flatten()
            if i < len(self.tasks) - 1:
                next_req = self.tasks[i+1].req
            else:
                next_req = np.zeros(3) # No more tasks
            next_state = np.concatenate([next_req, next_current_frees])
            
            agent.update(state, action, reward, next_state, False)
            
        return self.evaluate(alloc, time.time() - start)

    # --- Proposed: TPS-RAS ---
    def compute_scarcity_weight_dynamic(self, N_cap, N_free, tiers, eta=5.0):
            """
            Fix: Node-Level Scarcity Perception
            Instead of a global average, we calculate utilization and weights for EACH node independently.
            This allows the algorithm to detect local bottlenecks (e.g., Trap Nodes with CPU=100%, RAM=0%)
            and assign high penalty weights specific to that node.
            """
            # 1. Calculate utilization for EACH node independently
            # Formula: U_j,k = (C_total - C_free) / C_total
            # Shape: (num_nodes, 3)
            # Added epsilon to N_cap to prevent division by zero (though capacity is static > 0)
            U_node = 1.0 - (N_free / (N_cap + 1e-9))
        
            # 2. Softmax per node (axis=1 is the resource dimension: CPU, RAM, BW)
            # Higher utilization -> Higher weight -> Higher penalty in fitness score
            # w_j,k = exp(eta * U_j,k) / sum(exp(eta * U_j,k'))
            exp_vals = np.exp(eta * U_node)
            w_k = exp_vals / np.sum(exp_vals, axis=1, keepdims=True)
        
            return w_k

    def run_full_micro_gwo(self, S, num_nodes, pop_size, max_iter):
        num_tasks = S.shape[0]
        wolves = np.random.rand(pop_size, num_tasks) * num_nodes
        wolves[0] = np.argmax(S, axis=1).astype(float) # Elite
        
        alpha_pos, beta_pos, delta_pos = wolves[0].copy(), wolves[0].copy(), wolves[0].copy()
        alpha_score, beta_score, delta_score = -1e9, -1e9, -1e9

        for t in range(max_iter):
            fitness = np.zeros(pop_size)
            idx = np.clip(np.round(wolves), 0, num_nodes-1).astype(int)
            for i in range(pop_size):
                scores = S[np.arange(num_tasks), idx[i]]
                if np.any(scores < -9000): fitness[i] = -1e9 
                else: fitness[i] = np.sum(scores)
            
            sorted_indices = np.argsort(fitness)[::-1]
            if fitness[sorted_indices[0]] > alpha_score:
                alpha_score, alpha_pos = fitness[sorted_indices[0]], wolves[sorted_indices[0]].copy()
            if pop_size > 1 and fitness[sorted_indices[1]] > beta_score:
                beta_score, beta_pos = fitness[sorted_indices[1]], wolves[sorted_indices[1]].copy()
            if pop_size > 2 and fitness[sorted_indices[2]] > delta_score:
                delta_score, delta_pos = fitness[sorted_indices[2]], wolves[sorted_indices[2]].copy()

            a = 2 - t * (2.0 / max_iter)
            for i in range(pop_size):
                r1, r2 = np.random.rand(3, num_tasks), np.random.rand(3, num_tasks)
                A = 2 * a * r1 - a 
                C = 2 * r2         
                D_alpha = np.abs(C[0] * alpha_pos - wolves[i])
                D_beta  = np.abs(C[1] * beta_pos  - wolves[i])
                D_delta = np.abs(C[2] * delta_pos - wolves[i])
                X1 = alpha_pos - A[0] * D_alpha
                X2 = beta_pos  - A[1] * D_beta
                X3 = delta_pos - A[2] * D_delta
                wolves[i] = (X1 + X2 + X3) / 3.0

        return np.clip(np.round(alpha_pos), 0, num_nodes-1).astype(int)

    def run_tps_ras(self, eta=5.0, pop_size=8, max_iter=5):
        start = time.time()
        self.reset_nodes()
        alloc = {}
        
        # Prepare Matrices
        T_mat = np.array([t.req for t in self.tasks])
        N_cap = np.array([n.capacity for n in self.nodes])
        N_free = np.array([n.current_free for n in self.nodes])
        Tiers = np.array([n.tier for n in self.nodes])
        
        pending = np.ones(len(self.tasks), dtype=bool)
        step_count = 0
        
        while np.any(pending) and step_count < len(self.tasks)*2:
            step_count += 1
            batch_size = 20
            batch_indices = np.where(pending)[0][:batch_size]
            curr_T = T_mat[batch_indices]
            
            # [Fix 1] Calculate Node-Level Dynamic Weights
            # omega shape: (num_nodes, 3)
            omega = self.compute_scarcity_weight_dynamic(N_cap, N_free, Tiers, eta)

            # [Fix 2] Calculate Weighted Resource Margin
            # fits shape: (Batch, Nodes)
            fits = np.all(curr_T[:, None, :] <= N_free[None, :, :], axis=2)
            
            # M_matrix shape: (Batch, Nodes, 3)
            M_matrix = 1 - (curr_T[:, None, :] / (N_free[None, :, :] + 1.0))
            
            # Weighted Sum with Node-Specific Weights
            # We broadcast omega from (Nodes, 3) to (1, Nodes, 3) to match M_matrix
            # Element-wise multiply: Margin[i,j,k] * Weight[j,k]
            # Sum over k (axis 2) -> Score[i,j]
            U_res = np.sum(M_matrix * omega[None, :, :], axis=2)
            
            # Penalties and Constraints
            score_net = np.zeros_like(Tiers, dtype=float)
            score_net[Tiers == 1] = 0.1
            score_net[Tiers == 2] = 0.5
            score_net[Tiers == 3] = 10.0
            
            S = U_res - score_net[None, :]
            S[~fits] = -9999
            
            # Micro-GWO Optimization
            if len(batch_indices) > 1:
                best_nodes = self.run_full_micro_gwo(S, len(self.nodes), pop_size, max_iter)
            else:
                best_nodes = np.argmax(S, axis=1)
            
            # Atomic Allocation
            allocated_in_batch = False
            for i, t_idx in enumerate(batch_indices):
                n_idx = best_nodes[i]
                if S[i, n_idx] > -9000 and np.all(T_mat[t_idx] <= N_free[n_idx]):
                    N_free[n_idx] -= T_mat[t_idx]
                    self.nodes[n_idx].current_free = N_free[n_idx]
                    alloc[self.tasks[t_idx].id] = n_idx
                    pending[t_idx] = False
                    allocated_in_batch = True
            
            # Cloud Fallback
            if not allocated_in_batch:
                cloud_indices = np.where(Tiers == 3)[0]
                if len(cloud_indices) > 0:
                    c_idx = cloud_indices[0]
                    t_idx = batch_indices[0]
                    if np.all(T_mat[t_idx] <= N_free[c_idx]):
                        N_free[c_idx] -= T_mat[t_idx]
                        self.nodes[c_idx].current_free = N_free[c_idx]
                        alloc[self.tasks[t_idx].id] = c_idx
                        pending[t_idx] = False
        
        return self.evaluate(alloc, time.time() - start)

# ==========================================
# 4. Main Execution
# ==========================================

if __name__ == "__main__":
    print(f"=== FULL STATISTICAL EVALUATION (30 Runs per Scenario) ===")
    all_scenario_data = []
    workloads = [100, 300, 500]
    repeat_times = 30
    
    for num_tasks in workloads:
        print(f"\n>>> Running Scenario: N={num_tasks} Tasks <<<")
        
        # Data Containers
        data = {
            'Edgewards': {'Util': [], 'NetCost': [], 'Time': []},
            'KubeDSM':       {'Util': [], 'NetCost': [], 'Time': []},
            'DRL':       {'Util': [], 'NetCost': [], 'Time': []},
            'TPS-RAS':   {'Util': [], 'NetCost': [], 'Time': []}
        }
        
        for i in range(repeat_times):
            if i % 10 == 0: print(f"  Repetition {i}/{repeat_times}")
            
            nodes, tasks = generate_heterogeneous_env(num_tasks)
            sim = Simulator(nodes, tasks)
            
            # 1. Edgewards
            r1 = sim.run_edgewards()
            # 2. KubeDSM
            r2 = sim.run_kubedsm()
            # 3. DRL-A2C
            # DRL needs stability in random seeds for reproducibility if desired,
            # but for statistical analysis, randomizing initialization is actually better.
            # Here we let it run naturally.
            r3 = sim.run_drl() 
            # 4. TPS-RAS
            r4 = sim.run_tps_ras(eta=5.0, pop_size=8, max_iter=5)
            
            for name, r in zip(['Edgewards', 'KubeDSM', 'DRL', 'TPS-RAS'], [r1, r2, r3, r4]):
                data[name]['Util'].append(r['Util%'])
                data[name]['NetCost'].append(r['NetCost'])
                data[name]['Time'].append(r['Time(s)'])

        print(f"\n[Results for N={num_tasks} | Mean ± Std | p-value vs TPS-RAS]")
        print(f"{'Algo':<12} | {'Util%':<20} | {'NetCost':<20} | {'Time(s)':<20} | {'p-val (Util)':<12}")
        print("-" * 95)
        
        tps_util = data['TPS-RAS']['Util']
        
        for name in ['Edgewards', 'KubeDSM', 'DRL', 'TPS-RAS']:
            u_mean, u_std = np.mean(data[name]['Util']), np.std(data[name]['Util'])
            n_mean, n_std = np.mean(data[name]['NetCost']), np.std(data[name]['NetCost'])
            t_mean, t_std = np.mean(data[name]['Time']), np.std(data[name]['Time'])
            
            if name == 'TPS-RAS':
                p_str = "-"
            else:
                # Two-sided t-test
                # Handle cases with zero variance (like K8s always being same)
                try:
                    t_stat, p_val = stats.ttest_ind(tps_util, data[name]['Util'], equal_var=False)
                    if p_val < 0.001: p_str = "<0.001"
                    else: p_str = f"{p_val:.3f}"
                except:
                    p_str = "N/A"
            
            print(f"{name:<12} | {u_mean:.2f} ± {u_std:.2f}      | {n_mean:.2f} ± {n_std:.2f}      | {t_mean:.4f} ± {t_std:.4f}   | {p_str:<12}")
            all_scenario_data.append(copy.deepcopy(data))
            
plot_box_results(all_scenario_data, workloads)