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

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, LinearConstraint
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd

# Set up the problem parameters
np.random.seed(42)  # For reproducible results

class SanctionsOptimizer:
    def __init__(self, n_sectors=5, n_neutral_countries=3):
        self.n_sectors = n_sectors
        self.n_neutral_countries = n_neutral_countries
        
        # Define sectors (example: Oil, Banking, Technology, Manufacturing, Agriculture)
        self.sectors = ['Oil & Energy', 'Banking', 'Technology', 'Manufacturing', 'Agriculture']
        
        # Economic weights of each sector in target country (normalized)
        self.weights = np.array([0.3, 0.25, 0.2, 0.15, 0.1])
        
        # Effectiveness coefficients (how much damage each sector creates)
        self.effectiveness = np.array([0.8, 0.9, 0.7, 0.6, 0.4])
        
        # Cost coefficients (cost to sanctioning country)
        self.costs = np.array([0.4, 0.2, 0.6, 0.3, 0.1])
        
        # Neutral country penalty coefficients
        self.neutral_penalties = np.array([0.3, 0.2, 0.1])
        
        # Trade interdependency matrix (how much neutral countries care about each sector)
        self.trade_matrix = np.random.rand(self.n_neutral_countries, self.n_sectors) * 0.5
        
        # Budget constraint (total sanctions capacity)
        self.budget = 3.0
    
    def objective_function(self, x):
        """
        Objective function to maximize:
        Benefits - Costs - Neutral country penalties
        
        We minimize the negative of this function
        """
        # Benefits: weighted effectiveness of sanctions
        benefits = np.sum(self.weights * x * self.effectiveness)
        
        # Costs: quadratic cost function (increasing marginal costs)
        costs = np.sum(self.costs * x**2)
        
        # Neutral country penalties (response to sanctions affecting their trade)
        neutral_responses = 0
        for j in range(self.n_neutral_countries):
            # Each neutral country responds based on trade exposure
            trade_exposure = np.sum(self.trade_matrix[j] * x)
            neutral_responses += self.neutral_penalties[j] * trade_exposure**1.5
        
        # Return negative (since we're minimizing)
        return -(benefits - costs - neutral_responses)
    
    def optimize_sanctions(self):
        """
        Solve the optimization problem using scipy
        """
        # Initial guess
        x0 = np.ones(self.n_sectors) * 0.3
        
        # Bounds: sanctions intensity between 0 and 1
        bounds = [(0, 1) for _ in range(self.n_sectors)]
        
        # Budget constraint: sum of x_i <= budget
        constraint = LinearConstraint(np.ones(self.n_sectors), 0, self.budget)
        
        # Solve optimization
        result = minimize(
            self.objective_function,
            x0,
            method='SLSQP',
            bounds=bounds,
            constraints=constraint,
            options={'disp': True}
        )
        
        return result
    
    def analyze_sensitivity(self, base_solution):
        """
        Perform sensitivity analysis on key parameters
        """
        # Test different budget levels
        budgets = np.linspace(1.0, 5.0, 20)
        optimal_values = []
        
        original_budget = self.budget
        
        for budget in budgets:
            self.budget = budget
            constraint = LinearConstraint(np.ones(self.n_sectors), 0, self.budget)
            
            result = minimize(
                self.objective_function,
                base_solution.x,
                method='SLSQP',
                bounds=[(0, 1) for _ in range(self.n_sectors)],
                constraints=constraint,
                options={'disp': False}
            )
            
            optimal_values.append(-result.fun)  # Convert back to positive
        
        self.budget = original_budget  # Restore original budget
        return budgets, optimal_values
    
    def monte_carlo_simulation(self, n_simulations=1000):
        """
        Monte Carlo simulation to account for uncertainty in parameters
        """
        results = []
        
        for _ in range(n_simulations):
            # Add random noise to parameters
            noise_factor = 0.1
            temp_weights = self.weights * (1 + np.random.normal(0, noise_factor, self.n_sectors))
            temp_effectiveness = self.effectiveness * (1 + np.random.normal(0, noise_factor, self.n_sectors))
            temp_costs = self.costs * (1 + np.random.normal(0, noise_factor, self.n_sectors))
            
            # Store original values
            orig_weights, orig_effectiveness, orig_costs = self.weights, self.effectiveness, self.costs
            
            # Update with noisy values
            self.weights, self.effectiveness, self.costs = temp_weights, temp_effectiveness, temp_costs
            
            # Solve optimization
            result = self.optimize_sanctions()
            if result.success:
                results.append(result.x)
            
            # Restore original values
            self.weights, self.effectiveness, self.costs = orig_weights, orig_effectiveness, orig_costs
        
        return np.array(results)

# Initialize and solve the optimization problem
print("=== Economic Sanctions Strategy Optimization ===\n")
optimizer = SanctionsOptimizer()

print("Problem Setup:")
print(f"Sectors: {optimizer.sectors}")
print(f"Economic weights: {optimizer.weights}")
print(f"Effectiveness coefficients: {optimizer.effectiveness}")
print(f"Cost coefficients: {optimizer.costs}")
print(f"Budget constraint: {optimizer.budget}")
print("\n")

# Solve the main optimization problem
print("Solving optimization problem...")
optimal_result = optimizer.optimize_sanctions()

print(f"\nOptimization Result:")
print(f"Success: {optimal_result.success}")
print(f"Optimal objective value: {-optimal_result.fun:.4f}")
print(f"Optimal sanctions allocation:")

for i, (sector, allocation) in enumerate(zip(optimizer.sectors, optimal_result.x)):
    print(f"  {sector}: {allocation:.3f}")

# Perform sensitivity analysis
print("\nPerforming sensitivity analysis...")
budgets, optimal_values = optimizer.analyze_sensitivity(optimal_result)

# Monte Carlo simulation
print("Running Monte Carlo simulation...")
mc_results = optimizer.monte_carlo_simulation(500)

print(f"Monte Carlo completed with {len(mc_results)} successful runs")

# Create comprehensive visualizations
fig = plt.figure(figsize=(20, 15))

# Plot 1: Optimal allocation by sector
ax1 = plt.subplot(3, 3, 1)
bars = plt.bar(optimizer.sectors, optimal_result.x, 
               color=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'])
plt.title('Optimal Sanctions Allocation by Sector', fontsize=14, fontweight='bold')
plt.ylabel('Sanctions Intensity')
plt.xticks(rotation=45, ha='right')
plt.ylim(0, 1)

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

# Plot 2: Economic impact breakdown
ax2 = plt.subplot(3, 3, 2)
impact_data = optimizer.weights * optimal_result.x * optimizer.effectiveness
plt.pie(impact_data, labels=optimizer.sectors, autopct='%1.1f%%', 
        colors=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'])
plt.title('Economic Impact Distribution', fontsize=14, fontweight='bold')

# Plot 3: Cost-benefit analysis
ax3 = plt.subplot(3, 3, 3)
benefits = optimizer.weights * optimal_result.x * optimizer.effectiveness
costs = optimizer.costs * optimal_result.x**2
x_pos = np.arange(len(optimizer.sectors))

plt.bar(x_pos - 0.2, benefits, 0.4, label='Benefits', alpha=0.8, color='green')
plt.bar(x_pos + 0.2, costs, 0.4, label='Costs', alpha=0.8, color='red')
plt.title('Cost-Benefit Analysis by Sector', fontsize=14, fontweight='bold')
plt.xlabel('Sectors')
plt.ylabel('Value')
plt.xticks(x_pos, optimizer.sectors, rotation=45, ha='right')
plt.legend()

# Plot 4: Sensitivity analysis
ax4 = plt.subplot(3, 3, 4)
plt.plot(budgets, optimal_values, 'b-', linewidth=2, marker='o', markersize=4)
plt.axvline(x=optimizer.budget, color='r', linestyle='--', alpha=0.7, 
           label=f'Current Budget ({optimizer.budget})')
plt.title('Sensitivity Analysis: Budget vs Optimal Value', fontsize=14, fontweight='bold')
plt.xlabel('Budget Constraint')
plt.ylabel('Optimal Objective Value')
plt.grid(True, alpha=0.3)
plt.legend()

# Plot 5: Monte Carlo results - box plots
ax5 = plt.subplot(3, 3, 5)
mc_df = pd.DataFrame(mc_results, columns=optimizer.sectors)
plt.boxplot([mc_df[sector] for sector in optimizer.sectors], 
            labels=optimizer.sectors, patch_artist=True)
plt.title('Monte Carlo Simulation Results', fontsize=14, fontweight='bold')
plt.ylabel('Sanctions Intensity')
plt.xticks(rotation=45, ha='right')

# Add mean lines
for i, sector in enumerate(optimizer.sectors):
    mean_val = mc_df[sector].mean()
    plt.axhline(y=mean_val, xmin=(i+0.6)/len(optimizer.sectors), 
                xmax=(i+1.4)/len(optimizer.sectors), color='red', alpha=0.7)

# Plot 6: Trade interdependency heatmap
ax6 = plt.subplot(3, 3, 6)
neutral_countries = [f'Country {i+1}' for i in range(optimizer.n_neutral_countries)]
sns.heatmap(optimizer.trade_matrix, 
            xticklabels=optimizer.sectors, 
            yticklabels=neutral_countries,
            annot=True, fmt='.2f', cmap='YlOrRd')
plt.title('Trade Interdependency Matrix', fontsize=14, fontweight='bold')
plt.xlabel('Sectors')
plt.ylabel('Neutral Countries')

# Plot 7: Risk-return scatter
ax7 = plt.subplot(3, 3, 7)
expected_return = optimizer.weights * optimizer.effectiveness
risk = optimizer.costs + np.sum(optimizer.trade_matrix, axis=0) * optimizer.neutral_penalties.mean()

scatter = plt.scatter(risk, expected_return, s=optimal_result.x*500, 
                     c=range(len(optimizer.sectors)), cmap='viridis', alpha=0.7)
plt.title('Risk-Return Analysis', fontsize=14, fontweight='bold')
plt.xlabel('Risk (Costs + Neutral Response)')
plt.ylabel('Expected Return (Effectiveness)')

# Add sector labels
for i, sector in enumerate(optimizer.sectors):
    plt.annotate(sector, (risk[i], expected_return[i]), 
                xytext=(5, 5), textcoords='offset points', fontsize=9)

# Plot 8: Cumulative impact over time simulation
ax8 = plt.subplot(3, 3, 8)
time_periods = np.arange(1, 25)
# Simulate diminishing effectiveness over time
effectiveness_decay = np.exp(-time_periods * 0.1)
cumulative_impact = np.cumsum(
    np.sum(optimizer.weights * optimal_result.x * optimizer.effectiveness) * effectiveness_decay
)

plt.plot(time_periods, cumulative_impact, 'g-', linewidth=3, marker='s', markersize=4)
plt.title('Projected Cumulative Economic Impact', fontsize=14, fontweight='bold')
plt.xlabel('Time Periods (quarters)')
plt.ylabel('Cumulative Impact')
plt.grid(True, alpha=0.3)

# Plot 9: Strategy comparison
ax9 = plt.subplot(3, 3, 9)
# Compare with uniform allocation strategy
uniform_allocation = np.ones(optimizer.n_sectors) * (optimizer.budget / optimizer.n_sectors)
uniform_objective = -optimizer.objective_function(uniform_allocation)
optimal_objective = -optimal_result.fun

strategies = ['Uniform\nAllocation', 'Optimized\nAllocation']
objectives = [uniform_objective, optimal_objective]
colors = ['lightcoral', 'lightgreen']

bars = plt.bar(strategies, objectives, color=colors, alpha=0.8, edgecolor='black')
plt.title('Strategy Comparison', fontsize=14, fontweight='bold')
plt.ylabel('Objective Value')

# Add improvement percentage
improvement = ((optimal_objective - uniform_objective) / uniform_objective) * 100
plt.text(1, optimal_objective + 0.01, f'+{improvement:.1f}%', 
         ha='center', va='bottom', fontweight='bold', fontsize=12, color='darkgreen')

for bar, value in zip(bars, objectives):
    plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005, 
             f'{value:.3f}', ha='center', va='bottom', fontweight='bold')

plt.tight_layout()
plt.show()

# Print detailed analysis
print("\n=== DETAILED ANALYSIS ===")
print(f"1. OPTIMAL STRATEGY:")
print(f"   Total sanctions intensity: {np.sum(optimal_result.x):.3f} (Budget: {optimizer.budget})")
print(f"   Budget utilization: {(np.sum(optimal_result.x)/optimizer.budget)*100:.1f}%")

print(f"\n2. SECTOR PRIORITIES:")
sector_priorities = list(zip(optimizer.sectors, optimal_result.x))
sector_priorities.sort(key=lambda x: x[1], reverse=True)
for i, (sector, intensity) in enumerate(sector_priorities):
    print(f"   {i+1}. {sector}: {intensity:.3f}")

print(f"\n3. EXPECTED OUTCOMES:")
total_benefit = np.sum(optimizer.weights * optimal_result.x * optimizer.effectiveness)
total_cost = np.sum(optimizer.costs * optimal_result.x**2)
neutral_penalty = sum([optimizer.neutral_penalties[j] * 
                      np.sum(optimizer.trade_matrix[j] * optimal_result.x)**1.5 
                      for j in range(optimizer.n_neutral_countries)])

print(f"   Expected economic impact on target: {total_benefit:.4f}")
print(f"   Cost to sanctioning country: {total_cost:.4f}")
print(f"   Neutral country backlash: {neutral_penalty:.4f}")
print(f"   Net strategic value: {total_benefit - total_cost - neutral_penalty:.4f}")

print(f"\n4. MONTE CARLO INSIGHTS:")
print(f"   Average allocation ranges:")
for i, sector in enumerate(optimizer.sectors):
    mean_val = np.mean(mc_results[:, i])
    std_val = np.std(mc_results[:, i])
    print(f"   {sector}: {mean_val:.3f} ± {std_val:.3f}")

print(f"\n5. STRATEGIC RECOMMENDATIONS:")
print(f"   - Focus on {sector_priorities[0][0]} and {sector_priorities[1][0]} sectors")

# Find closest budget value to current budget for comparison
current_budget_idx = np.argmin(np.abs(budgets - optimizer.budget))
current_value = optimal_values[current_budget_idx]
max_value = max(optimal_values)
optimal_budget = budgets[np.argmax(optimal_values)]

if max_value > current_value:
    improvement = ((max_value - current_value) / current_value * 100)
    print(f"   - Budget increase to {optimal_budget:.1f} could improve effectiveness by {improvement:.1f}%")
else:
    print(f"   - Current budget of {optimizer.budget:.1f} is near optimal")
print(f"   - Consider diplomatic engagement with neutral countries to minimize backlash")