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()

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

# Set style for better visualization
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

class SupplyChainOptimizer:
    def __init__(self):
        """
        Initialize the Supply Chain Geopolitical Optimizer
        """
        # Define regions and their geopolitical risk scores (1-10, higher = more risky)
        self.suppliers = {
            'China': {'capacity': 1000, 'cost_base': 1.0, 'risk': 6.5},
            'Vietnam': {'capacity': 600, 'cost_base': 1.2, 'risk': 4.0},
            'Taiwan': {'capacity': 800, 'cost_base': 1.1, 'risk': 7.0},
            'South Korea': {'capacity': 700, 'cost_base': 1.3, 'risk': 3.5},
            'India': {'capacity': 500, 'cost_base': 1.1, 'risk': 4.5}
        }
        
        self.customers = {
            'USA': {'demand': 800, 'region': 'Americas'},
            'Germany': {'demand': 600, 'region': 'Europe'},
            'Japan': {'demand': 500, 'region': 'Asia'},
            'UK': {'demand': 400, 'region': 'Europe'},
            'Brazil': {'demand': 300, 'region': 'Americas'}
        }
        
        self.facilities = {
            'Singapore Hub': {'capacity': 1200, 'region': 'Asia', 'cost_mult': 1.0},
            'Netherlands Hub': {'capacity': 800, 'region': 'Europe', 'cost_mult': 1.1},
            'USA Hub': {'capacity': 1000, 'region': 'Americas', 'cost_mult': 1.05}
        }
        
        # Distance matrix for transportation costs (normalized)
        self.distances = self._create_distance_matrix()
        
        # Weight parameters for multi-objective optimization
        self.alpha = 0.4  # Cost weight
        self.beta = 0.4   # Risk weight
        self.gamma = 0.2  # Diversification weight
        
    def _create_distance_matrix(self):
        """Create a normalized distance matrix for transportation costs"""
        # Simplified distance matrix (in practice, use actual distances)
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        # Create distance matrices
        sup_to_fac = np.random.uniform(0.8, 2.0, (len(suppliers), len(facilities)))
        fac_to_cust = np.random.uniform(0.5, 1.5, (len(facilities), len(customers)))
        
        # Set realistic distances based on geography
        # Asia suppliers to Singapore Hub should be cheaper
        sup_to_fac[0:3, 0] *= 0.6  # China, Vietnam, Taiwan to Singapore
        sup_to_fac[3, 0] *= 0.7    # South Korea to Singapore
        
        return {'sup_to_fac': sup_to_fac, 'fac_to_cust': fac_to_cust}
    
    def calculate_total_cost(self, solution):
        """
        Calculate total logistics cost for a given solution
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        total_cost = 0
        idx = 0
        
        for i, supplier in enumerate(suppliers):
            for j, customer in enumerate(customers):
                for k, facility in enumerate(facilities):
                    if idx < len(solution):
                        flow = solution[idx]
                        
                        # Calculate transportation cost
                        transport_cost = (self.distances['sup_to_fac'][i, k] + 
                                        self.distances['fac_to_cust'][k, j])
                        
                        # Add base production cost
                        production_cost = self.suppliers[supplier]['cost_base']
                        
                        # Add facility cost multiplier
                        facility_cost = self.facilities[facility]['cost_mult']
                        
                        total_cost += flow * transport_cost * production_cost * facility_cost
                        idx += 1
        
        return total_cost
    
    def calculate_risk_score(self, solution):
        """
        Calculate total geopolitical risk score
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        total_risk = 0
        idx = 0
        
        for i, supplier in enumerate(suppliers):
            for j, customer in enumerate(customers):
                for k, facility in enumerate(facilities):
                    if idx < len(solution):
                        flow = solution[idx]
                        risk = self.suppliers[supplier]['risk']
                        
                        # Add facility region risk (simplified)
                        facility_region_risk = {'Asia': 1.2, 'Europe': 0.8, 'Americas': 1.0}
                        fac_risk = facility_region_risk[self.facilities[facility]['region']]
                        
                        total_risk += flow * risk * fac_risk
                        idx += 1
        
        return total_risk
    
    def calculate_diversification(self, solution):
        """
        Calculate supply diversification score (higher = more diverse)
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        # Calculate Herfindahl-Hirschman Index for each customer
        total_diversification = 0
        idx = 0
        
        for j, customer in enumerate(customers):
            supplier_shares = np.zeros(len(suppliers))
            total_supply = 0
            
            for i, supplier in enumerate(suppliers):
                for k, facility in enumerate(facilities):
                    if idx < len(solution):
                        flow = solution[idx]
                        supplier_shares[i] += flow
                        total_supply += flow
                        idx += 1
            
            if total_supply > 0:
                supplier_shares /= total_supply
                # Calculate diversity index (1 - HHI)
                hhi = np.sum(supplier_shares ** 2)
                diversity = 1 - hhi
                total_diversification += diversity * self.customers[customer]['demand']
        
        return total_diversification
    
    def objective_function(self, solution):
        """
        Multi-objective function combining cost, risk, and diversification
        """
        cost = self.calculate_total_cost(solution)
        risk = self.calculate_risk_score(solution)
        diversification = self.calculate_diversification(solution)
        
        # Normalize components (simplified normalization)
        normalized_cost = cost / 10000
        normalized_risk = risk / 10000
        normalized_diversification = diversification / 1000
        
        objective = (self.alpha * normalized_cost + 
                    self.beta * normalized_risk - 
                    self.gamma * normalized_diversification)
        
        return objective
    
    def create_constraints(self):
        """
        Create constraint functions for the optimization
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        constraints = []
        
        # Demand satisfaction constraints
        for j, customer in enumerate(customers):
            def demand_constraint(solution, customer_idx=j):
                total_received = 0
                idx = 0
                for i in range(len(suppliers)):
                    for cust_idx in range(len(customers)):
                        for k in range(len(facilities)):
                            if cust_idx == customer_idx and idx < len(solution):
                                total_received += solution[idx]
                            idx += 1
                return self.customers[list(customers)[customer_idx]]['demand'] - total_received
            
            constraints.append({'type': 'eq', 'fun': demand_constraint})
        
        # Supply capacity constraints
        for i, supplier in enumerate(suppliers):
            def supply_constraint(solution, supplier_idx=i):
                total_supplied = 0
                idx = 0
                for sup_idx in range(len(suppliers)):
                    for j in range(len(customers)):
                        for k in range(len(facilities)):
                            if sup_idx == supplier_idx and idx < len(solution):
                                total_supplied += solution[idx]
                            idx += 1
                return self.suppliers[list(suppliers)[supplier_idx]]['capacity'] - total_supplied
            
            constraints.append({'type': 'ineq', 'fun': supply_constraint})
        
        return constraints
    
    def optimize(self):
        """
        Perform the optimization
        """
        n_vars = len(self.suppliers) * len(self.customers) * len(self.facilities)
        
        # Initial guess
        x0 = np.random.uniform(0, 100, n_vars)
        
        # Bounds (non-negative flows)
        bounds = [(0, None) for _ in range(n_vars)]
        
        # Constraints
        constraints = self.create_constraints()
        
        # Solve optimization
        result = minimize(
            self.objective_function,
            x0,
            method='SLSQP',
            bounds=bounds,
            constraints=constraints,
            options={'maxiter': 1000, 'disp': True}
        )
        
        return result
    
    def analyze_solution(self, solution):
        """
        Analyze and visualize the optimal solution
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        # Create solution matrix
        flow_matrix = np.zeros((len(suppliers), len(customers), len(facilities)))
        idx = 0
        
        for i in range(len(suppliers)):
            for j in range(len(customers)):
                for k in range(len(facilities)):
                    if idx < len(solution):
                        flow_matrix[i, j, k] = max(0, solution[idx])
                        idx += 1
        
        # Calculate metrics
        total_cost = self.calculate_total_cost(solution)
        risk_score = self.calculate_risk_score(solution)
        diversification = self.calculate_diversification(solution)
        
        print(f"Optimization Results:")
        print(f"Total Cost: {total_cost:.2f}")
        print(f"Risk Score: {risk_score:.2f}")
        print(f"Diversification Index: {diversification:.2f}")
        print(f"Combined Objective: {self.objective_function(solution):.2f}")
        
        return flow_matrix, {'cost': total_cost, 'risk': risk_score, 'diversification': diversification}
    
    def visualize_results(self, flow_matrix, metrics):
        """
        Create comprehensive visualizations of the results
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        fig, axes = plt.subplots(2, 3, figsize=(20, 12))
        fig.suptitle('Supply Chain Geopolitical Optimization Results', fontsize=16, fontweight='bold')
        
        # 1. Supplier utilization
        supplier_utilization = np.sum(flow_matrix, axis=(1, 2))
        supplier_capacity = [self.suppliers[s]['capacity'] for s in suppliers]
        utilization_rate = supplier_utilization / supplier_capacity * 100
        
        axes[0, 0].bar(suppliers, utilization_rate, color='skyblue', edgecolor='navy')
        axes[0, 0].set_title('Supplier Utilization Rate (%)')
        axes[0, 0].set_ylabel('Utilization (%)')
        axes[0, 0].tick_params(axis='x', rotation=45)
        
        # 2. Customer supply sources
        customer_supply = np.sum(flow_matrix, axis=2).T
        bottom = np.zeros(len(customers))
        colors = plt.cm.Set3(np.linspace(0, 1, len(suppliers)))
        
        for i, supplier in enumerate(suppliers):
            axes[0, 1].bar(customers, customer_supply[:, i], bottom=bottom, 
                          label=supplier, color=colors[i])
            bottom += customer_supply[:, i]
        
        axes[0, 1].set_title('Customer Supply Sources')
        axes[0, 1].set_ylabel('Supply Volume')
        axes[0, 1].legend(bbox_to_anchor=(1.05, 1), loc='upper left')
        axes[0, 1].tick_params(axis='x', rotation=45)
        
        # 3. Facility throughput
        facility_throughput = np.sum(flow_matrix, axis=(0, 1))
        facility_capacity = [self.facilities[f]['capacity'] for f in facilities]
        
        x = np.arange(len(facilities))
        width = 0.35
        axes[0, 2].bar(x - width/2, facility_throughput, width, label='Throughput', color='lightgreen')
        axes[0, 2].bar(x + width/2, facility_capacity, width, label='Capacity', color='lightcoral')
        axes[0, 2].set_title('Facility Throughput vs Capacity')
        axes[0, 2].set_ylabel('Volume')
        axes[0, 2].set_xticks(x)
        axes[0, 2].set_xticklabels(facilities, rotation=45)
        axes[0, 2].legend()
        
        # 4. Risk vs Cost scatter
        supplier_costs = []
        supplier_risks = []
        supplier_volumes = []
        
        for i, supplier in enumerate(suppliers):
            volume = np.sum(flow_matrix[i, :, :])
            cost = self.suppliers[supplier]['cost_base'] * volume
            risk = self.suppliers[supplier]['risk'] * volume
            
            supplier_costs.append(cost)
            supplier_risks.append(risk)
            supplier_volumes.append(volume)
        
        scatter = axes[1, 0].scatter(supplier_costs, supplier_risks, 
                                   s=[v/5 for v in supplier_volumes], 
                                   c=range(len(suppliers)), cmap='viridis', alpha=0.7)
        axes[1, 0].set_title('Risk vs Cost by Supplier')
        axes[1, 0].set_xlabel('Cost Impact')
        axes[1, 0].set_ylabel('Risk Impact')
        
        # Add supplier labels
        for i, supplier in enumerate(suppliers):
            axes[1, 0].annotate(supplier, (supplier_costs[i], supplier_risks[i]),
                              xytext=(5, 5), textcoords='offset points', fontsize=8)
        
        # 5. Geographic distribution
        facility_flows = {}
        for k, facility in enumerate(facilities):
            facility_flows[facility] = np.sum(flow_matrix[:, :, k])
        
        axes[1, 1].pie(facility_flows.values(), labels=facility_flows.keys(), autopct='%1.1f%%')
        axes[1, 1].set_title('Geographic Distribution of Flows')
        
        # 6. Metrics comparison
        metrics_names = ['Cost\n(normalized)', 'Risk\n(normalized)', 'Diversification\n(normalized)']
        metrics_values = [
            metrics['cost'] / 10000,
            metrics['risk'] / 10000,
            metrics['diversification'] / 1000
        ]
        
        colors_metrics = ['red', 'orange', 'green']
        bars = axes[1, 2].bar(metrics_names, metrics_values, color=colors_metrics, alpha=0.7)
        axes[1, 2].set_title('Optimization Metrics')
        axes[1, 2].set_ylabel('Normalized Value')
        
        # Add value labels on bars
        for bar, value in zip(bars, metrics_values):
            height = bar.get_height()
            axes[1, 2].text(bar.get_x() + bar.get_width()/2., height + 0.01,
                           f'{value:.3f}', ha='center', va='bottom')
        
        plt.tight_layout()
        plt.show()
        
        # Additional heatmap for flow matrix
        self.create_flow_heatmap(flow_matrix)
    
    def create_flow_heatmap(self, flow_matrix):
        """
        Create a detailed heatmap of flows
        """
        suppliers = list(self.suppliers.keys())
        facilities = list(self.facilities.keys())
        customers = list(self.customers.keys())
        
        fig, axes = plt.subplots(1, len(facilities), figsize=(18, 6))
        fig.suptitle('Flow Distribution Heatmaps by Facility', fontsize=14, fontweight='bold')
        
        for k, facility in enumerate(facilities):
            flow_data = flow_matrix[:, :, k]
            
            im = axes[k].imshow(flow_data, cmap='YlOrRd', aspect='auto')
            axes[k].set_title(f'{facility}')
            axes[k].set_xlabel('Customers')
            axes[k].set_ylabel('Suppliers')
            axes[k].set_xticks(range(len(customers)))
            axes[k].set_yticks(range(len(suppliers)))
            axes[k].set_xticklabels(customers, rotation=45)
            axes[k].set_yticklabels(suppliers)
            
            # Add text annotations
            for i in range(len(suppliers)):
                for j in range(len(customers)):
                    text = axes[k].text(j, i, f'{flow_data[i, j]:.0f}',
                                      ha="center", va="center", color="black", fontsize=8)
            
            plt.colorbar(im, ax=axes[k], shrink=0.8)
        
        plt.tight_layout()
        plt.show()

# Run the optimization
def main():
    print("Supply Chain Geopolitical Optimization")
    print("="*50)
    
    # Initialize optimizer
    optimizer = SupplyChainOptimizer()
    
    print("\nProblem Setup:")
    print(f"Suppliers: {list(optimizer.suppliers.keys())}")
    print(f"Customers: {list(optimizer.customers.keys())}")
    print(f"Facilities: {list(optimizer.facilities.keys())}")
    print(f"Optimization weights - Cost: {optimizer.alpha}, Risk: {optimizer.beta}, Diversification: {optimizer.gamma}")
    
    # Perform optimization
    print("\nPerforming optimization...")
    result = optimizer.optimize()
    
    if result.success:
        print(f"\nOptimization successful!")
        print(f"Iterations: {result.nit}")
        print(f"Function evaluations: {result.nfev}")
        
        # Analyze results
        flow_matrix, metrics = optimizer.analyze_solution(result.x)
        
        # Create visualizations
        print("\nGenerating visualizations...")
        optimizer.visualize_results(flow_matrix, metrics)
        
    else:
        print(f"\nOptimization failed: {result.message}")
        
    return optimizer, result

# Execute the main function
if __name__ == "__main__":
    optimizer, result = main()