import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from itertools import combinations, chain
import seaborn as sns
from scipy.optimize import minimize
import networkx as nx
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

class AllianceOptimizer:
    """
    A class to optimize alliance relationships using game theory principles.
    """
    
    def __init__(self, nations, security_weights, cost_coefficients, power_exponents):
        """
        Initialize the alliance optimizer.
        
        Parameters:
        - nations: list of nation names
        - security_weights: weights for security value calculation
        - cost_coefficients: alpha values for cost calculation
        - power_exponents: beta values for cost calculation
        """
        self.nations = nations
        self.n = len(nations)
        self.security_weights = np.array(security_weights)
        self.cost_coefficients = np.array(cost_coefficients)
        self.power_exponents = np.array(power_exponents)
        
        # Create synergy matrix (how well nations work together)
        np.random.seed(42)  # For reproducible results
        self.synergy_matrix = np.random.uniform(0.8, 1.5, (self.n, self.n))
        np.fill_diagonal(self.synergy_matrix, 1.0)
        # Make matrix symmetric
        self.synergy_matrix = (self.synergy_matrix + self.synergy_matrix.T) / 2
    
    def security_value(self, coalition_indices):
        """
        Calculate the security value for a given coalition.
        
        V_i(S) = w_i * sum(synergy_ij for j in S) * sqrt(|S|)
        """
        if len(coalition_indices) == 0:
            return 0
        
        coalition_size = len(coalition_indices)
        total_value = 0
        
        for i in coalition_indices:
            # Security comes from synergy with coalition members
            synergy_sum = sum(self.synergy_matrix[i][j] for j in coalition_indices)
            individual_value = self.security_weights[i] * synergy_sum * np.sqrt(coalition_size)
            total_value += individual_value
            
        return total_value
    
    def cost_function(self, nation_idx, coalition_size):
        """
        Calculate cost for a nation in a coalition of given size.
        
        C_i(S) = alpha_i * |S|^beta_i
        """
        return self.cost_coefficients[nation_idx] * (coalition_size ** self.power_exponents[nation_idx])
    
    def total_cost(self, coalition_indices):
        """Calculate total cost for all nations in a coalition."""
        coalition_size = len(coalition_indices)
        return sum(self.cost_function(i, coalition_size) for i in coalition_indices)
    
    def utility_function(self, coalition_indices):
        """
        Calculate utility (security value - cost) for a coalition.
        """
        if len(coalition_indices) == 0:
            return 0
        
        security = self.security_value(coalition_indices)
        cost = self.total_cost(coalition_indices)
        return security - cost
    
    def generate_all_coalitions(self):
        """Generate all possible non-empty coalitions."""
        all_coalitions = []
        for r in range(1, self.n + 1):
            for coalition in combinations(range(self.n), r):
                all_coalitions.append(list(coalition))
        return all_coalitions
    
    def find_optimal_partition(self):
        """
        Find the optimal partition of nations into coalitions.
        Uses a greedy approach for computational efficiency.
        """
        remaining_nations = set(range(self.n))
        optimal_coalitions = []
        total_utility = 0
        
        while remaining_nations:
            best_coalition = None
            best_utility = -float('inf')
            
            # Try all possible coalitions from remaining nations
            for r in range(1, len(remaining_nations) + 1):
                for coalition in combinations(remaining_nations, r):
                    coalition_list = list(coalition)
                    utility = self.utility_function(coalition_list)
                    
                    if utility > best_utility:
                        best_utility = utility
                        best_coalition = coalition_list
            
            if best_coalition and best_utility > 0:
                optimal_coalitions.append(best_coalition)
                total_utility += best_utility
                remaining_nations -= set(best_coalition)
            else:
                # If no profitable coalition, nations go alone
                for nation in remaining_nations:
                    utility = self.utility_function([nation])
                    optimal_coalitions.append([nation])
                    total_utility += max(0, utility)
                break
        
        return optimal_coalitions, total_utility
    
    def analyze_coalition_stability(self, coalition):
        """
        Analyze the stability of a coalition using Shapley values.
        """
        coalition_size = len(coalition)
        shapley_values = np.zeros(coalition_size)
        
        # Calculate Shapley values
        for i, player in enumerate(coalition):
            marginal_contributions = []
            
            # Generate all possible orderings
            other_players = [p for p in coalition if p != player]
            
            for r in range(coalition_size):
                for predecessors in combinations(other_players, r):
                    pred_list = list(predecessors)
                    with_player = pred_list + [player]
                    
                    value_with = self.utility_function(with_player)
                    value_without = self.utility_function(pred_list) if pred_list else 0
                    
                    marginal_contribution = value_with - value_without
                    marginal_contributions.append(marginal_contribution)
            
            shapley_values[i] = np.mean(marginal_contributions)
        
        return shapley_values

# Example: 5 nations alliance optimization
nations = ['USA', 'UK', 'Germany', 'France', 'Japan']
security_weights = [10, 8, 7, 7, 6]  # Military/economic strength
cost_coefficients = [2, 1.8, 1.5, 1.6, 1.4]  # Cost efficiency
power_exponents = [1.2, 1.3, 1.1, 1.25, 1.15]  # Scaling factors

# Create optimizer instance
optimizer = AllianceOptimizer(nations, security_weights, cost_coefficients, power_exponents)

print("=== ALLIANCE OPTIMIZATION ANALYSIS ===\n")
print("Nations:", nations)
print("Security Weights:", security_weights)
print("Cost Coefficients:", cost_coefficients)
print("Power Exponents:", power_exponents)

# Find optimal partition
optimal_coalitions, total_utility = optimizer.find_optimal_partition()

print(f"\n=== OPTIMAL COALITION STRUCTURE ===")
print(f"Total System Utility: {total_utility:.2f}")
print("\nOptimal Coalitions:")
for i, coalition in enumerate(optimal_coalitions):
    coalition_nations = [nations[j] for j in coalition]
    utility = optimizer.utility_function(coalition)
    security = optimizer.security_value(coalition)
    cost = optimizer.total_cost(coalition)
    
    print(f"Coalition {i+1}: {coalition_nations}")
    print(f"  Utility: {utility:.2f} (Security: {security:.2f} - Cost: {cost:.2f})")

# Analyze individual coalition utilities for different sizes
print(f"\n=== COALITION SIZE ANALYSIS ===")
coalition_data = []

for size in range(1, len(nations) + 1):
    for coalition in combinations(range(len(nations)), size):
        coalition_list = list(coalition)
        utility = optimizer.utility_function(coalition_list)
        security = optimizer.security_value(coalition_list)
        cost = optimizer.total_cost(coalition_list)
        
        coalition_data.append({
            'Size': size,
            'Nations': [nations[i] for i in coalition_list],
            'Utility': utility,
            'Security': security,
            'Cost': cost,
            'Efficiency': utility/size if size > 0 else 0
        })

# Convert to DataFrame for analysis
df = pd.DataFrame(coalition_data)

# Show top 10 most efficient coalitions
print("\nTop 10 Most Efficient Coalitions (by Utility per Nation):")
top_coalitions = df.nlargest(10, 'Efficiency')
for _, row in top_coalitions.iterrows():
    print(f"{row['Nations']} - Utility: {row['Utility']:.2f}, Efficiency: {row['Efficiency']:.2f}")

# Stability analysis for the grand coalition
print(f"\n=== STABILITY ANALYSIS (Grand Coalition) ===")
grand_coalition = list(range(len(nations)))
shapley_values = optimizer.analyze_coalition_stability(grand_coalition)

print("Shapley Values (Fair Payoff Distribution):")
for i, (nation, value) in enumerate(zip(nations, shapley_values)):
    print(f"{nation}: {value:.2f}")

# Calculate synergy matrix visualization data
print(f"\n=== SYNERGY MATRIX ===")
print("Synergy between nations (how well they work together):")
synergy_df = pd.DataFrame(optimizer.synergy_matrix, 
                         index=nations, 
                         columns=nations)
print(synergy_df.round(2))

V_i(S) = w_i * Σ(synergy_ij for j ∈ S) * √|S|

C_i(S) = α_i * |S|^β_i

# Create comprehensive visualizations
plt.style.use('seaborn-v0_8')
fig = plt.figure(figsize=(20, 16))

# 1. Synergy Heatmap
ax1 = plt.subplot(3, 3, 1)
sns.heatmap(optimizer.synergy_matrix, 
            xticklabels=nations, 
            yticklabels=nations, 
            annot=True, 
            fmt='.2f', 
            cmap='RdYlBu_r',
            center=1.0,
            ax=ax1)
ax1.set_title('Nation Synergy Matrix\n(How Well Nations Work Together)', fontweight='bold')
ax1.set_xlabel('Nation')
ax1.set_ylabel('Nation')

# 2. Coalition Size vs Average Utility
ax2 = plt.subplot(3, 3, 2)
size_utilities = df.groupby('Size')['Utility'].agg(['mean', 'std', 'max']).reset_index()
ax2.errorbar(size_utilities['Size'], size_utilities['mean'], 
             yerr=size_utilities['std'], fmt='o-', capsize=5, linewidth=2)
ax2.plot(size_utilities['Size'], size_utilities['max'], 's--', 
         color='red', alpha=0.7, label='Maximum Utility')
ax2.set_xlabel('Coalition Size')
ax2.set_ylabel('Utility')
ax2.set_title('Coalition Size vs Utility\n(Error Bars Show Standard Deviation)', fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 3. Security vs Cost Trade-off
ax3 = plt.subplot(3, 3, 3)
scatter = ax3.scatter(df['Cost'], df['Security'], 
                     c=df['Utility'], s=df['Size']*20,
                     cmap='viridis', alpha=0.7)
ax3.set_xlabel('Total Cost')
ax3.set_ylabel('Security Value')
ax3.set_title('Security-Cost Trade-off\n(Size=Bubble Size, Color=Utility)', fontweight='bold')
plt.colorbar(scatter, ax=ax3, label='Utility')

# 4. Efficiency Distribution by Coalition Size
ax4 = plt.subplot(3, 3, 4)
df_positive = df[df['Utility'] > 0]  # Only profitable coalitions
box_data = [df_positive[df_positive['Size'] == size]['Efficiency'].values 
            for size in range(1, 6)]
box_plot = ax4.boxplot(box_data, labels=range(1, 6), patch_artist=True)
colors = plt.cm.Set3(np.linspace(0, 1, 5))
for patch, color in zip(box_plot['boxes'], colors):
    patch.set_facecolor(color)
ax4.set_xlabel('Coalition Size')
ax4.set_ylabel('Efficiency (Utility per Nation)')
ax4.set_title('Efficiency Distribution by Coalition Size', fontweight='bold')

# 5. Shapley Values Bar Chart
ax5 = plt.subplot(3, 3, 5)
bars = ax5.bar(nations, shapley_values, color=plt.cm.viridis(np.linspace(0, 1, len(nations))))
ax5.set_xlabel('Nation')
ax5.set_ylabel('Shapley Value')
ax5.set_title('Fair Payoff Distribution\n(Shapley Values for Grand Coalition)', fontweight='bold')
ax5.tick_params(axis='x', rotation=45)

# Add value labels on bars
for bar, value in zip(bars, shapley_values):
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height + 0.01,
             f'{value:.1f}', ha='center', va='bottom', fontweight='bold')

# 6. Optimal Coalition Structure Network
ax6 = plt.subplot(3, 3, 6)
G = nx.Graph()

# Add nodes for each nation
for i, nation in enumerate(nations):
    G.add_node(nation, weight=security_weights[i])

# Add edges within each optimal coalition
colors = ['red', 'blue', 'green', 'orange', 'purple']
coalition_colors = {}

for i, coalition in enumerate(optimal_coalitions):
    color = colors[i % len(colors)]
    coalition_nations = [nations[j] for j in coalition]
    coalition_colors.update({nation: color for nation in coalition_nations})
    
    # Add edges between coalition members
    for j in range(len(coalition)):
        for k in range(j+1, len(coalition)):
            G.add_edge(nations[coalition[j]], nations[coalition[k]])

# Draw network
pos = nx.spring_layout(G, seed=42)
node_colors = [coalition_colors[nation] for nation in G.nodes()]
node_sizes = [security_weights[nations.index(nation)] * 100 for nation in G.nodes()]

nx.draw(G, pos, ax=ax6, node_color=node_colors, node_size=node_sizes,
        with_labels=True, font_size=10, font_weight='bold')
ax6.set_title('Optimal Coalition Structure\n(Colors Show Coalition Membership)', fontweight='bold')

# 7. Cost Structure Analysis
ax7 = plt.subplot(3, 3, 7)
coalition_sizes = range(1, 6)
for i, nation in enumerate(nations):
    costs = [optimizer.cost_function(i, size) for size in coalition_sizes]
    ax7.plot(coalition_sizes, costs, 'o-', label=nation, linewidth=2)

ax7.set_xlabel('Coalition Size')
ax7.set_ylabel('Individual Cost')
ax7.set_title('Cost Structure by Nation\n(How Costs Scale with Coalition Size)', fontweight='bold')
ax7.legend()
ax7.grid(True, alpha=0.3)

# 8. Utility Breakdown for Optimal Coalitions
ax8 = plt.subplot(3, 3, 8)
coalition_labels = []
utilities = []
securities = []
costs = []

for i, coalition in enumerate(optimal_coalitions):
    coalition_nations = [nations[j] for j in coalition]
    label = f"C{i+1}: {', '.join(coalition_nations)}"
    if len(label) > 20:
        label = f"C{i+1}: {len(coalition)} nations"
    
    coalition_labels.append(label)
    utilities.append(optimizer.utility_function(coalition))
    securities.append(optimizer.security_value(coalition))
    costs.append(optimizer.total_cost(coalition))

x = np.arange(len(coalition_labels))
width = 0.25

bars1 = ax8.bar(x - width, securities, width, label='Security', alpha=0.8)
bars2 = ax8.bar(x, costs, width, label='Cost', alpha=0.8)
bars3 = ax8.bar(x + width, utilities, width, label='Utility', alpha=0.8)

ax8.set_xlabel('Coalition')
ax8.set_ylabel('Value')
ax8.set_title('Utility Breakdown for Optimal Coalitions', fontweight='bold')
ax8.set_xticks(x)
ax8.set_xticklabels(coalition_labels, rotation=45, ha='right')
ax8.legend()

# 9. Marginal Utility Analysis
ax9 = plt.subplot(3, 3, 9)
marginal_utilities = []
base_coalition = [0]  # Start with USA

for i in range(1, len(nations)):
    current_utility = optimizer.utility_function(base_coalition)
    expanded_coalition = base_coalition + [i]
    expanded_utility = optimizer.utility_function(expanded_coalition)
    marginal_utility = expanded_utility - current_utility
    marginal_utilities.append(marginal_utility)
    base_coalition = expanded_coalition

nation_additions = nations[1:]
bars = ax9.bar(nation_additions, marginal_utilities, 
               color=plt.cm.plasma(np.linspace(0, 1, len(marginal_utilities))))
ax9.set_xlabel('Nation Added to Coalition')
ax9.set_ylabel('Marginal Utility')
ax9.set_title('Marginal Utility of Adding Nations\n(Starting with USA)', fontweight='bold')
ax9.tick_params(axis='x', rotation=45)

# Add value labels on bars
for bar, value in zip(bars, marginal_utilities):
    height = bar.get_height()
    ax9.text(bar.get_x() + bar.get_width()/2., height + (0.01 if height >= 0 else -0.05),
             f'{value:.1f}', ha='center', va='bottom' if height >= 0 else 'top', 
             fontweight='bold')

plt.tight_layout()
plt.show()

# Additional Analysis: Coalition Formation Process
print("\n" + "="*60)
print("DETAILED COALITION FORMATION ANALYSIS")
print("="*60)

# Show step-by-step coalition formation
print("\nStep-by-step Coalition Formation Process:")
remaining = set(range(len(nations)))
step = 1

temp_coalitions = []
while remaining:
    print(f"\nStep {step}:")
    print(f"Remaining nations: {[nations[i] for i in remaining]}")
    
    best_coalition = None
    best_utility = -float('inf')
    
    # Try all possible coalitions from remaining nations
    for r in range(1, len(remaining) + 1):
        for coalition in combinations(remaining, r):
            coalition_list = list(coalition)
            utility = optimizer.utility_function(coalition_list)
            
            if utility > best_utility:
                best_utility = utility
                best_coalition = coalition_list
    
    if best_coalition and best_utility > 0:
        coalition_nations = [nations[i] for i in best_coalition]
        print(f"Best coalition found: {coalition_nations}")
        print(f"Utility: {best_utility:.2f}")
        temp_coalitions.append(best_coalition)
        remaining -= set(best_coalition)
        step += 1
    else:
        print("No profitable coalitions found. Remaining nations form individual coalitions.")
        for nation in remaining:
            temp_coalitions.append([nation])
        break

# Performance metrics
print(f"\n" + "="*40)
print("PERFORMANCE METRICS")
print("="*40)

total_security = sum(optimizer.security_value(coalition) for coalition in optimal_coalitions)
total_cost = sum(optimizer.total_cost(coalition) for coalition in optimal_coalitions)
efficiency = total_utility / len(nations)

print(f"Total System Security: {total_security:.2f}")
print(f"Total System Cost: {total_cost:.2f}")
print(f"Net System Utility: {total_utility:.2f}")
print(f"Average Utility per Nation: {efficiency:.2f}")
print(f"System Efficiency: {(total_utility/total_security)*100:.1f}%")

# Compare with alternative structures
print(f"\nComparison with Alternative Structures:")

# All individual (no alliances)
individual_utility = sum(max(0, optimizer.utility_function([i])) for i in range(len(nations)))
print(f"All Individual Coalitions Utility: {individual_utility:.2f}")

# Grand coalition (everyone together)
grand_utility = optimizer.utility_function(list(range(len(nations))))
print(f"Grand Coalition Utility: {grand_utility:.2f}")

print(f"\nOptimal Structure Advantage:")
print(f"vs Individual: +{total_utility - individual_utility:.2f} ({((total_utility/individual_utility)-1)*100:.1f}% improvement)")
print(f"vs Grand Coalition: +{total_utility - grand_utility:.2f} ({((total_utility/grand_utility)-1)*100:.1f}% improvement)")

# Coalition stability analysis
print(f"\n" + "="*40)
print("COALITION STABILITY ANALYSIS")
print("="*40)

for i, coalition in enumerate(optimal_coalitions):
    if len(coalition) > 1:
        coalition_nations = [nations[j] for j in coalition]
        print(f"\nCoalition {i+1}: {coalition_nations}")
        
        shapley_vals = optimizer.analyze_coalition_stability(coalition)
        coalition_utility = optimizer.utility_function(coalition)
        
        print(f"Total Coalition Utility: {coalition_utility:.2f}")
        print("Fair payoff distribution (Shapley values):")
        for j, (nation_idx, value) in enumerate(zip(coalition, shapley_vals)):
            nation_name = nations[nation_idx]
            percentage = (value / coalition_utility) * 100 if coalition_utility > 0 else 0
            print(f"  {nation_name}: {value:.2f} ({percentage:.1f}%)")
            
        # Check for stability (no negative Shapley values)
        if all(val >= 0 for val in shapley_vals):
            print("  ✓ Coalition is stable (all members benefit)")
        else:
            unstable_members = [nations[coalition[j]] for j, val in enumerate(shapley_vals) if val < 0]
            print(f"  ⚠ Potentially unstable: {unstable_members} might leave")

print("\n" + "="*60)
print("ANALYSIS COMPLETE")
print("="*60)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
import seaborn as sns
from matplotlib.patches import Wedge
import pandas as pd

# Set up the optimization problem
# Objective function coefficients (we need to negate for maximization)
# Original: Maximize 0.35*x1 + 0.25*x2 + 0.20*x3 + 0.20*x4
c = [-0.35, -0.25, -0.20, -0.20]

# Inequality constraints (Ax <= b)
# Budget constraint: x1 + x2 + x3 + x4 <= 10
# Strategic balance: x1 <= 5
A_ub = [[1, 1, 1, 1],    # Total budget constraint
        [1, 0, 0, 0]]    # Air defense limit

b_ub = [10, 5]

# Bounds for each variable (min, max)
# x1 >= 2, x2 >= 1.5, x3 >= 1, x4 >= 0.5
bounds = [(2.0, None),    # Air Defense
          (1.5, None),    # Naval Forces  
          (1.0, None),    # Ground Forces
          (0.5, None)]    # Intelligence & Cyber

# Solve the linear programming problem
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

print("=== DEFENSE BUDGET OPTIMIZATION RESULTS ===\n")
print(f"Optimization Status: {result.message}")
print(f"Total Defense Effectiveness Score: {-result.fun:.3f}")
print("\nOptimal Budget Allocation (in billions USD):")
print(f"Air Defense Systems:     ${result.x[0]:.2f}B ({result.x[0]/10*100:.1f}%)")
print(f"Naval Forces:            ${result.x[1]:.2f}B ({result.x[1]/10*100:.1f}%)")
print(f"Ground Forces:           ${result.x[2]:.2f}B ({result.x[2]/10*100:.1f}%)")
print(f"Intelligence & Cyber:    ${result.x[3]:.2f}B ({result.x[3]/10*100:.1f}%)")
print(f"Total Budget Used:       ${sum(result.x):.2f}B")

# Store results for visualization
categories = ['Air Defense', 'Naval Forces', 'Ground Forces', 'Intelligence & Cyber']
allocations = result.x
effectiveness_coeffs = [0.35, 0.25, 0.20, 0.20]
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4']

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

# 1. Budget Allocation Pie Chart
ax1 = plt.subplot(2, 3, 1)
wedges, texts, autotexts = ax1.pie(allocations, labels=categories, colors=colors, 
                                   autopct='%1.1f%%', startangle=90, 
                                   textprops={'fontsize': 10})
ax1.set_title('Optimal Budget Allocation\n($10B Total)', fontsize=14, fontweight='bold')

# 2. Budget vs Minimum Requirements
ax2 = plt.subplot(2, 3, 2)
min_requirements = [2.0, 1.5, 1.0, 0.5]
x_pos = np.arange(len(categories))

bars1 = ax2.bar(x_pos - 0.2, min_requirements, 0.4, label='Minimum Required', 
                color='lightcoral', alpha=0.7)
bars2 = ax2.bar(x_pos + 0.2, allocations, 0.4, label='Optimal Allocation', 
                color=colors, alpha=0.8)

ax2.set_xlabel('Defense Categories')
ax2.set_ylabel('Budget (Billions USD)')
ax2.set_title('Optimal vs Minimum Requirements', fontweight='bold')
ax2.set_xticks(x_pos)
ax2.set_xticklabels(categories, rotation=45, ha='right')
ax2.legend()
ax2.grid(axis='y', alpha=0.3)

# Add value labels on bars
for bar in bars2:
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height + 0.05,
             f'${height:.2f}B', ha='center', va='bottom', fontsize=9)

# 3. Effectiveness Contribution
ax3 = plt.subplot(2, 3, 3)
effectiveness_contributions = [allocations[i] * effectiveness_coeffs[i] for i in range(4)]
bars3 = ax3.bar(categories, effectiveness_contributions, color=colors, alpha=0.8)
ax3.set_xlabel('Defense Categories')
ax3.set_ylabel('Effectiveness Contribution')
ax3.set_title('Defense Effectiveness by Category', fontweight='bold')
ax3.tick_params(axis='x', rotation=45)
ax3.grid(axis='y', alpha=0.3)

# Add value labels
for bar in bars3:
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height + 0.01,
             f'{height:.3f}', ha='center', va='bottom', fontsize=9)

# 4. Resource Efficiency Analysis
ax4 = plt.subplot(2, 3, 4)
efficiency_ratios = [effectiveness_contributions[i]/allocations[i] for i in range(4)]
bars4 = ax4.bar(categories, efficiency_ratios, color=colors, alpha=0.8)
ax4.set_xlabel('Defense Categories')
ax4.set_ylabel('Effectiveness per Dollar')
ax4.set_title('Resource Efficiency Analysis', fontweight='bold')
ax4.tick_params(axis='x', rotation=45)
ax4.grid(axis='y', alpha=0.3)

# Add value labels
for bar in bars4:
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height + 0.001,
             f'{height:.3f}', ha='center', va='bottom', fontsize=9)

# 5. Budget Utilization Summary
ax5 = plt.subplot(2, 3, 5)
budget_data = pd.DataFrame({
    'Category': categories,
    'Allocated': allocations,
    'Minimum': min_requirements,
    'Excess': [allocations[i] - min_requirements[i] for i in range(4)]
})

x_pos = np.arange(len(categories))
ax5.bar(x_pos, budget_data['Minimum'], color='lightgray', label='Minimum Required')
ax5.bar(x_pos, budget_data['Excess'], bottom=budget_data['Minimum'], 
        color=colors, alpha=0.8, label='Additional Allocation')

ax5.set_xlabel('Defense Categories')
ax5.set_ylabel('Budget (Billions USD)')
ax5.set_title('Budget Composition Analysis', fontweight='bold')
ax5.set_xticks(x_pos)
ax5.set_xticklabels(categories, rotation=45, ha='right')
ax5.legend()
ax5.grid(axis='y', alpha=0.3)

# 6. Sensitivity Analysis - What if budget changes?
ax6 = plt.subplot(2, 3, 6)
budget_scenarios = np.linspace(7, 15, 20)
total_effectiveness = []

for budget in budget_scenarios:
    # Solve for each budget scenario
    A_ub_temp = [[1, 1, 1, 1], [1, 0, 0, 0]]
    b_ub_temp = [budget, min(budget * 0.5, 5)]
    
    try:
        result_temp = linprog(c, A_ub=A_ub_temp, b_ub=b_ub_temp, bounds=bounds, method='highs')
        if result_temp.success:
            total_effectiveness.append(-result_temp.fun)
        else:
            total_effectiveness.append(np.nan)
    except:
        total_effectiveness.append(np.nan)

ax6.plot(budget_scenarios, total_effectiveness, 'o-', color='darkblue', linewidth=2, markersize=4)
ax6.axvline(x=10, color='red', linestyle='--', alpha=0.7, label='Current Budget')
ax6.set_xlabel('Total Budget (Billions USD)')
ax6.set_ylabel('Total Defense Effectiveness')
ax6.set_title('Budget Sensitivity Analysis', fontweight='bold')
ax6.grid(True, alpha=0.3)
ax6.legend()

plt.tight_layout()
plt.show()

# Print detailed analysis
print("\n=== DETAILED ANALYSIS ===")
print("\n1. RESOURCE EFFICIENCY:")
for i, category in enumerate(categories):
    efficiency = effectiveness_contributions[i] / allocations[i]
    print(f"{category:20}: {efficiency:.3f} effectiveness per $1B")

print("\n2. BUDGET UTILIZATION:")
total_minimum = sum(min_requirements)
excess_budget = 10 - total_minimum
print(f"Total Minimum Required: ${total_minimum:.1f}B ({total_minimum/10*100:.1f}%)")
print(f"Excess Budget Available: ${excess_budget:.1f}B ({excess_budget/10*100:.1f}%)")

print("\n3. STRATEGIC INSIGHTS:")
print(f"• Air Defense receives highest allocation due to effectiveness coefficient (0.35)")
print(f"• Intelligence & Cyber gets minimal funding above requirement")
print(f"• Ground Forces receives moderate boost above minimum")
print(f"• Naval Forces allocation balances cost and strategic importance")

print(f"\n4. CONSTRAINT STATUS:")
print(f"• Budget Constraint: ${sum(allocations):.2f}B used of $10B available")
print(f"• Air Defense Limit: ${allocations[0]:.2f}B allocated (max ${min(10*0.5, 5):.1f}B)")
for i, category in enumerate(categories):
    surplus = allocations[i] - min_requirements[i]
    print(f"• {category}: ${surplus:.2f}B above minimum requirement")

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
from scipy.optimize import minimize
import seaborn as sns
from itertools import combinations
import warnings
warnings.filterwarnings('ignore')

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

class MilitaryBaseOptimizer:
    def __init__(self, potential_bases, critical_areas, base_costs, max_distance=50):
        """
        Initialize the military base optimization problem
        
        Parameters:
        - potential_bases: List of (x, y) coordinates for potential base locations
        - critical_areas: List of (x, y) coordinates for areas that need coverage
        - base_costs: List of costs for establishing each base
        - max_distance: Maximum effective range for base coverage
        """
        self.potential_bases = np.array(potential_bases)
        self.critical_areas = np.array(critical_areas)
        self.base_costs = np.array(base_costs)
        self.max_distance = max_distance
        
        # Calculate distance matrix between bases and critical areas
        self.distance_matrix = cdist(self.potential_bases, self.critical_areas)
        
        # Create coverage matrix (1 if base can cover area, 0 otherwise)
        self.coverage_matrix = (self.distance_matrix <= max_distance).astype(int)
        
        self.n_bases = len(potential_bases)
        self.n_areas = len(critical_areas)
        
    def greedy_solution(self):
        """
        Solve using greedy algorithm: select base that covers most uncovered areas
        """
        selected_bases = []
        covered_areas = set()
        total_cost = 0
        
        while len(covered_areas) < self.n_areas:
            best_base = -1
            best_new_coverage = 0
            best_efficiency = 0
            
            for i in range(self.n_bases):
                if i in selected_bases:
                    continue
                    
                # Count new areas this base would cover
                coverable_areas = set(np.where(self.coverage_matrix[i] == 1)[0])
                new_coverage = len(coverable_areas - covered_areas)
                
                # Calculate efficiency (coverage per cost)
                if new_coverage > 0:
                    efficiency = new_coverage / self.base_costs[i]
                    if efficiency > best_efficiency:
                        best_efficiency = efficiency
                        best_base = i
                        best_new_coverage = new_coverage
            
            if best_base == -1:
                break
                
            selected_bases.append(best_base)
            covered_areas.update(np.where(self.coverage_matrix[best_base] == 1)[0])
            total_cost += self.base_costs[best_base]
        
        return selected_bases, total_cost, covered_areas
    
    def brute_force_optimal(self):
        """
        Find optimal solution by checking all possible combinations (for small instances)
        """
        best_cost = float('inf')
        best_solution = None
        
        # Try all possible combinations of bases
        for r in range(1, min(self.n_bases + 1, 8)):  # Limit to prevent exponential explosion
            for base_combo in combinations(range(self.n_bases), r):
                # Check if this combination covers all areas
                covered = set()
                for base_idx in base_combo:
                    covered.update(np.where(self.coverage_matrix[base_idx] == 1)[0])
                
                if len(covered) == self.n_areas:
                    # Calculate total cost
                    cost = sum(self.base_costs[i] for i in base_combo)
                    if cost < best_cost:
                        best_cost = cost
                        best_solution = list(base_combo)
        
        return best_solution, best_cost
    
    def plot_solution(self, selected_bases, title="Military Base Deployment"):
        """
        Visualize the base deployment solution
        """
        plt.figure(figsize=(12, 8))
        
        # Plot critical areas
        plt.scatter(self.critical_areas[:, 0], self.critical_areas[:, 1], 
                   c='red', s=100, alpha=0.7, marker='s', label='Critical Areas')
        
        # Plot potential bases (not selected)
        unselected = [i for i in range(self.n_bases) if i not in selected_bases]
        if unselected:
            plt.scatter(self.potential_bases[unselected, 0], 
                       self.potential_bases[unselected, 1], 
                       c='lightgray', s=150, alpha=0.5, marker='^', 
                       label='Unselected Base Sites')
        
        # Plot selected bases
        if selected_bases:
            plt.scatter(self.potential_bases[selected_bases, 0], 
                       self.potential_bases[selected_bases, 1], 
                       c='blue', s=200, marker='^', label='Selected Bases')
            
            # Draw coverage circles
            for base_idx in selected_bases:
                circle = plt.Circle((self.potential_bases[base_idx, 0], 
                                   self.potential_bases[base_idx, 1]), 
                                  self.max_distance, 
                                  fill=False, color='blue', alpha=0.3, linestyle='--')
                plt.gca().add_patch(circle)
        
        plt.xlabel('X Coordinate (km)')
        plt.ylabel('Y Coordinate (km)')
        plt.title(title)
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.axis('equal')
        plt.tight_layout()
        plt.show()
    
    def analyze_coverage(self, selected_bases):
        """
        Analyze coverage statistics for the solution
        """
        coverage_stats = {}
        
        # Calculate which areas are covered by which bases
        area_coverage = {}
        for area_idx in range(self.n_areas):
            covering_bases = []
            for base_idx in selected_bases:
                if self.coverage_matrix[base_idx, area_idx] == 1:
                    covering_bases.append(base_idx)
            area_coverage[area_idx] = covering_bases
        
        # Coverage redundancy analysis
        redundancy_levels = [len(bases) for bases in area_coverage.values()]
        coverage_stats['avg_redundancy'] = np.mean(redundancy_levels)
        coverage_stats['min_redundancy'] = min(redundancy_levels) if redundancy_levels else 0
        coverage_stats['max_redundancy'] = max(redundancy_levels) if redundancy_levels else 0
        
        # Base utilization
        base_utilization = {}
        for base_idx in selected_bases:
            covered_areas = np.sum(self.coverage_matrix[base_idx])
            base_utilization[base_idx] = covered_areas
        
        coverage_stats['base_utilization'] = base_utilization
        
        return coverage_stats, area_coverage

# Generate sample data for our military base optimization problem
def generate_scenario():
    """Generate a realistic military base placement scenario"""
    
    # Define potential base locations (strategic points)
    potential_bases = [
        (10, 10), (30, 15), (50, 20), (70, 25),
        (15, 40), (35, 45), (55, 50), (75, 55),
        (20, 70), (40, 75), (60, 80), (80, 85)
    ]
    
    # Define critical areas that need protection
    critical_areas = [
        (25, 30), (45, 35), (65, 40),  # Population centers
        (20, 60), (50, 65), (70, 70),  # Industrial areas
        (35, 25), (55, 30),            # Infrastructure
        (40, 55), (60, 60)             # Strategic resources
    ]
    
    # Base establishment costs (in millions)
    base_costs = [100, 120, 95, 110, 130, 105, 125, 115, 140, 108, 122, 135]
    
    return potential_bases, critical_areas, base_costs

# Main execution
if __name__ == "__main__":
    # Generate scenario
    potential_bases, critical_areas, base_costs = generate_scenario()
    
    # Initialize optimizer
    optimizer = MilitaryBaseOptimizer(potential_bases, critical_areas, base_costs, max_distance=35)
    
    print("=== MILITARY BASE OPTIMIZATION PROBLEM ===")
    print(f"Number of potential base locations: {optimizer.n_bases}")
    print(f"Number of critical areas to protect: {optimizer.n_areas}")
    print(f"Maximum base coverage range: {optimizer.max_distance} km")
    print(f"Base establishment costs: {base_costs}")
    print()
    
    # Solve using greedy algorithm
    print("--- GREEDY ALGORITHM SOLUTION ---")
    greedy_bases, greedy_cost, greedy_covered = optimizer.greedy_solution()
    print(f"Selected bases (greedy): {greedy_bases}")
    print(f"Total cost (greedy): ${greedy_cost} million")
    print(f"Areas covered: {len(greedy_covered)}/{optimizer.n_areas}")
    print()
    
    # Solve using brute force (optimal for small instances)
    print("--- OPTIMAL SOLUTION (Brute Force) ---")
    optimal_bases, optimal_cost = optimizer.brute_force_optimal()
    if optimal_bases:
        print(f"Selected bases (optimal): {optimal_bases}")
        print(f"Total cost (optimal): ${optimal_cost} million")
        print(f"Cost savings vs greedy: ${greedy_cost - optimal_cost} million")
    else:
        print("No feasible solution found with brute force method")
    print()
    
    # Analyze coverage for optimal solution
    if optimal_bases:
        coverage_stats, area_coverage = optimizer.analyze_coverage(optimal_bases)
        print("--- COVERAGE ANALYSIS ---")
        print(f"Average coverage redundancy: {coverage_stats['avg_redundancy']:.2f}")
        print(f"Minimum coverage redundancy: {coverage_stats['min_redundancy']}")
        print(f"Maximum coverage redundancy: {coverage_stats['max_redundancy']}")
        print("\nBase utilization (areas covered per base):")
        for base_idx, util in coverage_stats['base_utilization'].items():
            print(f"  Base {base_idx}: covers {util} areas (cost: ${base_costs[base_idx]}M)")
    
    # Visualize both solutions
    optimizer.plot_solution(greedy_bases, "Greedy Algorithm Solution")
    if optimal_bases:
        optimizer.plot_solution(optimal_bases, "Optimal Solution")
    
    # Create comparison chart
    plt.figure(figsize=(10, 6))
    methods = ['Greedy', 'Optimal']
    costs = [greedy_cost, optimal_cost if optimal_bases else greedy_cost]
    colors = ['orange', 'green']
    
    bars = plt.bar(methods, costs, color=colors, alpha=0.7)
    plt.title('Cost Comparison: Greedy vs Optimal Solution')
    plt.ylabel('Total Cost (Million $)')
    plt.grid(True, alpha=0.3)
    
    # Add value labels on bars
    for bar, cost in zip(bars, costs):
        plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2,
                f'${cost}M', ha='center', va='bottom', fontweight='bold')
    
    plt.tight_layout()
    plt.show()
    
    # Distance matrix heatmap
    plt.figure(figsize=(12, 8))
    plt.subplot(1, 2, 1)
    sns.heatmap(optimizer.distance_matrix, annot=True, fmt='.1f', cmap='YlOrRd',
                xticklabels=[f'Area {i}' for i in range(optimizer.n_areas)],
                yticklabels=[f'Base {i}' for i in range(optimizer.n_bases)])
    plt.title('Distance Matrix (km)')
    plt.ylabel('Potential Bases')
    plt.xlabel('Critical Areas')
    
    plt.subplot(1, 2, 2)
    sns.heatmap(optimizer.coverage_matrix, annot=True, fmt='d', cmap='RdYlGn',
                xticklabels=[f'Area {i}' for i in range(optimizer.n_areas)],
                yticklabels=[f'Base {i}' for i in range(optimizer.n_bases)])
    plt.title('Coverage Matrix (1=Can Cover, 0=Cannot)')
    plt.ylabel('Potential Bases')
    plt.xlabel('Critical Areas')
    
    plt.tight_layout()
    plt.show()

efficiency = new_coverage / self.base_costs[i]

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from scipy.optimize import minimize
from itertools import combinations
import pandas as pd
import seaborn as sns
from sklearn.cluster import KMeans
import warnings
warnings.filterwarnings('ignore')

class EconomicZoneOptimizer:
    """
    A class to optimize economic zone development through mathematical modeling.
    This includes facility location, network design, and resource allocation.
    """
    
    def __init__(self, cities, populations, gdp_per_capita):
        """
        Initialize the optimizer with city data.
        
        Parameters:
        cities (list): List of city names
        populations (list): Population of each city
        gdp_per_capita (list): GDP per capita for each city
        """
        self.cities = cities
        self.populations = np.array(populations)
        self.gdp_per_capita = np.array(gdp_per_capita)
        self.n_cities = len(cities)
        
        # Generate random coordinates for cities (in practice, use real coordinates)
        np.random.seed(42)
        self.coordinates = np.random.rand(self.n_cities, 2) * 100
        
        # Calculate distance matrix
        self.distance_matrix = self._calculate_distances()
        
        # Calculate economic potential for each city
        self.economic_potential = self.populations * self.gdp_per_capita
        
    def _calculate_distances(self):
        """Calculate Euclidean distances between all pairs of cities."""
        distances = np.zeros((self.n_cities, self.n_cities))
        for i in range(self.n_cities):
            for j in range(self.n_cities):
                distances[i, j] = np.sqrt(
                    (self.coordinates[i, 0] - self.coordinates[j, 0])**2 +
                    (self.coordinates[i, 1] - self.coordinates[j, 1])**2
                )
        return distances
    
    def calculate_infrastructure_cost(self, distance, population_i, population_j):
        """
        Calculate infrastructure cost between two cities.
        Cost increases with distance but decreases with economic benefit.
        """
        base_cost = distance * 1000  # Base cost per unit distance
        economic_factor = np.sqrt(population_i * population_j) / 10000  # Economic benefit
        return base_cost / (1 + economic_factor)
    
    def optimize_hub_selection(self, max_hubs=5):
        """
        Optimize hub selection using K-means clustering and economic potential.
        """
        # Use K-means to identify geographical clusters
        kmeans = KMeans(n_clusters=min(max_hubs, self.n_cities), random_state=42)
        clusters = kmeans.fit_predict(self.coordinates)
        
        # Select hub for each cluster based on economic potential
        hubs = []
        hub_indices = []
        
        for cluster_id in range(max(clusters) + 1):
            cluster_cities = np.where(clusters == cluster_id)[0]
            # Select city with highest economic potential in this cluster
            best_city = cluster_cities[np.argmax(self.economic_potential[cluster_cities])]
            hubs.append(self.cities[best_city])
            hub_indices.append(best_city)
        
        return hub_indices, hubs
    
    def optimize_network_topology(self, hub_indices):
        """
        Optimize network connections between hubs using minimum spanning tree approach.
        """
        n_hubs = len(hub_indices)
        
        # Calculate cost matrix for hub connections
        hub_costs = np.zeros((n_hubs, n_hubs))
        for i in range(n_hubs):
            for j in range(n_hubs):
                if i != j:
                    city_i, city_j = hub_indices[i], hub_indices[j]
                    hub_costs[i, j] = self.calculate_infrastructure_cost(
                        self.distance_matrix[city_i, city_j],
                        self.populations[city_i],
                        self.populations[city_j]
                    )
        
        # Create graph and find minimum spanning tree
        G = nx.Graph()
        for i in range(n_hubs):
            for j in range(i + 1, n_hubs):
                G.add_edge(i, j, weight=hub_costs[i, j])
        
        # Find minimum spanning tree
        mst = nx.minimum_spanning_tree(G, weight='weight')
        
        return mst, hub_costs
    
    def calculate_trade_efficiency(self, hub_indices, mst):
        """
        Calculate trade efficiency metrics for the optimized network.
        """
        # Calculate average path length in the network
        path_lengths = []
        total_trade_volume = 0
        
        for i in range(len(hub_indices)):
            for j in range(i + 1, len(hub_indices)):
                try:
                    path_length = nx.shortest_path_length(mst, i, j, weight='weight')
                    path_lengths.append(path_length)
                    
                    # Estimate trade volume based on economic potential and distance
                    city_i, city_j = hub_indices[i], hub_indices[j]
                    trade_volume = (self.economic_potential[city_i] * 
                                  self.economic_potential[city_j]) / path_length
                    total_trade_volume += trade_volume
                except nx.NetworkXNoPath:
                    # No path exists
                    path_lengths.append(float('inf'))
        
        avg_path_length = np.mean([p for p in path_lengths if p != float('inf')])
        efficiency_score = total_trade_volume / (avg_path_length + 1)
        
        return {
            'average_path_length': avg_path_length,
            'total_trade_volume': total_trade_volume,
            'efficiency_score': efficiency_score,
            'connectivity': len([p for p in path_lengths if p != float('inf')]) / len(path_lengths)
        }
    
    def run_optimization(self, max_hubs=5):
        """
        Run the complete optimization process.
        """
        print("🏗️ Starting Economic Zone Optimization...")
        
        # Step 1: Optimize hub selection
        print("📍 Optimizing hub selection...")
        hub_indices, hub_cities = self.optimize_hub_selection(max_hubs)
        
        # Step 2: Optimize network topology
        print("🌐 Optimizing network topology...")
        mst, hub_costs = self.optimize_network_topology(hub_indices)
        
        # Step 3: Calculate efficiency metrics
        print("📊 Calculating efficiency metrics...")
        efficiency_metrics = self.calculate_trade_efficiency(hub_indices, mst)
        
        # Calculate total infrastructure cost
        total_cost = sum([mst[u][v]['weight'] for u, v in mst.edges()])
        
        results = {
            'hub_indices': hub_indices,
            'hub_cities': hub_cities,
            'network': mst,
            'hub_costs': hub_costs,
            'total_cost': total_cost,
            'efficiency_metrics': efficiency_metrics
        }
        
        print("✅ Optimization complete!")
        return results
    
    def visualize_results(self, results):
        """
        Create comprehensive visualizations of the optimization results.
        """
        fig = plt.figure(figsize=(20, 15))
        
        # Plot 1: Economic Zone Network Map
        ax1 = plt.subplot(2, 3, 1)
        
        # Plot all cities
        scatter = ax1.scatter(self.coordinates[:, 0], self.coordinates[:, 1], 
                             c=self.economic_potential, s=self.populations/1000,
                             alpha=0.6, cmap='viridis', edgecolors='black')
        
        # Highlight hubs
        hub_coords = self.coordinates[results['hub_indices']]
        ax1.scatter(hub_coords[:, 0], hub_coords[:, 1], 
                   c='red', s=200, marker='*', edgecolors='black', linewidth=2,
                   label='Economic Hubs')
        
        # Draw network connections
        mst = results['network']
        for edge in mst.edges():
            i, j = results['hub_indices'][edge[0]], results['hub_indices'][edge[1]]
            ax1.plot([self.coordinates[i, 0], self.coordinates[j, 0]],
                    [self.coordinates[i, 1], self.coordinates[j, 1]], 
                    'r-', linewidth=2, alpha=0.7)
        
        # Add city labels
        for i, city in enumerate(self.cities):
            ax1.annotate(city, (self.coordinates[i, 0], self.coordinates[i, 1]), 
                        xytext=(5, 5), textcoords='offset points', fontsize=8)
        
        ax1.set_title('Optimized Economic Zone Network', fontsize=14, fontweight='bold')
        ax1.set_xlabel('X Coordinate (km)')
        ax1.set_ylabel('Y Coordinate (km)')
        plt.colorbar(scatter, ax=ax1, label='Economic Potential')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # Plot 2: Cost Analysis
        ax2 = plt.subplot(2, 3, 2)
        
        # Extract edge costs for visualization
        edge_costs = [mst[u][v]['weight'] for u, v in mst.edges()]
        edge_labels = [f"Hub {u+1}-{v+1}" for u, v in mst.edges()]
        
        bars = ax2.bar(range(len(edge_costs)), edge_costs, color='skyblue', edgecolor='navy')
        ax2.set_title('Infrastructure Costs by Connection', fontsize=14, fontweight='bold')
        ax2.set_xlabel('Network Connections')
        ax2.set_ylabel('Cost (Million $)')
        ax2.set_xticks(range(len(edge_labels)))
        ax2.set_xticklabels(edge_labels, rotation=45)
        
        # Add value labels on bars
        for bar, cost in zip(bars, edge_costs):
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + max(edge_costs)*0.01,
                    f'${cost:.1f}M', ha='center', va='bottom', fontweight='bold')
        
        ax2.grid(True, alpha=0.3, axis='y')
        
        # Plot 3: Economic Potential Comparison
        ax3 = plt.subplot(2, 3, 3)
        
        hub_potential = self.economic_potential[results['hub_indices']]
        non_hub_potential = np.delete(self.economic_potential, results['hub_indices'])
        
        ax3.hist([hub_potential, non_hub_potential], bins=10, alpha=0.7, 
                label=['Hub Cities', 'Non-Hub Cities'], color=['red', 'blue'])
        ax3.set_title('Economic Potential Distribution', fontsize=14, fontweight='bold')
        ax3.set_xlabel('Economic Potential')
        ax3.set_ylabel('Number of Cities')
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # Plot 4: Efficiency Metrics Dashboard
        ax4 = plt.subplot(2, 3, 4)
        
        metrics = results['efficiency_metrics']
        metric_names = ['Avg Path Length', 'Trade Volume\n(×10⁶)', 'Efficiency Score\n(×10⁻³)', 'Connectivity']
        metric_values = [
            metrics['average_path_length'],
            metrics['total_trade_volume'] / 1e6,
            metrics['efficiency_score'] / 1e3,
            metrics['connectivity'] * 100  # Convert to percentage
        ]
        
        bars = ax4.bar(metric_names, metric_values, color=['orange', 'green', 'purple', 'brown'])
        ax4.set_title('Network Efficiency Metrics', fontsize=14, fontweight='bold')
        ax4.set_ylabel('Metric Values')
        
        # Add value labels
        for bar, value in zip(bars, metric_values):
            ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + max(metric_values)*0.01,
                    f'{value:.2f}', ha='center', va='bottom', fontweight='bold')
        
        ax4.grid(True, alpha=0.3, axis='y')
        
        # Plot 5: City Rankings
        ax5 = plt.subplot(2, 3, 5)
        
        # Create ranking based on economic potential
        city_data = pd.DataFrame({
            'City': self.cities,
            'Population': self.populations,
            'GDP_per_capita': self.gdp_per_capita,
            'Economic_Potential': self.economic_potential,
            'Is_Hub': [i in results['hub_indices'] for i in range(self.n_cities)]
        })
        city_data = city_data.sort_values('Economic_Potential', ascending=True)
        
        colors = ['red' if is_hub else 'skyblue' for is_hub in city_data['Is_Hub']]
        bars = ax5.barh(city_data['City'], city_data['Economic_Potential'], color=colors)
        ax5.set_title('Cities Ranked by Economic Potential', fontsize=14, fontweight='bold')
        ax5.set_xlabel('Economic Potential')
        ax5.grid(True, alpha=0.3, axis='x')
        
        # Plot 6: Cost-Benefit Analysis
        ax6 = plt.subplot(2, 3, 6)
        
        # Calculate ROI for each hub
        hub_benefits = []
        hub_costs_individual = []
        
        for i, hub_idx in enumerate(results['hub_indices']):
            # Estimate benefit as economic potential of hub
            benefit = self.economic_potential[hub_idx]
            hub_benefits.append(benefit)
            
            # Estimate individual cost (proportional to connections)
            individual_cost = sum([mst[i][j]['weight'] for j in mst.neighbors(i)]) / 2
            hub_costs_individual.append(individual_cost)
        
        # Create scatter plot
        scatter = ax6.scatter(hub_costs_individual, hub_benefits, 
                            s=200, alpha=0.7, c=range(len(results['hub_indices'])), 
                            cmap='plasma', edgecolors='black')
        
        # Add hub labels
        for i, (cost, benefit) in enumerate(zip(hub_costs_individual, hub_benefits)):
            ax6.annotate(f"Hub {i+1}\n({results['hub_cities'][i]})", 
                        (cost, benefit), xytext=(10, 10), 
                        textcoords='offset points', fontsize=9,
                        bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.7))
        
        ax6.set_title('Hub Cost-Benefit Analysis', fontsize=14, fontweight='bold')
        ax6.set_xlabel('Infrastructure Cost (Million $)')
        ax6.set_ylabel('Economic Benefit')
        ax6.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        return fig

# Example usage and demonstration
def run_economic_zone_example():
    """
    Run a comprehensive example of economic zone optimization.
    """
    print("🌟 Economic Zone Development Optimization Example")
    print("=" * 60)
    
    # Define example cities with their characteristics
    cities = ['Tokyo', 'Osaka', 'Nagoya', 'Yokohama', 'Kobe', 'Kyoto', 'Fukuoka', 'Sendai']
    populations = [13929000, 2691000, 2295000, 3726000, 1538000, 1474000, 1581000, 1082000]  # Population
    gdp_per_capita = [48000, 42000, 45000, 46000, 43000, 41000, 38000, 40000]  # GDP per capita in USD
    
    # Initialize optimizer
    optimizer = EconomicZoneOptimizer(cities, populations, gdp_per_capita)
    
    # Run optimization
    results = optimizer.run_optimization(max_hubs=4)
    
    # Display results
    print(f"\n📊 OPTIMIZATION RESULTS")
    print(f"=" * 40)
    print(f"🏆 Selected Economic Hubs: {', '.join(results['hub_cities'])}")
    print(f"💰 Total Infrastructure Cost: ${results['total_cost']:.1f} Million")
    print(f"📈 Network Efficiency Score: {results['efficiency_metrics']['efficiency_score']:.2f}")
    print(f"🔗 Network Connectivity: {results['efficiency_metrics']['connectivity']:.1%}")
    print(f"📍 Average Path Length: {results['efficiency_metrics']['average_path_length']:.1f} units")
    print(f"💼 Total Trade Volume: ${results['efficiency_metrics']['total_trade_volume']:.0f}")
    
    # Create visualizations
    print(f"\n📊 Generating comprehensive visualizations...")
    fig = optimizer.visualize_results(results)
    
    # Additional analysis
    print(f"\n🔍 DETAILED ANALYSIS")
    print(f"=" * 40)
    
    # Hub efficiency analysis
    print("Hub Efficiency Rankings:")
    hub_data = []
    for i, hub_idx in enumerate(results['hub_indices']):
        efficiency = optimizer.economic_potential[hub_idx] / (optimizer.populations[hub_idx] / 1000)
        hub_data.append({
            'Hub': results['hub_cities'][i],
            'Population': optimizer.populations[hub_idx],
            'Economic_Potential': optimizer.economic_potential[hub_idx],
            'Efficiency_Ratio': efficiency
        })
    
    hub_df = pd.DataFrame(hub_data).sort_values('Efficiency_Ratio', ascending=False)
    for _, row in hub_df.iterrows():
        print(f"  🏙️ {row['Hub']}: Efficiency Ratio = {row['Efficiency_Ratio']:.2f}")
    
    return optimizer, results

# Run the example
if __name__ == "__main__":
    optimizer, results = run_economic_zone_example()