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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

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

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

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

# Set up the naval strategy optimization problem
print("=== Naval Fleet Deployment Optimization ===\n")

# Define ship types and zones
ship_types = ['Destroyer', 'Cruiser', 'Aircraft Carrier']
zones = ['North Pacific', 'South Pacific', 'Atlantic']

# Strategic value matrix v_ij (value of deploying ship type i to zone j)
# Higher values indicate better strategic positioning
strategic_value = np.array([
    [8.5, 7.2, 6.8],  # Destroyer values for each zone
    [9.1, 8.6, 8.9],  # Cruiser values for each zone
    [12.5, 11.8, 13.2]  # Aircraft Carrier values for each zone
])

# Cost matrix c_ij (operational cost of deploying ship type i to zone j)
# Includes fuel, maintenance, and operational expenses (in millions)
operational_cost = np.array([
    [2.1, 2.8, 3.2],  # Destroyer costs for each zone
    [4.5, 5.1, 5.8],  # Cruiser costs for each zone
    [8.2, 9.1, 10.5]  # Aircraft Carrier costs for each zone
])

# Available ships of each type
available_ships = np.array([15, 8, 4])  # [Destroyers, Cruisers, Aircraft Carriers]

# Total budget (in millions)
total_budget = 180

print("Strategic Value Matrix (higher is better):")
value_df = pd.DataFrame(strategic_value, index=ship_types, columns=zones)
print(value_df)
print(f"\nOperational Cost Matrix (millions USD):")
cost_df = pd.DataFrame(operational_cost, index=ship_types, columns=zones)
print(cost_df)
print(f"\nAvailable Ships: {dict(zip(ship_types, available_ships))}")
print(f"Total Budget: ${total_budget} million")

# Convert to linear programming format
# We need to flatten the 2D decision variables into 1D
# x = [x_00, x_01, x_02, x_10, x_11, x_12, x_20, x_21, x_22]
# where x_ij represents ships of type i deployed to zone j

# Objective function coefficients (negative because linprog minimizes)
c = -strategic_value.flatten()

# Inequality constraint matrix A_ub and bounds b_ub
# Ship availability constraints: sum over zones for each ship type
A_ub = []
b_ub = []

# Ship availability constraints
for i in range(3):  # for each ship type
    constraint = np.zeros(9)
    for j in range(3):  # for each zone
        constraint[i*3 + j] = 1
    A_ub.append(constraint)
    b_ub.append(available_ships[i])

# Budget constraint
budget_constraint = operational_cost.flatten()
A_ub.append(budget_constraint)
b_ub.append(total_budget)

A_ub = np.array(A_ub)
b_ub = np.array(b_ub)

# Bounds for variables (non-negativity)
bounds = [(0, None) for _ in range(9)]

print(f"\n=== Solving Optimization Problem ===")
print("Constraints:")
for i, ship_type in enumerate(ship_types):
    print(f"- {ship_type} availability: ≤ {available_ships[i]} ships")
print(f"- Total budget: ≤ ${total_budget} million")

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

if result.success:
    print(f"\nOptimization Status: SUCCESS")
    print(f"Maximum Strategic Value: {-result.fun:.2f}")
    
    # Reshape solution back to 2D matrix
    optimal_deployment = result.x.reshape(3, 3)
    
    print(f"\nOptimal Fleet Deployment:")
    deployment_df = pd.DataFrame(optimal_deployment, index=ship_types, columns=zones)
    print(deployment_df.round(2))
    
    # Calculate total cost and resource usage
    total_cost = np.sum(operational_cost * optimal_deployment)
    total_ships_used = np.sum(optimal_deployment, axis=1)
    
    print(f"\nResource Utilization:")
    for i, ship_type in enumerate(ship_types):
        print(f"- {ship_type}: {total_ships_used[i]:.1f}/{available_ships[i]} ships ({100*total_ships_used[i]/available_ships[i]:.1f}%)")
    
    print(f"\nBudget Utilization: ${total_cost:.1f}/{total_budget} million ({100*total_cost/total_budget:.1f}%)")
    
else:
    print("Optimization failed!")
    print(result.message)

# Visualizations
plt.style.use('default')
fig = plt.figure(figsize=(18, 12))

# 1. Strategic Value Heatmap
plt.subplot(2, 3, 1)
sns.heatmap(strategic_value, annot=True, fmt='.1f', cmap='YlOrRd', 
            xticklabels=zones, yticklabels=ship_types, cbar_kws={'label': 'Strategic Value'})
plt.title('Strategic Value Matrix\n(Higher values = Better positioning)', fontsize=12, fontweight='bold')
plt.ylabel('Ship Types')

# 2. Operational Cost Heatmap
plt.subplot(2, 3, 2)
sns.heatmap(operational_cost, annot=True, fmt='.1f', cmap='Blues', 
            xticklabels=zones, yticklabels=ship_types, cbar_kws={'label': 'Cost (Millions USD)'})
plt.title('Operational Cost Matrix\n(Cost per ship deployment)', fontsize=12, fontweight='bold')
plt.ylabel('Ship Types')

# 3. Optimal Deployment Heatmap
plt.subplot(2, 3, 3)
sns.heatmap(optimal_deployment, annot=True, fmt='.1f', cmap='Greens', 
            xticklabels=zones, yticklabels=ship_types, cbar_kws={'label': 'Number of Ships'})
plt.title('Optimal Fleet Deployment\n(Number of ships per zone)', fontsize=12, fontweight='bold')
plt.ylabel('Ship Types')

# 4. Ship Utilization Bar Chart
plt.subplot(2, 3, 4)
utilization_pct = 100 * total_ships_used / available_ships
bars = plt.bar(ship_types, utilization_pct, color=['skyblue', 'lightcoral', 'lightgreen'])
plt.ylabel('Utilization (%)')
plt.title('Ship Type Utilization', fontsize=12, fontweight='bold')
plt.xticks(rotation=45)
for i, bar in enumerate(bars):
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height + 1,
             f'{height:.1f}%\n({total_ships_used[i]:.1f}/{available_ships[i]})',
             ha='center', va='bottom')
plt.ylim(0, 105)

# 5. Zone Deployment Distribution
plt.subplot(2, 3, 5)
zone_totals = np.sum(optimal_deployment, axis=0)
colors = ['#FF9999', '#66B2FF', '#99FF99']
bars = plt.bar(zones, zone_totals, color=colors)
plt.ylabel('Total Ships Deployed')
plt.title('Ships per Strategic Zone', fontsize=12, fontweight='bold')
plt.xticks(rotation=45)
for bar, total in zip(bars, zone_totals):
    plt.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.1,
             f'{total:.1f}', ha='center', va='bottom')

# 6. Value vs Cost Efficiency Analysis
plt.subplot(2, 3, 6)
efficiency = strategic_value / operational_cost
zone_efficiency = np.mean(efficiency, axis=0)
ship_efficiency = np.mean(efficiency, axis=1)

x_pos = np.arange(len(zones))
width = 0.35
bars1 = plt.bar(x_pos - width/2, zone_efficiency, width, label='Zone Efficiency', color='lightblue')
bars2 = plt.bar(x_pos + width/2, ship_efficiency, width, label='Ship Efficiency', color='orange')

plt.ylabel('Value/Cost Ratio')
plt.title('Strategic Efficiency Analysis', fontsize=12, fontweight='bold')
plt.xticks(x_pos, zones)
plt.legend()

for bars in [bars1, bars2]:
    for bar in bars:
        height = bar.get_height()
        plt.text(bar.get_x() + bar.get_width()/2., height + 0.02,
                 f'{height:.2f}', ha='center', va='bottom', fontsize=8)

plt.tight_layout()
plt.show()

# Detailed Analysis
print(f"\n=== Strategic Analysis ===")
print(f"\n1. DEPLOYMENT SUMMARY:")
for i, ship_type in enumerate(ship_types):
    for j, zone in enumerate(zones):
        if optimal_deployment[i, j] > 0.1:  # Only show significant deployments
            value_contrib = optimal_deployment[i, j] * strategic_value[i, j]
            cost_contrib = optimal_deployment[i, j] * operational_cost[i, j]
            print(f"   • {optimal_deployment[i, j]:.1f} {ship_type}(s) → {zone}")
            print(f"     Strategic Value: {value_contrib:.1f}, Cost: ${cost_contrib:.1f}M")

print(f"\n2. EFFICIENCY METRICS:")
total_strategic_value = -result.fun
roi = total_strategic_value / total_cost
print(f"   • Strategic ROI: {roi:.2f} value points per million USD")
print(f"   • Average cost per strategic point: ${total_cost/total_strategic_value:.2f}M")

print(f"\n3. STRATEGIC RECOMMENDATIONS:")
# Find the most and least efficient zones
zone_values = np.sum(optimal_deployment * strategic_value, axis=0)
zone_costs = np.sum(optimal_deployment * operational_cost, axis=0)
zone_roi = zone_values / (zone_costs + 0.001)  # Add small value to avoid division by zero

best_zone = zones[np.argmax(zone_roi)]
worst_zone = zones[np.argmin(zone_roi)]

print(f"   • Highest ROI Zone: {best_zone} ({zone_roi[np.argmax(zone_roi)]:.2f})")
print(f"   • Focus Area: Consider increasing presence in high-value, low-cost zones")
print(f"   • Resource Optimization: {100*total_cost/total_budget:.1f}% budget utilization indicates {'efficient' if total_cost/total_budget > 0.9 else 'potential for expansion'}")

# Sensitivity Analysis Preview
print(f"\n4. SENSITIVITY INSIGHTS:")
print(f"   • Critical Resource: {ship_types[np.argmax(total_ships_used/available_ships)]} (highest utilization)")
print(f"   • Strategic Priority: {zones[np.argmax(zone_totals)]} (most ships deployed)")
print(f"   • Cost Driver: ${np.max(operational_cost):.1f}M max single deployment cost")

strategic_value = np.array([
    [8.5, 7.2, 6.8],  # Destroyer effectiveness per zone
    [9.1, 8.6, 8.9],  # Cruiser effectiveness per zone  
    [12.5, 11.8, 13.2]  # Aircraft Carrier effectiveness per zone
])

c = -strategic_value.flatten()  # Negative because linprog minimizes

for i in range(3):  # for each ship type
    constraint = np.zeros(9)
    for j in range(3):  # for each zone
        constraint[i*3 + j] = 1  # Sum across zones for ship type i
    A_ub.append(constraint)
    b_ub.append(available_ships[i])

result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

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

# Set up the cultural programs and their parameters
programs = {
    'Cultural Exchanges': {'efficiency': 2.5, 'min_investment': 0.5},
    'Language Education': {'efficiency': 3.0, 'min_investment': 1.0}, 
    'Arts & Film Festivals': {'efficiency': 2.2, 'min_investment': 0.3},
    'Digital Content Creation': {'efficiency': 3.5, 'min_investment': 0.8},
    'Academic Partnerships': {'efficiency': 2.8, 'min_investment': 1.2},
    'Sports Diplomacy': {'efficiency': 2.0, 'min_investment': 0.4},
    'Tourism Promotion': {'efficiency': 1.8, 'min_investment': 0.6},
    'Technology Showcases': {'efficiency': 3.2, 'min_investment': 0.9}
}

# Parameters
total_budget = 10.0  # $10 million
n_programs = len(programs)
program_names = list(programs.keys())

# Extract efficiency rates and minimum investments
efficiency_rates = [programs[name]['efficiency'] for name in program_names]
min_investments = [programs[name]['min_investment'] for name in program_names]

print("Cultural Soft Power Investment Optimization")
print("=" * 50)
print(f"Total Budget: ${total_budget} million")
print(f"Number of Programs: {n_programs}")
print("\nProgram Details:")
for i, name in enumerate(program_names):
    print(f"{i+1}. {name}: Efficiency={efficiency_rates[i]:.1f}, Min=${min_investments[i]:.1f}M")

# Define the objective function (negative because we want to maximize)
def objective_function(x):
    """
    Objective function: Maximize total soft power influence
    Uses square root to model diminishing returns
    """
    total_influence = 0
    for i in range(len(x)):
        if x[i] > 0:
            total_influence += efficiency_rates[i] * np.sqrt(x[i])
    return -total_influence  # Negative for minimization

# Define constraints
def budget_constraint(x):
    """Budget constraint: total investment <= total budget"""
    return total_budget - np.sum(x)

def min_investment_constraints(x):
    """Minimum investment constraints for each program"""
    constraints = []
    for i in range(len(x)):
        if x[i] > 0:  # If we invest in program i
            constraints.append(x[i] - min_investments[i])
    return np.array(constraints) if constraints else np.array([0])

# Set up optimization constraints
constraints = [
    {'type': 'ineq', 'fun': budget_constraint},
]

# Bounds for each variable (0 to total budget)
bounds = [(0, total_budget) for _ in range(n_programs)]

# Initial guess (equal distribution)
x0 = np.array([total_budget / n_programs] * n_programs)

# Solve the optimization problem
print("\nSolving optimization problem...")
result = minimize(objective_function, x0, method='SLSQP', bounds=bounds, constraints=constraints)

if result.success:
    optimal_allocation = result.x
    max_influence = -result.fun
    
    print("\n" + "="*50)
    print("OPTIMIZATION RESULTS")
    print("="*50)
    print(f"Maximum Total Soft Power Influence: {max_influence:.2f}")
    print(f"Total Budget Used: ${np.sum(optimal_allocation):.2f} million")
    print(f"Budget Utilization: {np.sum(optimal_allocation)/total_budget*100:.1f}%")
    
    print("\nOptimal Investment Allocation:")
    results_df = pd.DataFrame({
        'Program': program_names,
        'Investment ($M)': optimal_allocation,
        'Investment (%)': optimal_allocation/total_budget*100,
        'Efficiency Rate': efficiency_rates,
        'Expected Influence': [efficiency_rates[i] * np.sqrt(optimal_allocation[i]) if optimal_allocation[i] > 0 else 0 
                              for i in range(n_programs)]
    })
    
    # Sort by investment amount (descending)
    results_df = results_df.sort_values('Investment ($M)', ascending=False)
    
    for _, row in results_df.iterrows():
        print(f"{row['Program']:<25}: ${row['Investment ($M)']:>6.2f}M ({row['Investment (%)']:>5.1f}%) -> Influence: {row['Expected Influence']:>5.2f}")
    
else:
    print("Optimization failed!")
    print(result.message)

# Calculate some additional metrics for analysis
print("\n" + "="*50)
print("INVESTMENT ANALYSIS")
print("="*50)

# Calculate ROI (Return on Investment) for each program
roi_analysis = []
for i, name in enumerate(program_names):
    if optimal_allocation[i] > 0:
        roi = (efficiency_rates[i] * np.sqrt(optimal_allocation[i])) / optimal_allocation[i]
        roi_analysis.append({
            'Program': name,
            'Investment': optimal_allocation[i],
            'ROI': roi,
            'Influence_per_Dollar': efficiency_rates[i] * np.sqrt(optimal_allocation[i]) / optimal_allocation[i]
        })

roi_df = pd.DataFrame(roi_analysis).sort_values('ROI', ascending=False)
print("\nReturn on Investment (ROI) Analysis:")
for _, row in roi_df.iterrows():
    print(f"{row['Program']:<25}: ROI = {row['ROI']:>5.2f}, Influence/$ = {row['Influence_per_Dollar']:>5.2f}")

# Compare with equal allocation strategy
equal_allocation = total_budget / n_programs
equal_total_influence = sum(efficiency_rates[i] * np.sqrt(equal_allocation) for i in range(n_programs))
improvement = (max_influence - equal_total_influence) / equal_total_influence * 100

print(f"\nComparison with Equal Allocation Strategy:")
print(f"Equal allocation influence: {equal_total_influence:.2f}")
print(f"Optimized allocation influence: {max_influence:.2f}")
print(f"Improvement: {improvement:.1f}%")

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

# 1. Investment Allocation Pie Chart
ax1 = plt.subplot(3, 3, 1)
# Filter out zero investments for cleaner pie chart
non_zero_investments = [(name, inv) for name, inv in zip(program_names, optimal_allocation) if inv > 0.01]
if non_zero_investments:
    names_nz, investments_nz = zip(*non_zero_investments)
    colors = plt.cm.Set3(np.linspace(0, 1, len(names_nz)))
    wedges, texts, autotexts = ax1.pie(investments_nz, labels=names_nz, autopct='%1.1f%%', 
                                      colors=colors, startangle=90)
    ax1.set_title('Optimal Investment Allocation\n($10M Total Budget)', fontsize=12, fontweight='bold')

# 2. Investment vs Efficiency Scatter Plot
ax2 = plt.subplot(3, 3, 2)
scatter = ax2.scatter(efficiency_rates, optimal_allocation, 
                     s=[inv*50 for inv in optimal_allocation], 
                     c=optimal_allocation, cmap='viridis', alpha=0.7)
ax2.set_xlabel('Program Efficiency Rate')
ax2.set_ylabel('Investment Amount ($M)')
ax2.set_title('Investment vs Efficiency\n(Bubble size = Investment)', fontsize=12, fontweight='bold')
plt.colorbar(scatter, ax=ax2, label='Investment ($M)')

# Add program labels
for i, name in enumerate(program_names):
    if optimal_allocation[i] > 0.1:  # Only label significant investments
        ax2.annotate(name.split()[0], (efficiency_rates[i], optimal_allocation[i]), 
                    xytext=(5, 5), textcoords='offset points', fontsize=8)

# 3. Expected Influence by Program
ax3 = plt.subplot(3, 3, 3)
expected_influence = [efficiency_rates[i] * np.sqrt(optimal_allocation[i]) if optimal_allocation[i] > 0 else 0 
                     for i in range(n_programs)]
bars = ax3.bar(range(len(program_names)), expected_influence, color='skyblue', edgecolor='navy')
ax3.set_xlabel('Programs')
ax3.set_ylabel('Expected Influence')
ax3.set_title('Expected Soft Power Influence\nby Program', fontsize=12, fontweight='bold')
ax3.set_xticks(range(len(program_names)))
ax3.set_xticklabels([name.split()[0] for name in program_names], rotation=45, ha='right')

# 4. ROI Analysis
ax4 = plt.subplot(3, 3, 4)
roi_values = [expected_influence[i]/optimal_allocation[i] if optimal_allocation[i] > 0 else 0 
              for i in range(n_programs)]
bars = ax4.bar(range(len(program_names)), roi_values, color='lightcoral', edgecolor='darkred')
ax4.set_xlabel('Programs')
ax4.set_ylabel('ROI (Influence per $M)')
ax4.set_title('Return on Investment\nby Program', fontsize=12, fontweight='bold')
ax4.set_xticks(range(len(program_names)))
ax4.set_xticklabels([name.split()[0] for name in program_names], rotation=45, ha='right')

# 5. Diminishing Returns Visualization
ax5 = plt.subplot(3, 3, 5)
investment_range = np.linspace(0.1, 5, 100)
for i in range(min(4, len(program_names))):  # Show top 4 programs
    influence_curve = efficiency_rates[i] * np.sqrt(investment_range)
    ax5.plot(investment_range, influence_curve, label=program_names[i].split()[0], linewidth=2)
    # Mark optimal point
    if optimal_allocation[i] > 0.1:
        optimal_influence = efficiency_rates[i] * np.sqrt(optimal_allocation[i])
        ax5.plot(optimal_allocation[i], optimal_influence, 'ro', markersize=8)

ax5.set_xlabel('Investment Amount ($M)')
ax5.set_ylabel('Soft Power Influence')
ax5.set_title('Diminishing Returns Curves\n(Red dots = Optimal points)', fontsize=12, fontweight='bold')
ax5.legend()
ax5.grid(True, alpha=0.3)

# 6. Comparison: Optimal vs Equal Allocation
ax6 = plt.subplot(3, 3, 6)
equal_allocation_array = np.full(n_programs, total_budget/n_programs)
equal_influence = [efficiency_rates[i] * np.sqrt(equal_allocation_array[i]) for i in range(n_programs)]

x_pos = np.arange(len(program_names))
width = 0.35

bars1 = ax6.bar(x_pos - width/2, expected_influence, width, label='Optimized', color='lightblue', alpha=0.8)
bars2 = ax6.bar(x_pos + width/2, equal_influence, width, label='Equal Allocation', color='lightgreen', alpha=0.8)

ax6.set_xlabel('Programs')
ax6.set_ylabel('Expected Influence')
ax6.set_title('Optimized vs Equal Allocation\nStrategy Comparison', fontsize=12, fontweight='bold')
ax6.set_xticks(x_pos)
ax6.set_xticklabels([name.split()[0] for name in program_names], rotation=45, ha='right')
ax6.legend()

# 7. Cumulative Investment Distribution
ax7 = plt.subplot(3, 3, 7)
sorted_indices = np.argsort(optimal_allocation)[::-1]
cumulative_investment = np.cumsum([optimal_allocation[i] for i in sorted_indices])
cumulative_percentage = cumulative_investment / total_budget * 100

ax7.plot(range(1, len(cumulative_percentage)+1), cumulative_percentage, 'bo-', linewidth=2, markersize=6)
ax7.fill_between(range(1, len(cumulative_percentage)+1), cumulative_percentage, alpha=0.3)
ax7.set_xlabel('Number of Programs (Ranked by Investment)')
ax7.set_ylabel('Cumulative Investment (%)')
ax7.set_title('Cumulative Investment Distribution\n(Pareto Analysis)', fontsize=12, fontweight='bold')
ax7.grid(True, alpha=0.3)
ax7.set_ylim(0, 100)

# Add 80% line for Pareto principle
ax7.axhline(y=80, color='red', linestyle='--', alpha=0.7, label='80% threshold')
ax7.legend()

# 8. Efficiency vs Minimum Investment Requirements
ax8 = plt.subplot(3, 3, 8)
scatter = ax8.scatter(min_investments, efficiency_rates, s=optimal_allocation*50, 
                     c=optimal_allocation, cmap='plasma', alpha=0.7, edgecolors='black', linewidth=1)
ax8.set_xlabel('Minimum Investment Requirement ($M)')
ax8.set_ylabel('Program Efficiency Rate')
ax8.set_title('Efficiency vs Requirements\n(Bubble size = Optimal investment)', fontsize=12, fontweight='bold')
plt.colorbar(scatter, ax=ax8, label='Optimal Investment ($M)')

# Add program labels
for i, name in enumerate(program_names):
    ax8.annotate(name.split()[0], (min_investments[i], efficiency_rates[i]), 
                xytext=(5, 5), textcoords='offset points', fontsize=8)

# 9. Summary Statistics
ax9 = plt.subplot(3, 3, 9)
ax9.axis('off')

# Calculate summary statistics
total_influence = sum(expected_influence)
avg_roi = np.mean([roi for roi in roi_values if roi > 0])
programs_funded = sum(1 for inv in optimal_allocation if inv > 0.01)
top_program_idx = np.argmax(optimal_allocation)
top_program = program_names[top_program_idx]

summary_text = f"""
OPTIMIZATION SUMMARY

Total Soft Power Influence: {total_influence:.2f}
Budget Utilization: {np.sum(optimal_allocation)/total_budget*100:.1f}%
Programs Funded: {programs_funded}/{n_programs}
Average ROI: {avg_roi:.2f}

TOP INVESTMENT:
{top_program}
${optimal_allocation[top_program_idx]:.2f}M ({optimal_allocation[top_program_idx]/total_budget*100:.1f}%)

IMPROVEMENT OVER EQUAL ALLOCATION:
{improvement:.1f}%

KEY INSIGHTS:
• Higher efficiency programs receive more funding
• Diminishing returns limit individual investments
• Portfolio diversification maintains 
  balanced influence across domains
• ROI optimization drives allocation decisions
"""

ax9.text(0.05, 0.95, summary_text, transform=ax9.transAxes, fontsize=10, 
         verticalalignment='top', fontfamily='monospace',
         bbox=dict(boxstyle='round,pad=0.5', facecolor='lightgray', alpha=0.8))

plt.tight_layout()
plt.savefig('cultural_investment_optimization.png', dpi=300, bbox_inches='tight')
plt.show()

# Additional Analysis: Sensitivity Analysis
print("\n" + "="*50)
print("SENSITIVITY ANALYSIS")
print("="*50)

# Test different budget scenarios
budget_scenarios = [5, 7.5, 10, 12.5, 15]  # Different budget levels in millions
scenario_results = []

for budget in budget_scenarios:
    # Re-solve optimization for each budget level
    def budget_constraint_scenario(x):
        return budget - np.sum(x)
    
    constraints_scenario = [{'type': 'ineq', 'fun': budget_constraint_scenario}]
    bounds_scenario = [(0, budget) for _ in range(n_programs)]
    x0_scenario = np.array([budget / n_programs] * n_programs)
    
    result_scenario = minimize(objective_function, x0_scenario, method='SLSQP', 
                              bounds=bounds_scenario, constraints=constraints_scenario)
    
    if result_scenario.success:
        total_influence_scenario = -result_scenario.fun
        scenario_results.append({
            'Budget': budget,
            'Total_Influence': total_influence_scenario,
            'Influence_per_Dollar': total_influence_scenario / budget
        })

# Plot sensitivity analysis
plt.figure(figsize=(12, 5))

# Budget vs Total Influence
plt.subplot(1, 2, 1)
budgets = [r['Budget'] for r in scenario_results]
influences = [r['Total_Influence'] for r in scenario_results]
plt.plot(budgets, influences, 'bo-', linewidth=2, markersize=8)
plt.xlabel('Budget ($M)')
plt.ylabel('Total Soft Power Influence')
plt.title('Budget vs Total Influence\n(Sensitivity Analysis)', fontweight='bold')
plt.grid(True, alpha=0.3)

# Budget vs Influence per Dollar (Efficiency)
plt.subplot(1, 2, 2)
efficiency_curve = [r['Influence_per_Dollar'] for r in scenario_results]
plt.plot(budgets, efficiency_curve, 'ro-', linewidth=2, markersize=8)
plt.xlabel('Budget ($M)')
plt.ylabel('Influence per Dollar')
plt.title('Budget Efficiency Curve\n(Diminishing Returns at Scale)', fontweight='bold')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("Sensitivity Analysis Results:")
for result in scenario_results:
    print(f"Budget: ${result['Budget']:>5.1f}M -> Total Influence: {result['Total_Influence']:>6.2f}, Efficiency: {result['Influence_per_Dollar']:>4.2f}")

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

class MultilateralDiplomacy:
    """
    A class to model and optimize multilateral diplomatic strategies
    using game-theoretic principles.
    """
    
    def __init__(self, countries, alpha, beta, gamma):
        """
        Initialize the diplomatic game.
        
        Parameters:
        - countries: list of country names
        - alpha: array of benefit coefficients from others' cooperation
        - beta: array of cost coefficients of own cooperation  
        - gamma: array of synergy coefficients
        """
        self.countries = countries
        self.n_countries = len(countries)
        self.alpha = np.array(alpha)
        self.beta = np.array(beta)
        self.gamma = np.array(gamma)
        
    def payoff(self, x, country_idx):
        """
        Calculate payoff for a specific country given all cooperation levels.
        
        Payoff function: U_i = α_i * Σ(x_j) - β_i * x_i^2 + γ_i * x_i * Σ(x_j)
        where j ≠ i
        """
        others_cooperation = np.sum(x) - x[country_idx]
        own_cooperation = x[country_idx]
        
        payoff = (self.alpha[country_idx] * others_cooperation - 
                 self.beta[country_idx] * own_cooperation**2 + 
                 self.gamma[country_idx] * own_cooperation * others_cooperation)
        
        return payoff
    
    def best_response(self, x, country_idx):
        """
        Calculate best response for country i given others' strategies.
        
        Taking derivative of payoff w.r.t. x_i and setting to zero:
        dU_i/dx_i = -2β_i * x_i + γ_i * Σ(x_j) = 0
        
        Optimal x_i = γ_i * Σ(x_j) / (2β_i)
        """
        others_cooperation = np.sum(x) - x[country_idx]
        
        if self.beta[country_idx] == 0:
            return 1.0  # Maximum cooperation if no cost
        
        optimal_x = (self.gamma[country_idx] * others_cooperation) / (2 * self.beta[country_idx])
        
        # Ensure cooperation level is between 0 and 1
        return np.clip(optimal_x, 0, 1)
    
    def nash_equilibrium_system(self, x):
        """
        System of equations for Nash equilibrium.
        Each country plays their best response given others' strategies.
        """
        equations = []
        for i in range(self.n_countries):
            best_resp = self.best_response(x, i)
            equations.append(x[i] - best_resp)
        return equations
    
    def find_nash_equilibrium(self, initial_guess=None):
        """
        Find Nash equilibrium using iterative best response method.
        """
        if initial_guess is None:
            x = np.random.rand(self.n_countries) * 0.5
        else:
            x = np.array(initial_guess)
        
        # Use scipy's fsolve to find the equilibrium
        equilibrium = fsolve(self.nash_equilibrium_system, x)
        
        # Ensure all values are between 0 and 1
        equilibrium = np.clip(equilibrium, 0, 1)
        
        return equilibrium
    
    def social_welfare_optimization(self):
        """
        Find socially optimal solution by maximizing total welfare.
        
        Social welfare = Σ U_i(x)
        """
        def negative_social_welfare(x):
            total_welfare = sum(self.payoff(x, i) for i in range(self.n_countries))
            return -total_welfare
        
        # Constraints: cooperation levels between 0 and 1
        constraints = [{'type': 'ineq', 'fun': lambda x: x},
                      {'type': 'ineq', 'fun': lambda x: 1 - x}]
        
        # Initial guess
        x0 = np.ones(self.n_countries) * 0.5
        
        result = minimize(negative_social_welfare, x0, method='SLSQP', 
                         constraints=constraints)
        
        return result.x if result.success else None
    
    def analyze_stability(self, equilibrium):
        """
        Analyze the stability of the equilibrium by computing payoff changes
        from unilateral deviations.
        """
        stability_analysis = {}
        
        for i, country in enumerate(self.countries):
            # Calculate current payoff
            current_payoff = self.payoff(equilibrium, i)
            
            # Test deviations
            deviations = np.linspace(0, 1, 21)
            deviation_payoffs = []
            
            for dev in deviations:
                x_dev = equilibrium.copy()
                x_dev[i] = dev
                deviation_payoffs.append(self.payoff(x_dev, i))
            
            stability_analysis[country] = {
                'current_payoff': current_payoff,
                'deviations': deviations,
                'deviation_payoffs': deviation_payoffs,
                'is_stable': all(p <= current_payoff + 1e-6 for p in deviation_payoffs)
            }
        
        return stability_analysis

# Define the diplomatic scenario
countries = ['USA', 'EU', 'China', 'Japan']

# Parameters representing different country characteristics
# Alpha: benefit from others' cooperation (trade opportunities)
alpha = [0.8, 0.9, 0.7, 0.85]  # EU most benefits from others' cooperation

# Beta: cost of own cooperation (domestic adjustment costs)
beta = [0.6, 0.4, 0.8, 0.5]    # China has highest domestic adjustment costs

# Gamma: synergy effects (network benefits)
gamma = [0.3, 0.4, 0.2, 0.35]  # EU gets most synergy benefits

# Create the diplomatic game
game = MultilateralDiplomacy(countries, alpha, beta, gamma)

print("=== Multilateral Trade Agreement Negotiation Analysis ===\n")

# Find Nash Equilibrium
nash_eq = game.find_nash_equilibrium()
print("Nash Equilibrium Cooperation Levels:")
for i, country in enumerate(countries):
    print(f"{country}: {nash_eq[i]:.3f}")

print(f"\nNash Equilibrium Payoffs:")
nash_payoffs = [game.payoff(nash_eq, i) for i in range(len(countries))]
for i, country in enumerate(countries):
    print(f"{country}: {nash_payoffs[i]:.3f}")

# Find Social Optimum
social_opt = game.social_welfare_optimization()
print(f"\nSocial Optimum Cooperation Levels:")
for i, country in enumerate(countries):
    print(f"{country}: {social_opt[i]:.3f}")

print(f"\nSocial Optimum Payoffs:")
social_payoffs = [game.payoff(social_opt, i) for i in range(len(countries))]
for i, country in enumerate(countries):
    print(f"{country}: {social_payoffs[i]:.3f}")

# Calculate welfare metrics
nash_total_welfare = sum(nash_payoffs)
social_total_welfare = sum(social_payoffs)
welfare_gap = social_total_welfare - nash_total_welfare

print(f"\nWelfare Analysis:")
print(f"Nash Equilibrium Total Welfare: {nash_total_welfare:.3f}")
print(f"Social Optimum Total Welfare: {social_total_welfare:.3f}")
print(f"Welfare Gap (Price of Anarchy): {welfare_gap:.3f}")
print(f"Efficiency Loss: {(welfare_gap/social_total_welfare)*100:.2f}%")

# Stability Analysis
stability = game.analyze_stability(nash_eq)
print(f"\nStability Analysis:")
for country, analysis in stability.items():
    status = "Stable" if analysis['is_stable'] else "Unstable"
    print(f"{country}: {status} (Current payoff: {analysis['current_payoff']:.3f})")

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

# 1. Cooperation Levels Comparison
ax1 = plt.subplot(2, 3, 1)
x_pos = np.arange(len(countries))
width = 0.35

plt.bar(x_pos - width/2, nash_eq, width, label='Nash Equilibrium', 
        alpha=0.8, color='steelblue')
plt.bar(x_pos + width/2, social_opt, width, label='Social Optimum', 
        alpha=0.8, color='coral')

plt.xlabel('Countries')
plt.ylabel('Cooperation Level')
plt.title('Cooperation Levels: Nash vs Social Optimum')
plt.xticks(x_pos, countries, rotation=45)
plt.legend()
plt.grid(True, alpha=0.3)

# 2. Payoff Comparison
ax2 = plt.subplot(2, 3, 2)
plt.bar(x_pos - width/2, nash_payoffs, width, label='Nash Equilibrium', 
        alpha=0.8, color='steelblue')
plt.bar(x_pos + width/2, social_payoffs, width, label='Social Optimum', 
        alpha=0.8, color='coral')

plt.xlabel('Countries')
plt.ylabel('Payoff')
plt.title('Individual Payoffs: Nash vs Social Optimum')
plt.xticks(x_pos, countries, rotation=45)
plt.legend()
plt.grid(True, alpha=0.3)

# 3. Stability Analysis Heatmap
ax3 = plt.subplot(2, 3, 3)
stability_data = []
for country in countries:
    stability_data.append(stability[country]['deviation_payoffs'])

stability_df = pd.DataFrame(stability_data, 
                           columns=[f'{d:.2f}' for d in stability[countries[0]]['deviations']], 
                           index=countries)

sns.heatmap(stability_df, annot=False, cmap='RdYlBu_r', center=0, ax=ax3)
plt.title('Payoff Landscape for Unilateral Deviations')
plt.xlabel('Cooperation Level')
plt.ylabel('Countries')

# 4. Individual Deviation Analysis
ax4 = plt.subplot(2, 3, 4)
for i, country in enumerate(countries):
    analysis = stability[country]
    plt.plot(analysis['deviations'], analysis['deviation_payoffs'], 
             label=country, marker='o', linewidth=2, markersize=4)
    # Mark current equilibrium point
    current_coop = nash_eq[i]
    current_payoff = analysis['current_payoff']
    plt.scatter([current_coop], [current_payoff], 
               s=100, marker='*', color='red', zorder=5)

plt.xlabel('Cooperation Level')
plt.ylabel('Payoff')
plt.title('Payoff Functions and Equilibrium Points')
plt.legend()
plt.grid(True, alpha=0.3)

# 5. Parameter Sensitivity Analysis
ax5 = plt.subplot(2, 3, 5)
# Vary beta (cost parameter) and see how equilibrium changes
beta_range = np.linspace(0.1, 1.5, 20)
equilibrium_paths = {country: [] for country in countries}

for beta_val in beta_range:
    # Create temporary game with modified beta for all countries
    temp_beta = [beta_val] * len(countries)
    temp_game = MultilateralDiplomacy(countries, alpha, temp_beta, gamma)
    temp_eq = temp_game.find_nash_equilibrium()
    
    for i, country in enumerate(countries):
        equilibrium_paths[country].append(temp_eq[i])

for country in countries:
    plt.plot(beta_range, equilibrium_paths[country], 
             label=country, marker='o', linewidth=2)

plt.xlabel('Cost Parameter (β)')
plt.ylabel('Equilibrium Cooperation Level')
plt.title('Sensitivity to Cooperation Costs')
plt.legend()
plt.grid(True, alpha=0.3)

# 6. 3D Visualization of Two-Player Interaction
ax6 = plt.subplot(2, 3, 6, projection='3d')

# Create a simplified 2-player version for 3D visualization
x1_range = np.linspace(0, 1, 30)
x2_range = np.linspace(0, 1, 30)
X1, X2 = np.meshgrid(x1_range, x2_range)

# Calculate payoffs for USA (player 1) across different strategy combinations
Z = np.zeros_like(X1)
for i in range(len(x1_range)):
    for j in range(len(x2_range)):
        x_temp = np.array([X1[i,j], X2[i,j], nash_eq[2], nash_eq[3]])
        Z[i,j] = game.payoff(x_temp, 0)  # USA's payoff

surface = ax6.plot_surface(X1, X2, Z, cmap='viridis', alpha=0.7)
ax6.scatter([nash_eq[0]], [nash_eq[1]], [nash_payoffs[0]], 
           color='red', s=100, label='Nash Equilibrium')

ax6.set_xlabel('USA Cooperation')
ax6.set_ylabel('EU Cooperation') 
ax6.set_zlabel('USA Payoff')
ax6.set_title('USA Payoff Landscape\n(China & Japan at Nash levels)')

plt.tight_layout()
plt.show()

# Additional Analysis: Coalition Formation
print("\n=== Coalition Analysis ===")

# Analyze potential bilateral coalitions
coalitions = [('USA', 'EU'), ('USA', 'China'), ('USA', 'Japan'), 
              ('EU', 'China'), ('EU', 'Japan'), ('China', 'Japan')]

coalition_benefits = {}
for coalition in coalitions:
    idx1, idx2 = countries.index(coalition[0]), countries.index(coalition[1])
    
    # Calculate joint optimization for the coalition
    def coalition_welfare(x):
        x_full = nash_eq.copy()
        x_full[idx1], x_full[idx2] = x[0], x[1]
        return -(game.payoff(x_full, idx1) + game.payoff(x_full, idx2))
    
    from scipy.optimize import minimize
    result = minimize(coalition_welfare, [nash_eq[idx1], nash_eq[idx2]], 
                     bounds=[(0, 1), (0, 1)])
    
    if result.success:
        optimal_coop = result.x
        x_coalition = nash_eq.copy()
        x_coalition[idx1], x_coalition[idx2] = optimal_coop[0], optimal_coop[1]
        
        coalition_payoff1 = game.payoff(x_coalition, idx1)
        coalition_payoff2 = game.payoff(x_coalition, idx2)
        
        benefit1 = coalition_payoff1 - nash_payoffs[idx1]
        benefit2 = coalition_payoff2 - nash_payoffs[idx2]
        
        coalition_benefits[coalition] = {
            'cooperation': optimal_coop,
            'benefits': (benefit1, benefit2),
            'total_benefit': benefit1 + benefit2
        }

print("\nPotential Coalition Benefits:")
for coalition, data in sorted(coalition_benefits.items(), 
                             key=lambda x: x[1]['total_benefit'], reverse=True):
    print(f"{coalition[0]}-{coalition[1]}: "
          f"Benefits = ({data['benefits'][0]:.3f}, {data['benefits'][1]:.3f}), "
          f"Total = {data['total_benefit']:.3f}")