import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, quad
from scipy.optimize import minimize
import warnings
warnings.filterwarnings('ignore')

# Set random seed for reproducibility
np.random.seed(42)

class PopulationOptimizer:
    """
    A class to optimize population harvesting strategies under carrying capacity constraints.
    """
    
    def __init__(self, r=0.1, K=1000, N_min=100, T=50):
        """
        Initialize the population optimizer.
        
        Parameters:
        r (float): Intrinsic growth rate
        K (float): Carrying capacity
        N_min (float): Minimum viable population
        T (float): Time horizon
        """
        self.r = r
        self.K = K
        self.N_min = N_min
        self.T = T
        self.time_points = np.linspace(0, T, 1000)
    
    def logistic_growth_with_harvest(self, t, N, harvest_func):
        """
        Differential equation for logistic growth with harvesting.
        
        Parameters:
        t (float): Time
        N (array): Population size
        harvest_func (function): Harvesting function h(t)
        
        Returns:
        float: Rate of change of population
        """
        if N[0] < 0:
            return [0]  # Prevent negative population
        
        h_t = harvest_func(t)
        # Ensure harvesting doesn't exceed current population
        h_t = min(h_t, N[0] * 0.9)  # Leave at least 10% of population
        
        dN_dt = self.r * N[0] * (1 - N[0] / self.K) - h_t
        return [dN_dt]
    
    def simulate_population(self, N0, harvest_func):
        """
        Simulate population dynamics with given harvesting strategy.
        
        Parameters:
        N0 (float): Initial population
        harvest_func (function): Harvesting function
        
        Returns:
        tuple: (time_points, population_trajectory, harvest_trajectory)
        """
        def harvest_wrapper(t):
            return max(0, harvest_func(t))  # Ensure non-negative harvest
        
        # Solve the differential equation
        sol = solve_ivp(
            lambda t, N: self.logistic_growth_with_harvest(t, N, harvest_wrapper),
            [0, self.T],
            [N0],
            t_eval=self.time_points,
            method='RK45',
            rtol=1e-8
        )
        
        # Calculate harvest trajectory
        harvest_trajectory = np.array([harvest_wrapper(t) for t in self.time_points])
        
        return self.time_points, sol.y[0], harvest_trajectory
    
    def constant_harvest_strategy(self, h_constant):
        """Constant harvesting strategy."""
        return lambda t: h_constant
    
    def proportional_harvest_strategy(self, alpha):
        """Proportional harvesting strategy: h(t) = α * N(t)."""
        def harvest_func(t):
            # This requires knowing N(t), so we'll approximate
            # In practice, this would be implemented with feedback control
            N_approx = self.K / 2  # Rough approximation
            return alpha * N_approx
        return harvest_func
    
    def optimal_constant_harvest(self, N0):
        """
        Find optimal constant harvesting rate using the analytical solution.
        
        For constant harvesting, the equilibrium condition is:
        h* = r * N* * (1 - N*/K)
        
        Where N* is the equilibrium population that maximizes harvest.
        """
        # Maximum sustainable yield occurs at N* = K/2
        N_star = self.K / 2
        h_optimal = self.r * N_star * (1 - N_star / self.K)
        return h_optimal
    
    def calculate_total_harvest(self, harvest_trajectory, dt):
        """Calculate total harvest over the time period."""
        return np.trapz(harvest_trajectory, dx=dt)
    
    def plot_comparison(self, N0=800):
        """
        Plot comparison of different harvesting strategies.
        
        Parameters:
        N0 (float): Initial population
        """
        # Define different strategies
        strategies = {
            'No Harvest': self.constant_harvest_strategy(0),
            'Low Constant Harvest': self.constant_harvest_strategy(5),
            'Moderate Constant Harvest': self.constant_harvest_strategy(15),
            'High Constant Harvest': self.constant_harvest_strategy(30),
            'Optimal Constant Harvest': self.constant_harvest_strategy(
                self.optimal_constant_harvest(N0)
            )
        }
        
        # Set up the plots
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle('Population Growth Optimization: Harvesting Strategies Comparison', 
                    fontsize=16, fontweight='bold')
        
        colors = ['blue', 'green', 'orange', 'red', 'purple']
        total_harvests = {}
        
        # Simulate each strategy
        for i, (name, strategy) in enumerate(strategies.items()):
            time, population, harvest = self.simulate_population(N0, strategy)
            dt = time[1] - time[0]
            total_harvest = self.calculate_total_harvest(harvest, dt)
            total_harvests[name] = total_harvest
            
            # Plot population trajectory
            ax1.plot(time, population, label=name, color=colors[i], linewidth=2)
            
            # Plot harvest trajectory
            ax2.plot(time, harvest, label=name, color=colors[i], linewidth=2)
        
        # Configure population plot
        ax1.set_xlabel('Time')
        ax1.set_ylabel('Population Size')
        ax1.set_title('Population Trajectories')
        ax1.axhline(y=self.K, color='black', linestyle='--', alpha=0.7, 
                   label=f'Carrying Capacity (K={self.K})')
        ax1.axhline(y=self.N_min, color='red', linestyle='--', alpha=0.7, 
                   label=f'Minimum Viable Pop. ({self.N_min})')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # Configure harvest plot
        ax2.set_xlabel('Time')
        ax2.set_ylabel('Harvest Rate')
        ax2.set_title('Harvest Trajectories')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        # Plot total harvest comparison
        strategy_names = list(total_harvests.keys())
        harvest_values = list(total_harvests.values())
        bars = ax3.bar(range(len(strategy_names)), harvest_values, 
                      color=colors[:len(strategy_names)], alpha=0.7)
        ax3.set_xlabel('Harvesting Strategy')
        ax3.set_ylabel('Total Harvest')
        ax3.set_title('Total Harvest Comparison')
        ax3.set_xticks(range(len(strategy_names)))
        ax3.set_xticklabels(strategy_names, rotation=45, ha='right')
        ax3.grid(True, alpha=0.3)
        
        # Add value labels on bars
        for i, bar in enumerate(bars):
            height = bar.get_height()
            ax3.text(bar.get_x() + bar.get_width()/2., height,
                    f'{height:.1f}',
                    ha='center', va='bottom', fontweight='bold')
        
        # Plot phase portrait (Population vs dN/dt)
        N_range = np.linspace(self.N_min, self.K * 1.2, 100)
        
        # Without harvesting
        dN_dt_no_harvest = self.r * N_range * (1 - N_range / self.K)
        ax4.plot(N_range, dN_dt_no_harvest, 'b-', linewidth=2, 
                label='No Harvest')
        
        # With optimal constant harvest
        h_opt = self.optimal_constant_harvest(N0)
        dN_dt_with_harvest = self.r * N_range * (1 - N_range / self.K) - h_opt
        ax4.plot(N_range, dN_dt_with_harvest, 'purple', linewidth=2, 
                label=f'Optimal Harvest (h={h_opt:.2f})')
        
        ax4.axhline(y=0, color='black', linestyle='-', alpha=0.5)
        ax4.axvline(x=self.K/2, color='red', linestyle='--', alpha=0.7, 
                   label='MSY Population (K/2)')
        ax4.set_xlabel('Population Size (N)')
        ax4.set_ylabel('Population Growth Rate (dN/dt)')
        ax4.set_title('Phase Portrait: Growth Rate vs Population')
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        # Print numerical results
        print("\n" + "="*60)
        print("POPULATION GROWTH OPTIMIZATION RESULTS")
        print("="*60)
        print(f"Model Parameters:")
        print(f"  • Intrinsic growth rate (r): {self.r}")
        print(f"  • Carrying capacity (K): {self.K}")
        print(f"  • Initial population (N₀): {N0}")
        print(f"  • Time horizon (T): {self.T} years")
        print(f"  • Minimum viable population: {self.N_min}")
        
        print(f"\nOptimal Strategy Analysis:")
        optimal_h = self.optimal_constant_harvest(N0)
        print(f"  • Optimal constant harvest rate: {optimal_h:.2f} individuals/year")
        print(f"  • Equilibrium population at optimal harvest: {self.K/2:.0f} individuals")
        print(f"  • Maximum sustainable yield: {optimal_h:.2f} individuals/year")
        
        print(f"\nTotal Harvest Comparison (over {self.T} years):")
        for name, harvest in sorted(total_harvests.items(), 
                                  key=lambda x: x[1], reverse=True):
            print(f"  • {name}: {harvest:.1f} individuals")
        
        return total_harvests

# Advanced Analysis: Time-varying Optimal Control
class AdvancedPopulationOptimizer(PopulationOptimizer):
    """
    Advanced optimizer using optimal control theory for time-varying harvest rates.
    """
    
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        
    def bang_bang_harvest(self, N0, switch_time=None):
        """
        Implement bang-bang control strategy.
        
        Parameters:
        N0 (float): Initial population
        switch_time (float): Time to switch from max to min harvest
        """
        if switch_time is None:
            switch_time = self.T / 2
            
        def harvest_func(t):
            if t < switch_time:
                return 25  # High harvest initially
            else:
                return 5   # Low harvest later
        
        return harvest_func
    
    def adaptive_harvest(self, N0):
        """
        Adaptive harvest strategy that adjusts based on population level.
        """
        def harvest_func(t):
            # This is a simplified adaptive strategy
            # In practice, would need real-time population monitoring
            base_harvest = 10
            time_factor = 1 + 0.5 * np.sin(2 * np.pi * t / (self.T / 4))
            return base_harvest * time_factor
        
        return harvest_func
    
    def demonstrate_advanced_strategies(self, N0=800):
        """Demonstrate advanced harvesting strategies."""
        
        print("\n" + "="*60)
        print("ADVANCED POPULATION OPTIMIZATION STRATEGIES")
        print("="*60)
        
        # Advanced strategies
        strategies = {
            'Bang-Bang Control': self.bang_bang_harvest(N0),
            'Adaptive Harvest': self.adaptive_harvest(N0),
            'Optimal Constant': self.constant_harvest_strategy(
                self.optimal_constant_harvest(N0)
            )
        }
        
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle('Advanced Population Optimization Strategies', 
                    fontsize=16, fontweight='bold')
        
        colors = ['red', 'green', 'blue']
        results = {}
        
        for i, (name, strategy) in enumerate(strategies.items()):
            time, population, harvest = self.simulate_population(N0, strategy)
            dt = time[1] - time[0]
            total_harvest = self.calculate_total_harvest(harvest, dt)
            results[name] = {
                'time': time,
                'population': population,
                'harvest': harvest,
                'total_harvest': total_harvest
            }
            
            # Plot trajectories
            ax1.plot(time, population, label=name, color=colors[i], linewidth=2)
            ax2.plot(time, harvest, label=name, color=colors[i], linewidth=2)
        
        # Configure plots
        ax1.set_xlabel('Time')
        ax1.set_ylabel('Population Size')
        ax1.set_title('Population Trajectories - Advanced Strategies')
        ax1.axhline(y=self.K, color='black', linestyle='--', alpha=0.7)
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        ax2.set_xlabel('Time')
        ax2.set_ylabel('Harvest Rate')
        ax2.set_title('Harvest Trajectories - Advanced Strategies')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        # Cumulative harvest over time
        for i, (name, result) in enumerate(results.items()):
            cumulative_harvest = np.cumsum(result['harvest']) * (time[1] - time[0])
            ax3.plot(time, cumulative_harvest, label=name, color=colors[i], linewidth=2)
        
        ax3.set_xlabel('Time')
        ax3.set_ylabel('Cumulative Harvest')
        ax3.set_title('Cumulative Harvest Over Time')
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # Efficiency analysis (harvest per unit population)
        for i, (name, result) in enumerate(results.items()):
            efficiency = result['harvest'] / (result['population'] + 1e-10)  # Avoid division by zero
            ax4.plot(time, efficiency, label=name, color=colors[i], linewidth=2)
        
        ax4.set_xlabel('Time')
        ax4.set_ylabel('Harvest Efficiency (harvest/population)')
        ax4.set_title('Harvesting Efficiency Over Time')
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        # Summary statistics
        print(f"\nAdvanced Strategy Performance:")
        for name, result in results.items():
            avg_pop = np.mean(result['population'])
            min_pop = np.min(result['population'])
            max_harvest_rate = np.max(result['harvest'])
            print(f"  • {name}:")
            print(f"    - Total harvest: {result['total_harvest']:.1f}")
            print(f"    - Average population: {avg_pop:.1f}")
            print(f"    - Minimum population: {min_pop:.1f}")
            print(f"    - Maximum harvest rate: {max_harvest_rate:.1f}")
        
        return results

# Create and run the analysis
print("POPULATION GROWTH OPTIMIZATION ANALYSIS")
print("="*50)
print("Analyzing optimal harvesting strategies under carrying capacity constraints...")

# Basic optimization
optimizer = PopulationOptimizer(r=0.08, K=1000, N_min=50, T=50)
basic_results = optimizer.plot_comparison(N0=800)

# Advanced optimization
advanced_optimizer = AdvancedPopulationOptimizer(r=0.08, K=1000, N_min=50, T=50)
advanced_results = advanced_optimizer.demonstrate_advanced_strategies(N0=800)

print("\nAnalysis complete!")
print("\nKey Insights:")
print("1. Optimal constant harvest rate maximizes total yield")
print("2. Harvesting above optimal rate leads to population collapse")  
print("3. Advanced strategies can outperform constant harvest in specific scenarios")
print("4. Population must be maintained above minimum viable level")
print("5. Maximum Sustainable Yield occurs at K/2 population level")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import seaborn as sns
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

# Set up the cybersecurity game theory model
class CyberSecurityGame:
    def __init__(self, num_nodes=5, node_values=None, total_defense_budget=100, total_attack_budget=80):
        """
        Initialize the cybersecurity game
        
        Parameters:
        - num_nodes: Number of network nodes to defend/attack
        - node_values: Importance/value of each node (if None, random values assigned)
        - total_defense_budget: Total resources available for defense
        - total_attack_budget: Total resources available for attack
        """
        self.num_nodes = num_nodes
        self.total_defense_budget = total_defense_budget
        self.total_attack_budget = total_attack_budget
        
        # Assign node values (importance of each network component)
        if node_values is None:
            # Different types of network components with varying importance
            self.node_values = np.array([100, 80, 60, 40, 20])  # Critical to low priority
            self.node_names = ['Database Server', 'Web Server', 'Email Server', 'File Server', 'Test Server']
        else:
            self.node_values = np.array(node_values)
            self.node_names = [f'Node {i+1}' for i in range(num_nodes)]
    
    def attack_success_probability(self, attack_allocation, defense_allocation, epsilon=1):
        """
        Calculate probability of successful attack for each node
        
        P_attack(i) = A_i / (A_i + D_i + ε)
        """
        return attack_allocation / (attack_allocation + defense_allocation + epsilon)
    
    def attacker_utility(self, attack_allocation, defense_allocation):
        """
        Calculate total expected utility for attacker
        
        U_attacker = Σ V_i * P_attack(i)
        """
        success_probs = self.attack_success_probability(attack_allocation, defense_allocation)
        return np.sum(self.node_values * success_probs)
    
    def defender_utility(self, attack_allocation, defense_allocation):
        """
        Calculate utility for defender (negative of attacker utility)
        """
        return -self.attacker_utility(attack_allocation, defense_allocation)
    
    def optimize_defense_strategy(self, attack_allocation):
        """
        Find optimal defense allocation given attacker's strategy
        """
        def objective(defense_allocation):
            return self.attacker_utility(attack_allocation, defense_allocation)
        
        # Constraint: total defense budget
        constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - self.total_defense_budget}
        bounds = [(0, self.total_defense_budget) for _ in range(self.num_nodes)]
        
        # Initial guess: uniform allocation
        x0 = np.full(self.num_nodes, self.total_defense_budget / self.num_nodes)
        
        result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)
        return result.x if result.success else x0
    
    def optimize_attack_strategy(self, defense_allocation):
        """
        Find optimal attack allocation given defender's strategy
        """
        def objective(attack_allocation):
            return -self.attacker_utility(attack_allocation, defense_allocation)
        
        # Constraint: total attack budget
        constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - self.total_attack_budget}
        bounds = [(0, self.total_attack_budget) for _ in range(self.num_nodes)]
        
        # Initial guess: uniform allocation
        x0 = np.full(self.num_nodes, self.total_attack_budget / self.num_nodes)
        
        result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)
        return result.x if result.success else x0
    
    def find_nash_equilibrium(self, max_iterations=100, tolerance=1e-6):
        """
        Find Nash equilibrium using iterative best response
        """
        # Initialize with uniform allocations
        defense_allocation = np.full(self.num_nodes, self.total_defense_budget / self.num_nodes)
        attack_allocation = np.full(self.num_nodes, self.total_attack_budget / self.num_nodes)
        
        history = {
            'defense': [defense_allocation.copy()],
            'attack': [attack_allocation.copy()],
            'utilities': []
        }
        
        for iteration in range(max_iterations):
            # Optimize attack given current defense
            new_attack = self.optimize_attack_strategy(defense_allocation)
            
            # Optimize defense given new attack
            new_defense = self.optimize_defense_strategy(new_attack)
            
            # Calculate utilities
            attacker_util = self.attacker_utility(new_attack, new_defense)
            defender_util = self.defender_utility(new_attack, new_defense)
            
            history['attack'].append(new_attack.copy())
            history['defense'].append(new_defense.copy())
            history['utilities'].append((attacker_util, defender_util))
            
            # Check for convergence
            if (np.linalg.norm(new_attack - attack_allocation) < tolerance and 
                np.linalg.norm(new_defense - defense_allocation) < tolerance):
                print(f"Converged after {iteration + 1} iterations")
                break
            
            attack_allocation = new_attack
            defense_allocation = new_defense
        
        return {
            'optimal_defense': defense_allocation,
            'optimal_attack': attack_allocation,
            'attacker_utility': attacker_util,
            'defender_utility': defender_util,
            'history': history
        }

# Initialize and solve the game
game = CyberSecurityGame()
print("=== Cybersecurity Resource Allocation Game ===")
print(f"Number of nodes: {game.num_nodes}")
print(f"Node values: {game.node_values}")
print(f"Node names: {game.node_names}")
print(f"Defense budget: {game.total_defense_budget}")
print(f"Attack budget: {game.total_attack_budget}")

# Find Nash equilibrium
print("\nFinding Nash equilibrium...")
equilibrium = game.find_nash_equilibrium()

print("\n=== RESULTS ===")
print("Optimal Defense Allocation:")
for i, (name, allocation) in enumerate(zip(game.node_names, equilibrium['optimal_defense'])):
    print(f"  {name}: {allocation:.2f} resources ({allocation/game.total_defense_budget*100:.1f}%)")

print("\nOptimal Attack Allocation:")
for i, (name, allocation) in enumerate(zip(game.node_names, equilibrium['optimal_attack'])):
    print(f"  {name}: {allocation:.2f} resources ({allocation/game.total_attack_budget*100:.1f}%)")

print(f"\nAttacker Expected Utility: {equilibrium['attacker_utility']:.2f}")
print(f"Defender Utility: {equilibrium['defender_utility']:.2f}")

# Calculate success probabilities
success_probs = game.attack_success_probability(
    equilibrium['optimal_attack'], 
    equilibrium['optimal_defense']
)

print("\nAttack Success Probabilities:")
for i, (name, prob) in enumerate(zip(game.node_names, success_probs)):
    print(f"  {name}: {prob:.3f} ({prob*100:.1f}%)")

# Create comprehensive visualizations
fig = plt.figure(figsize=(16, 12))

# 1. Resource Allocation Comparison
ax1 = plt.subplot(2, 3, 1)
x = np.arange(len(game.node_names))
width = 0.35

bars1 = ax1.bar(x - width/2, equilibrium['optimal_defense'], width, 
               label='Defense Allocation', color='#2E86AB', alpha=0.8)
bars2 = ax1.bar(x + width/2, equilibrium['optimal_attack'], width,
               label='Attack Allocation', color='#A23B72', alpha=0.8)

ax1.set_xlabel('Network Nodes')
ax1.set_ylabel('Resource Allocation')
ax1.set_title('Optimal Resource Allocation\nat Nash Equilibrium')
ax1.set_xticks(x)
ax1.set_xticklabels([name.split()[0] for name in game.node_names], rotation=45)
ax1.legend()
ax1.grid(True, alpha=0.3)

# Add value labels on bars
for bar in bars1:
    height = bar.get_height()
    ax1.annotate(f'{height:.1f}',
                xy=(bar.get_x() + bar.get_width() / 2, height),
                xytext=(0, 3),
                textcoords="offset points",
                ha='center', va='bottom', fontsize=8)

for bar in bars2:
    height = bar.get_height()
    ax1.annotate(f'{height:.1f}',
                xy=(bar.get_x() + bar.get_width() / 2, height),
                xytext=(0, 3),
                textcoords="offset points",
                ha='center', va='bottom', fontsize=8)

# 2. Attack Success Probabilities
ax2 = plt.subplot(2, 3, 2)
colors = plt.cm.Reds(success_probs / max(success_probs))
bars = ax2.bar(range(len(game.node_names)), success_probs, color=colors, alpha=0.8)
ax2.set_xlabel('Network Nodes')
ax2.set_ylabel('Success Probability')
ax2.set_title('Attack Success Probabilities')
ax2.set_xticks(range(len(game.node_names)))
ax2.set_xticklabels([name.split()[0] for name in game.node_names], rotation=45)
ax2.grid(True, alpha=0.3)

# Add percentage labels
for i, (bar, prob) in enumerate(zip(bars, success_probs)):
    ax2.annotate(f'{prob*100:.1f}%',
                xy=(bar.get_x() + bar.get_width() / 2, bar.get_height()),
                xytext=(0, 3),
                textcoords="offset points",
                ha='center', va='bottom', fontsize=9)

# 3. Node Value vs Resource Allocation
ax3 = plt.subplot(2, 3, 3)
ax3.scatter(game.node_values, equilibrium['optimal_defense'], 
           s=100, alpha=0.7, color='#2E86AB', label='Defense')
ax3.scatter(game.node_values, equilibrium['optimal_attack'], 
           s=100, alpha=0.7, color='#A23B72', label='Attack')

for i, name in enumerate(game.node_names):
    ax3.annotate(name.split()[0], 
                xy=(game.node_values[i], equilibrium['optimal_defense'][i]),
                xytext=(5, 5), textcoords='offset points', fontsize=8)

ax3.set_xlabel('Node Value (Importance)')
ax3.set_ylabel('Resource Allocation')
ax3.set_title('Value vs Resource Allocation')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Convergence Analysis
if len(equilibrium['history']['utilities']) > 1:
    ax4 = plt.subplot(2, 3, 4)
    iterations = range(1, len(equilibrium['history']['utilities']) + 1)
    attacker_utils = [u[0] for u in equilibrium['history']['utilities']]
    defender_utils = [u[1] for u in equilibrium['history']['utilities']]
    
    ax4.plot(iterations, attacker_utils, 'o-', color='#A23B72', label='Attacker Utility', linewidth=2)
    ax4.plot(iterations, defender_utils, 's-', color='#2E86AB', label='Defender Utility', linewidth=2)
    ax4.set_xlabel('Iteration')
    ax4.set_ylabel('Utility')
    ax4.set_title('Convergence to Nash Equilibrium')
    ax4.legend()
    ax4.grid(True, alpha=0.3)

# 5. Risk-Return Analysis
ax5 = plt.subplot(2, 3, 5)
expected_damage = game.node_values * success_probs
total_expected_damage = np.sum(expected_damage)

wedges, texts, autotexts = ax5.pie(expected_damage, labels=[name.split()[0] for name in game.node_names],
                                  autopct='%1.1f%%', startangle=90, colors=plt.cm.Reds_r(np.linspace(0.3, 0.8, len(game.node_names))))
ax5.set_title(f'Expected Damage Distribution\nTotal: {total_expected_damage:.1f}')

# 6. Security Investment Efficiency
ax6 = plt.subplot(2, 3, 6)
defense_per_value = equilibrium['optimal_defense'] / game.node_values
attack_per_value = equilibrium['optimal_attack'] / game.node_values

x = np.arange(len(game.node_names))
bars1 = ax6.bar(x - width/2, defense_per_value, width, 
               label='Defense/Value Ratio', color='#2E86AB', alpha=0.8)
bars2 = ax6.bar(x + width/2, attack_per_value, width,
               label='Attack/Value Ratio', color='#A23B72', alpha=0.8)

ax6.set_xlabel('Network Nodes')
ax6.set_ylabel('Resource per Unit Value')
ax6.set_title('Resource Efficiency Analysis')
ax6.set_xticks(x)
ax6.set_xticklabels([name.split()[0] for name in game.node_names], rotation=45)
ax6.legend()
ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print additional analysis
print("\n=== STRATEGIC ANALYSIS ===")
print(f"Total Expected Damage: {np.sum(expected_damage):.2f}")
print(f"Most Vulnerable Node: {game.node_names[np.argmax(success_probs)]} ({success_probs[np.argmax(success_probs)]*100:.1f}% success rate)")
print(f"Best Protected Node: {game.node_names[np.argmin(success_probs)]} ({success_probs[np.argmin(success_probs)]*100:.1f}% success rate)")

defense_efficiency = equilibrium['optimal_defense'] / game.node_values
attack_focus = equilibrium['optimal_attack'] / game.total_attack_budget

print(f"\nDefense Efficiency (resources per unit value):")
for i, (name, eff) in enumerate(zip(game.node_names, defense_efficiency)):
    print(f"  {name}: {eff:.3f}")

print(f"\nAttack Focus Distribution:")
for i, (name, focus) in enumerate(zip(game.node_names, attack_focus)):
    print(f"  {name}: {focus*100:.1f}% of total attack budget")

def attack_success_probability(self, attack_allocation, defense_allocation, epsilon=1):
    return attack_allocation / (attack_allocation + defense_allocation + epsilon)

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')

# Set random seed for reproducibility
np.random.seed(42)

class ArcticResourceOptimizer:
    def __init__(self):
        """
        Initialize the Arctic Resource Optimization model
        """
        # Define countries and territories
        self.countries = ['USA', 'Canada', 'Russia', 'Norway', 'Denmark']
        self.territories = ['Sector_A', 'Sector_B', 'Sector_C', 'Sector_D', 'Sector_E', 'Sector_F']
        
        # Resource values in each territory (billions USD)
        self.resource_values = np.array([150, 200, 180, 120, 160, 140])
        
        # Extraction efficiency matrix (countries x territories)
        # Higher values indicate better extraction capabilities
        self.efficiency_matrix = np.array([
            [0.8, 0.7, 0.6, 0.9, 0.5, 0.7],  # USA
            [0.9, 0.8, 0.5, 0.7, 0.6, 0.8],  # Canada
            [0.7, 0.6, 0.9, 0.6, 0.8, 0.7],  # Russia
            [0.6, 0.7, 0.7, 0.8, 0.9, 0.6],  # Norway
            [0.5, 0.6, 0.6, 0.7, 0.7, 0.9]   # Denmark
        ])
        
        # Development costs for each territory (billions USD)
        self.development_costs = np.array([80, 100, 90, 60, 85, 75])
        
        # Country development budgets (billions USD)
        self.country_budgets = np.array([300, 250, 280, 200, 180])
        
        self.n_countries = len(self.countries)
        self.n_territories = len(self.territories)
    
    def setup_optimization_problem(self):
        """
        Set up the linear programming problem for resource allocation
        """
        # Decision variables: x_ij represents allocation of territory j to country i
        # We flatten the matrix into a vector for linear programming
        n_vars = self.n_countries * self.n_territories
        
        # Objective function coefficients (maximize benefits)
        # We need to negate for minimization in scipy.optimize.linprog
        c = []
        for i in range(self.n_countries):
            for j in range(self.n_territories):
                benefit = self.resource_values[j] * self.efficiency_matrix[i][j]
                c.append(-benefit)  # Negative for maximization
        
        c = np.array(c)
        
        # Constraint matrices
        A_ub = []  # Inequality constraints (<=)
        b_ub = []
        
        # Territory allocation constraints: sum of allocations for each territory <= 1
        for j in range(self.n_territories):
            constraint = np.zeros(n_vars)
            for i in range(self.n_countries):
                idx = i * self.n_territories + j
                constraint[idx] = 1
            A_ub.append(constraint)
            b_ub.append(1.0)
        
        # Budget constraints for each country
        for i in range(self.n_countries):
            constraint = np.zeros(n_vars)
            for j in range(self.n_territories):
                idx = i * self.n_territories + j
                constraint[idx] = self.development_costs[j]
            A_ub.append(constraint)
            b_ub.append(self.country_budgets[i])
        
        A_ub = np.array(A_ub)
        b_ub = np.array(b_ub)
        
        # Bounds: all variables >= 0, <= 1
        bounds = [(0, 1) for _ in range(n_vars)]
        
        return c, A_ub, b_ub, bounds
    
    def solve_optimization(self):
        """
        Solve the optimization problem
        """
        c, A_ub, b_ub, bounds = self.setup_optimization_problem()
        
        # Solve using linear programming
        result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')
        
        if result.success:
            # Reshape solution back to matrix form
            solution_matrix = result.x.reshape(self.n_countries, self.n_territories)
            total_benefit = -result.fun  # Convert back from minimization
            
            return solution_matrix, total_benefit, result
        else:
            print("Optimization failed!")
            return None, None, result
    
    def analyze_results(self, solution_matrix, total_benefit):
        """
        Analyze and interpret the optimization results
        """
        print("=== ARCTIC RESOURCE ALLOCATION OPTIMIZATION RESULTS ===\n")
        print(f"Total System Benefit: ${total_benefit:.2f} billion\n")
        
        # Create results DataFrame
        results_df = pd.DataFrame(
            solution_matrix,
            index=self.countries,
            columns=self.territories
        )
        
        print("Allocation Matrix (fraction of territory allocated to each country):")
        print(results_df.round(3))
        print()
        
        # Calculate benefits per country
        country_benefits = []
        for i in range(self.n_countries):
            benefit = 0
            for j in range(self.n_territories):
                benefit += (solution_matrix[i][j] * 
                           self.resource_values[j] * 
                           self.efficiency_matrix[i][j])
            country_benefits.append(benefit)
        
        benefits_df = pd.DataFrame({
            'Country': self.countries,
            'Total_Benefit': country_benefits,
            'Budget_Used': [np.sum(solution_matrix[i] * self.development_costs) 
                           for i in range(self.n_countries)]
        })
        
        print("Country Benefits and Budget Usage:")
        print(benefits_df.round(2))
        
        return results_df, benefits_df
    
    def create_visualizations(self, solution_matrix, benefits_df):
        """
        Create comprehensive visualizations of the results
        """
        fig = plt.figure(figsize=(20, 15))
        
        # 1. Territory Allocation Heatmap
        plt.subplot(2, 3, 1)
        sns.heatmap(solution_matrix, 
                   xticklabels=self.territories,
                   yticklabels=self.countries,
                   annot=True, 
                   fmt='.3f',
                   cmap='YlOrRd',
                   cbar_kws={'label': 'Allocation Fraction'})
        plt.title('Territory Allocation Matrix\n(Fraction of Each Territory by Country)', 
                 fontsize=12, fontweight='bold')
        plt.xlabel('Arctic Territories')
        plt.ylabel('Countries')
        
        # 2. Country Benefits Bar Chart
        plt.subplot(2, 3, 2)
        bars = plt.bar(benefits_df['Country'], benefits_df['Total_Benefit'], 
                      color=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'])
        plt.title('Total Benefits by Country', fontsize=12, fontweight='bold')
        plt.xlabel('Countries')
        plt.ylabel('Benefits (Billion USD)')
        plt.xticks(rotation=45)
        
        # Add value labels on bars
        for bar in bars:
            height = bar.get_height()
            plt.text(bar.get_x() + bar.get_width()/2., height + 2,
                    f'${height:.1f}B', ha='center', va='bottom')
        
        # 3. Budget Utilization
        plt.subplot(2, 3, 3)
        budget_utilization = benefits_df['Budget_Used'] / self.country_budgets * 100
        bars = plt.bar(benefits_df['Country'], budget_utilization,
                      color=['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4', '#FECA57'])
        plt.title('Budget Utilization by Country', fontsize=12, fontweight='bold')
        plt.xlabel('Countries')
        plt.ylabel('Budget Utilization (%)')
        plt.xticks(rotation=45)
        plt.ylim(0, 100)
        
        for bar in bars:
            height = bar.get_height()
            plt.text(bar.get_x() + bar.get_width()/2., height + 1,
                    f'{height:.1f}%', ha='center', va='bottom')
        
        # 4. Resource Value vs Territory Allocation
        plt.subplot(2, 3, 4)
        territory_allocation = np.sum(solution_matrix, axis=0)
        plt.scatter(self.resource_values, territory_allocation, 
                   s=100, alpha=0.7, c=range(len(self.territories)), cmap='viridis')
        plt.xlabel('Resource Value (Billion USD)')
        plt.ylabel('Total Territory Allocation')
        plt.title('Resource Value vs Territory Allocation', fontsize=12, fontweight='bold')
        
        # Add territory labels
        for i, territory in enumerate(self.territories):
            plt.annotate(territory, (self.resource_values[i], territory_allocation[i]),
                        xytext=(5, 5), textcoords='offset points', fontsize=8)
        
        # 5. Efficiency vs Allocation 3D Plot
        ax = fig.add_subplot(2, 3, 5, projection='3d')
        
        # Create meshgrid for surface plot
        x_country = np.arange(len(self.countries))
        y_territory = np.arange(len(self.territories))
        X, Y = np.meshgrid(x_country, y_territory)
        Z_efficiency = self.efficiency_matrix.T
        Z_allocation = solution_matrix.T
        
        # Plot efficiency surface
        surf1 = ax.plot_surface(X, Y, Z_efficiency, alpha=0.3, cmap='coolwarm', label='Efficiency')
        
        # Plot allocation points
        for i in range(len(self.countries)):
            for j in range(len(self.territories)):
                if solution_matrix[i][j] > 0.1:  # Only show significant allocations
                    ax.scatter(i, j, solution_matrix[i][j], 
                             s=solution_matrix[i][j]*500, 
                             c=self.efficiency_matrix[i][j], 
                             cmap='viridis', alpha=0.8)
        
        ax.set_xlabel('Countries')
        ax.set_ylabel('Territories')
        ax.set_zlabel('Value')
        ax.set_title('3D: Efficiency vs Allocation', fontsize=12, fontweight='bold')
        ax.set_xticks(range(len(self.countries)))
        ax.set_xticklabels(self.countries, rotation=45)
        ax.set_yticks(range(len(self.territories)))
        ax.set_yticklabels(self.territories, rotation=45)
        
        # 6. Pie Chart of Total Benefits
        plt.subplot(2, 3, 6)
        colors = ['#FF9999', '#66B2FF', '#99FF99', '#FFB366', '#FF99FF']
        plt.pie(benefits_df['Total_Benefit'], labels=benefits_df['Country'], 
                autopct='%1.1f%%', colors=colors, startangle=90)
        plt.title('Distribution of Total Benefits', fontsize=12, fontweight='bold')
        
        plt.tight_layout()
        plt.show()
        
        return fig

# Create and run the optimization model
print("Initializing Arctic Resource Optimization Model...")
optimizer = ArcticResourceOptimizer()

print("Setting up optimization problem...")
solution_matrix, total_benefit, result = optimizer.solve_optimization()

if solution_matrix is not None:
    print("Optimization completed successfully!")
    print("Analyzing results...\n")
    
    results_df, benefits_df = optimizer.analyze_results(solution_matrix, total_benefit)
    
    print("\n" + "="*60)
    print("CREATING VISUALIZATIONS...")
    print("="*60)
    
    fig = optimizer.create_visualizations(solution_matrix, benefits_df)
    
    # Additional analysis
    print("\n" + "="*60)
    print("ADDITIONAL STRATEGIC INSIGHTS")
    print("="*60)
    
    # Territory competition analysis
    print("\nTerritory Competition Analysis:")
    territory_competition = np.sum(solution_matrix > 0.01, axis=0)
    for i, territory in enumerate(optimizer.territories):
        competing_countries = territory_competition[i]
        total_allocation = np.sum(solution_matrix[:, i])
        print(f"{territory}: {competing_countries} countries competing, "
              f"{total_allocation:.1%} total allocation")
    
    # Efficiency analysis
    print("\nEfficiency Analysis:")
    for i, country in enumerate(optimizer.countries):
        avg_efficiency = np.mean(optimizer.efficiency_matrix[i])
        allocated_territories = np.sum(solution_matrix[i] > 0.01)
        print(f"{country}: Avg efficiency {avg_efficiency:.3f}, "
              f"Active in {allocated_territories} territories")
else:
    print("Optimization failed. Please check constraints and parameters.")

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
from scipy.optimize import minimize
import pandas as pd
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

# Set up the optimization problem for border management
class BorderManagementOptimizer:
    def __init__(self):
        # Define passenger types
        self.passenger_types = ['Regular', 'Business', 'Diplomatic', 'Transit']
        self.checkpoint_types = ['Standard', 'Expedited', 'Security-Enhanced', 'Automated']
        
        # Passenger demand (arrivals per hour)
        self.demand = np.array([150, 80, 20, 60])  # Regular, Business, Diplomatic, Transit
        
        # Processing times (minutes per passenger) - passenger_type x checkpoint_type
        self.processing_times = np.array([
            [3.5, 2.0, 5.0, 1.5],  # Regular passengers
            [2.5, 1.5, 4.0, 1.2],  # Business passengers  
            [1.5, 1.0, 2.5, 1.0],  # Diplomatic passengers
            [2.0, 1.8, 3.0, 1.1]   # Transit passengers
        ])
        
        # Processing costs ($ per passenger) - passenger_type x checkpoint_type
        self.costs = np.array([
            [5.0, 8.0, 12.0, 3.0],   # Regular passengers
            [5.5, 8.5, 12.5, 3.2],  # Business passengers
            [4.0, 6.0, 10.0, 2.5],  # Diplomatic passengers
            [4.5, 7.0, 11.0, 2.8]   # Transit passengers
        ])
        
        # Checkpoint capacities (passengers per hour)
        self.capacities = np.array([100, 60, 40, 120])  # Standard, Expedited, Security-Enhanced, Automated
        
        # Weight factors for different passenger types (priority weights)
        self.weights = np.array([1.0, 2.0, 3.0, 1.5])  # Higher weight = higher priority
        
    def solve_optimization(self):
        """
        Solve the border management optimization problem using linear programming
        """
        n_passengers = len(self.passenger_types)
        n_checkpoints = len(self.checkpoint_types)
        
        # Decision variables: x[i,j] = number of passenger type i assigned to checkpoint j
        # Flattened to 1D array for linprog
        
        # Objective function coefficients (minimize total cost + weighted processing time)
        c = []
        for i in range(n_passengers):
            for j in range(n_checkpoints):
                # Combined cost and time factor
                cost_factor = self.costs[i,j]
                time_factor = self.weights[i] * self.processing_times[i,j] * 0.1  # Scale time cost
                c.append(cost_factor + time_factor)
        
        # Inequality constraints: Ax <= b
        A_ub = []
        b_ub = []
        
        # Capacity constraints: sum over passenger types <= checkpoint capacity
        for j in range(n_checkpoints):
            constraint = [0] * (n_passengers * n_checkpoints)
            for i in range(n_passengers):
                constraint[i * n_checkpoints + j] = 1
            A_ub.append(constraint)
            b_ub.append(self.capacities[j])
        
        # Equality constraints: Ax = b (demand satisfaction)
        A_eq = []
        b_eq = []
        
        # Demand constraints: sum over checkpoints = passenger demand
        for i in range(n_passengers):
            constraint = [0] * (n_passengers * n_checkpoints)
            for j in range(n_checkpoints):
                constraint[i * n_checkpoints + j] = 1
            A_eq.append(constraint)
            b_eq.append(self.demand[i])
        
        # Bounds: all variables >= 0
        bounds = [(0, None) for _ in range(n_passengers * n_checkpoints)]
        
        # Solve the linear program
        result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, 
                        bounds=bounds, method='highs')
        
        if result.success:
            # Reshape solution back to matrix form
            solution = result.x.reshape(n_passengers, n_checkpoints)
            return solution, result.fun
        else:
            return None, None
    
    def analyze_solution(self, solution):
        """
        Analyze the optimization solution and calculate metrics
        """
        n_passengers, n_checkpoints = solution.shape
        
        # Calculate utilization rates
        utilization = np.sum(solution, axis=0) / self.capacities
        
        # Calculate average processing times by passenger type
        avg_times = []
        for i in range(n_passengers):
            total_time = np.sum(solution[i] * self.processing_times[i])
            total_passengers = np.sum(solution[i])
            if total_passengers > 0:
                avg_times.append(total_time / total_passengers)
            else:
                avg_times.append(0)
        
        # Calculate total costs by checkpoint
        checkpoint_costs = []
        for j in range(n_checkpoints):
            cost = np.sum(solution[:, j] * self.costs[:, j])
            checkpoint_costs.append(cost)
        
        return {
            'utilization': utilization,
            'avg_processing_times': avg_times,
            'checkpoint_costs': checkpoint_costs,
            'total_cost': np.sum(solution * self.costs),
            'total_passengers': np.sum(solution)
        }

# Create and solve the optimization problem
optimizer = BorderManagementOptimizer()
solution, objective_value = optimizer.solve_optimization()

print("=== BORDER MANAGEMENT OPTIMIZATION RESULTS ===\n")

if solution is not None:
    print(f"Optimization successful!")
    print(f"Total objective value: ${objective_value:.2f}\n")
    
    # Analyze the solution
    analysis = optimizer.analyze_solution(solution)
    
    print("PASSENGER ALLOCATION MATRIX:")
    print("Rows: Passenger Types, Columns: Checkpoint Types")
    df_solution = pd.DataFrame(solution, 
                              index=optimizer.passenger_types,
                              columns=optimizer.checkpoint_types)
    print(df_solution.round(1))
    print()
    
    print("CHECKPOINT UTILIZATION:")
    for i, checkpoint in enumerate(optimizer.checkpoint_types):
        util_pct = analysis['utilization'][i] * 100
        print(f"{checkpoint}: {util_pct:.1f}% ({analysis['utilization'][i] * optimizer.capacities[i]:.0f}/{optimizer.capacities[i]} passengers/hour)")
    print()
    
    print("AVERAGE PROCESSING TIMES BY PASSENGER TYPE:")
    for i, ptype in enumerate(optimizer.passenger_types):
        print(f"{ptype}: {analysis['avg_processing_times'][i]:.2f} minutes")
    print()
    
    print("COSTS BY CHECKPOINT:")
    for i, checkpoint in enumerate(optimizer.checkpoint_types):
        print(f"{checkpoint}: ${analysis['checkpoint_costs'][i]:.2f}/hour")
    
    print(f"\nTotal hourly processing cost: ${analysis['total_cost']:.2f}")
    print(f"Total passengers processed: {analysis['total_passengers']:.0f}/hour")
    print(f"Average cost per passenger: ${analysis['total_cost']/analysis['total_passengers']:.2f}")

else:
    print("Optimization failed - no feasible solution found")
    
# Create comprehensive visualizations
plt.style.use('default')
fig = plt.figure(figsize=(20, 15))

if solution is not None:
    # 1. Passenger allocation heatmap
    ax1 = plt.subplot(2, 3, 1)
    sns.heatmap(solution, annot=True, fmt='.1f', cmap='YlOrRd', 
                xticklabels=optimizer.checkpoint_types,
                yticklabels=optimizer.passenger_types,
                cbar_kws={'label': 'Passengers/hour'})
    plt.title('Optimal Passenger Allocation\n(Passengers per Hour)', fontsize=14, fontweight='bold')
    plt.xlabel('Checkpoint Types')
    plt.ylabel('Passenger Types')
    
    # 2. Checkpoint utilization
    ax2 = plt.subplot(2, 3, 2)
    utilization_pct = analysis['utilization'] * 100
    bars = plt.bar(optimizer.checkpoint_types, utilization_pct, 
                   color=['#2E86AB', '#A23B72', '#F18F01', '#C73E1D'])
    plt.axhline(y=100, color='red', linestyle='--', alpha=0.7, label='Maximum Capacity')
    plt.title('Checkpoint Utilization Rates', fontsize=14, fontweight='bold')
    plt.xlabel('Checkpoint Types')
    plt.ylabel('Utilization (%)')
    plt.xticks(rotation=45)
    plt.legend()
    
    # Add value labels on bars
    for i, bar in enumerate(bars):
        height = bar.get_height()
        plt.text(bar.get_x() + bar.get_width()/2., height + 1,
                f'{height:.1f}%', ha='center', va='bottom', fontweight='bold')
    
    # 3. Average processing times
    ax3 = plt.subplot(2, 3, 3)
    colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4']
    bars = plt.bar(optimizer.passenger_types, analysis['avg_processing_times'], color=colors)
    plt.title('Average Processing Time by Passenger Type', fontsize=14, fontweight='bold')
    plt.xlabel('Passenger Types')
    plt.ylabel('Processing Time (minutes)')
    plt.xticks(rotation=45)
    
    # Add value labels
    for i, bar in enumerate(bars):
        height = bar.get_height()
        plt.text(bar.get_x() + bar.get_width()/2., height + 0.05,
                f'{height:.2f}', ha='center', va='bottom', fontweight='bold')
    
    # 4. Cost breakdown by checkpoint
    ax4 = plt.subplot(2, 3, 4)
    explode = (0.05, 0.05, 0.05, 0.05)
    plt.pie(analysis['checkpoint_costs'], labels=optimizer.checkpoint_types, autopct='%1.1f%%',
            explode=explode, colors=['#FF9999', '#66B2FF', '#99FF99', '#FFCC99'])
    plt.title('Cost Distribution by Checkpoint Type\n(Total: ${:.0f}/hour)'.format(analysis['total_cost']), 
              fontsize=14, fontweight='bold')
    
    # 5. Processing time comparison matrix
    ax5 = plt.subplot(2, 3, 5)
    sns.heatmap(optimizer.processing_times, annot=True, fmt='.1f', cmap='RdYlBu_r',
                xticklabels=optimizer.checkpoint_types,
                yticklabels=optimizer.passenger_types,
                cbar_kws={'label': 'Processing Time (minutes)'})
    plt.title('Processing Time Matrix\n(Base Times)', fontsize=14, fontweight='bold')
    plt.xlabel('Checkpoint Types')
    plt.ylabel('Passenger Types')
    
    # 6. Demand vs Capacity analysis
    ax6 = plt.subplot(2, 3, 6)
    x_pos = np.arange(len(optimizer.checkpoint_types))
    width = 0.35
    
    allocated = np.sum(solution, axis=0)
    capacity = optimizer.capacities
    
    bars1 = plt.bar(x_pos - width/2, allocated, width, label='Allocated', color='#FF7F7F')
    bars2 = plt.bar(x_pos + width/2, capacity, width, label='Capacity', color='#7F7FFF')
    
    plt.title('Allocation vs Capacity by Checkpoint', fontsize=14, fontweight='bold')
    plt.xlabel('Checkpoint Types')
    plt.ylabel('Passengers/hour')
    plt.xticks(x_pos, optimizer.checkpoint_types, rotation=45)
    plt.legend()
    
    # Add value labels
    for i, (bar1, bar2) in enumerate(zip(bars1, bars2)):
        plt.text(bar1.get_x() + bar1.get_width()/2., bar1.get_height() + 2,
                f'{allocated[i]:.0f}', ha='center', va='bottom', fontsize=10)
        plt.text(bar2.get_x() + bar2.get_width()/2., bar2.get_height() + 2,
                f'{capacity[i]:.0f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.show()

# Performance metrics summary
print("\n" + "="*60)
print("OPTIMIZATION PERFORMANCE METRICS")
print("="*60)

if solution is not None:
    # Calculate efficiency metrics
    total_capacity = np.sum(optimizer.capacities)
    total_demand = np.sum(optimizer.demand)
    overall_utilization = np.sum(solution) / total_capacity
    
    print(f"Overall system utilization: {overall_utilization*100:.1f}%")
    print(f"Total system capacity: {total_capacity} passengers/hour")
    print(f"Total demand served: {np.sum(solution):.0f} passengers/hour")
    print(f"Demand coverage: {(np.sum(solution)/total_demand)*100:.1f}%")
    
    # Calculate weighted average processing time
    total_weighted_time = np.sum([analysis['avg_processing_times'][i] * optimizer.weights[i] * np.sum(solution[i])
                                 for i in range(len(optimizer.passenger_types))])
    weighted_avg_time = total_weighted_time / np.sum(solution)
    print(f"Weighted average processing time: {weighted_avg_time:.2f} minutes")
    
    # Calculate cost efficiency
    cost_per_minute_saved = analysis['total_cost'] / (60 * np.sum(solution))  # Cost per passenger-hour
    print(f"Cost efficiency: ${cost_per_minute_saved:.4f} per passenger-minute")