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}")

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

# Set up the problem data
np.random.seed(42)

# Create realistic country data
countries = [
    'USA', 'China', 'Germany', 'Japan', 'UK', 'France', 'India', 'Italy',
    'Brazil', 'Canada', 'Russia', 'South Korea', 'Spain', 'Australia', 'Mexico'
]

n_countries = len(countries)

# Generate realistic data
# Voting weights (normalized to sum to 1)
voting_weights = np.array([0.15, 0.14, 0.08, 0.07, 0.06, 0.06, 0.05, 0.05,
                          0.04, 0.04, 0.04, 0.04, 0.03, 0.03, 0.02])

# Probability of support (0-1 scale)
support_probability = np.array([0.7, 0.4, 0.9, 0.8, 0.85, 0.9, 0.6, 0.8,
                               0.75, 0.9, 0.3, 0.7, 0.85, 0.8, 0.7])

# Cost to influence (arbitrary units, higher = more expensive)
influence_cost = np.array([100, 120, 60, 70, 50, 55, 80, 45,
                          40, 35, 90, 50, 40, 45, 35])

# Total budget
total_budget = 300

# Create DataFrame for better visualization
df = pd.DataFrame({
    'Country': countries,
    'Voting_Weight': voting_weights,
    'Support_Probability': support_probability,
    'Influence_Cost': influence_cost,
    'Expected_Value': voting_weights * support_probability,
    'Value_per_Cost': (voting_weights * support_probability) / influence_cost
})

print("=== INTERNATIONAL ORGANIZATION INFLUENCE OPTIMIZATION ===")
print("\nCountry Data:")
print(df.round(4))

# Method 1: Greedy Algorithm (Value per Cost)
def greedy_solution(weights, probs, costs, budget):
    """
    Greedy algorithm: select countries with highest value-per-cost ratio
    """
    n = len(weights)
    value_per_cost = (weights * probs) / costs
    
    # Sort by value per cost (descending)
    sorted_indices = np.argsort(-value_per_cost)
    
    selected = np.zeros(n, dtype=bool)
    total_cost = 0
    total_value = 0
    
    for idx in sorted_indices:
        if total_cost + costs[idx] <= budget:
            selected[idx] = True
            total_cost += costs[idx]
            total_value += weights[idx] * probs[idx]
    
    return selected, total_value, total_cost

# Method 2: Dynamic Programming (Exact solution for 0/1 Knapsack)
def dp_knapsack(weights, probs, costs, budget):
    """
    Dynamic programming solution for 0/1 knapsack problem
    """
    n = len(weights)
    values = weights * probs
    
    # Scale values to integers for DP
    scale_factor = 10000
    scaled_values = (values * scale_factor).astype(int)
    budget_int = int(budget)
    
    # DP table: dp[i][w] = maximum value using first i items with weight limit w
    dp = np.zeros((n + 1, budget_int + 1))
    
    for i in range(1, n + 1):
        for w in range(budget_int + 1):
            cost_i = int(costs[i-1])
            if cost_i <= w:
                # Max of including or not including item i-1
                dp[i][w] = max(dp[i-1][w], 
                              dp[i-1][w-cost_i] + scaled_values[i-1])
            else:
                dp[i][w] = dp[i-1][w]
    
    # Backtrack to find selected items
    selected = np.zeros(n, dtype=bool)
    w = budget_int
    for i in range(n, 0, -1):
        if dp[i][w] != dp[i-1][w]:
            selected[i-1] = True
            w -= int(costs[i-1])
    
    total_value = np.sum(values[selected])
    total_cost = np.sum(costs[selected])
    
    return selected, total_value, total_cost

# Method 3: Brute Force (for comparison, only feasible for small n)
def brute_force_solution(weights, probs, costs, budget):
    """
    Brute force: try all possible combinations (exponential time)
    Only use for small problems due to 2^n complexity
    """
    n = len(weights)
    if n > 20:  # Avoid exponential explosion
        return None, 0, 0
    
    best_value = 0
    best_selection = np.zeros(n, dtype=bool)
    best_cost = 0
    
    # Try all 2^n combinations
    for r in range(1, n + 1):
        for combo in combinations(range(n), r):
            total_cost = sum(costs[i] for i in combo)
            if total_cost <= budget:
                total_value = sum(weights[i] * probs[i] for i in combo)
                if total_value > best_value:
                    best_value = total_value
                    best_selection = np.zeros(n, dtype=bool)
                    best_selection[list(combo)] = True
                    best_cost = total_cost
    
    return best_selection, best_value, best_cost

print(f"\nBudget: {total_budget} units")
print("="*60)

# Solve using all methods
greedy_selected, greedy_value, greedy_cost = greedy_solution(
    voting_weights, support_probability, influence_cost, total_budget)

dp_selected, dp_value, dp_cost = dp_knapsack(
    voting_weights, support_probability, influence_cost, total_budget)

brute_selected, brute_value, brute_cost = brute_force_solution(
    voting_weights, support_probability, influence_cost, total_budget)

# Display results
methods = ['Greedy', 'Dynamic Programming', 'Brute Force']
values = [greedy_value, dp_value, brute_value]
costs = [greedy_cost, dp_cost, brute_cost]
selections = [greedy_selected, dp_selected, brute_selected]

print("\n=== OPTIMIZATION RESULTS ===")
for i, method in enumerate(methods):
    if selections[i] is not None:
        selected_countries = [countries[j] for j in range(n_countries) if selections[i][j]]
        print(f"\n{method}:")
        print(f"  Expected Influence: {values[i]:.4f}")
        print(f"  Total Cost: {costs[i]:.1f} / {total_budget}")
        print(f"  Selected Countries: {', '.join(selected_countries)}")
        print(f"  Number of Countries: {sum(selections[i])}")

# Create comprehensive visualizations
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('International Organization Influence Optimization Analysis', fontsize=16, fontweight='bold')

# 1. Country Characteristics Scatter Plot
ax1 = axes[0, 0]
scatter = ax1.scatter(influence_cost, voting_weights, 
                     c=support_probability, cmap='RdYlGn', 
                     s=200, alpha=0.7, edgecolors='black', linewidth=1)
for i, country in enumerate(countries):
    ax1.annotate(country, (influence_cost[i], voting_weights[i]), 
                xytext=(5, 5), textcoords='offset points', fontsize=8)
ax1.set_xlabel('Influence Cost')
ax1.set_ylabel('Voting Weight')
ax1.set_title('Country Characteristics\n(Color = Support Probability)')
plt.colorbar(scatter, ax=ax1, label='Support Probability')
ax1.grid(True, alpha=0.3)

# 2. Value per Cost Analysis
ax2 = axes[0, 1]
value_per_cost = df['Value_per_Cost'].values
colors = ['green' if greedy_selected[i] else 'lightgray' for i in range(n_countries)]
bars = ax2.bar(range(n_countries), value_per_cost, color=colors, alpha=0.7, edgecolor='black')
ax2.set_xlabel('Country Index')
ax2.set_ylabel('Expected Value per Cost')
ax2.set_title('Value per Cost Ratio\n(Green = Selected by Greedy)')
ax2.set_xticks(range(n_countries))
ax2.set_xticklabels([c[:3] for c in countries], rotation=45)
ax2.grid(True, alpha=0.3)

# 3. Method Comparison
ax3 = axes[0, 2]
method_names = ['Greedy', 'DP', 'Brute Force']
method_values = [greedy_value, dp_value, brute_value]
method_costs = [greedy_cost, dp_cost, brute_cost]

x = np.arange(len(method_names))
width = 0.35

bars1 = ax3.bar(x - width/2, method_values, width, label='Expected Influence', 
                color='skyblue', alpha=0.8, edgecolor='black')
bars2 = ax3.bar(x + width/2, np.array(method_costs)/total_budget, width, 
                label='Budget Utilization', color='lightcoral', alpha=0.8, edgecolor='black')

ax3.set_xlabel('Optimization Method')
ax3.set_ylabel('Normalized Value')
ax3.set_title('Method Performance Comparison')
ax3.set_xticks(x)
ax3.set_xticklabels(method_names)
ax3.legend()
ax3.grid(True, alpha=0.3)

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

# 4. Selection Comparison Heatmap
ax4 = axes[1, 0]
selection_matrix = np.array([greedy_selected.astype(int), 
                           dp_selected.astype(int), 
                           brute_selected.astype(int)])
sns.heatmap(selection_matrix, annot=True, fmt='d', cmap='RdYlGn', 
            xticklabels=[c[:3] for c in countries], 
            yticklabels=['Greedy', 'DP', 'Brute Force'],
            ax=ax4, cbar_kws={'label': 'Selected (1) / Not Selected (0)'})
ax4.set_title('Country Selection by Method')
ax4.set_xlabel('Countries')

# 5. Budget Allocation Analysis
ax5 = axes[1, 1]
# Show cost breakdown for DP solution (optimal)
dp_costs = influence_cost * dp_selected
selected_mask = dp_selected
selected_countries_short = [countries[i][:3] for i in range(n_countries) if selected_mask[i]]
selected_costs = dp_costs[selected_mask]

wedges, texts, autotexts = ax5.pie(selected_costs, labels=selected_countries_short, 
                                  autopct='%1.1f%%', startangle=90)
ax5.set_title(f'Budget Allocation (DP Solution)\nTotal: {dp_cost:.1f}/{total_budget}')

# 6. Influence Distribution
ax6 = axes[1, 2]
dp_influence = (voting_weights * support_probability) * dp_selected
selected_influence = dp_influence[selected_mask]

bars = ax6.bar(selected_countries_short, selected_influence, 
               color='lightgreen', alpha=0.8, edgecolor='black')
ax6.set_xlabel('Selected Countries')
ax6.set_ylabel('Expected Influence')
ax6.set_title('Expected Influence by Country\n(DP Solution)')
ax6.tick_params(axis='x', rotation=45)
ax6.grid(True, alpha=0.3)

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

plt.tight_layout()
plt.show()

# Additional Analysis: Sensitivity Analysis
print("\n=== SENSITIVITY ANALYSIS ===")
print("How does the solution change with different budget levels?")

budgets = np.arange(100, 501, 50)
sensitivity_results = []

for budget in budgets:
    _, value, cost = dp_knapsack(voting_weights, support_probability, influence_cost, budget)
    sensitivity_results.append({
        'Budget': budget,
        'Expected_Influence': value,
        'Utilization': cost/budget
    })

sensitivity_df = pd.DataFrame(sensitivity_results)
print("\nBudget Sensitivity:")
print(sensitivity_df.round(4))

# Plot sensitivity analysis
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

ax1.plot(sensitivity_df['Budget'], sensitivity_df['Expected_Influence'], 
         marker='o', linewidth=2, markersize=8, color='blue')
ax1.set_xlabel('Budget')
ax1.set_ylabel('Expected Influence')
ax1.set_title('Expected Influence vs Budget')
ax1.grid(True, alpha=0.3)

ax2.plot(sensitivity_df['Budget'], sensitivity_df['Utilization'], 
         marker='s', linewidth=2, markersize=8, color='red')
ax2.set_xlabel('Budget')
ax2.set_ylabel('Budget Utilization Rate')
ax2.set_title('Budget Utilization vs Budget Level')
ax2.grid(True, alpha=0.3)
ax2.set_ylim(0, 1.1)

plt.tight_layout()
plt.show()

print("\n=== KEY INSIGHTS ===")
print("1. Optimal Strategy: The Dynamic Programming solution provides the mathematically optimal allocation")
print("2. Greedy Performance: Greedy algorithm provides good approximation with much lower computational cost")
print("3. Budget Efficiency: Higher budgets show diminishing returns due to discrete country selection")
print("4. Strategic Focus: Target countries with high value-per-cost ratios first")
print(f"5. Coalition Size: Optimal coalition includes {sum(dp_selected)} out of {n_countries} countries")

# Calculate coalition strength
total_possible_influence = np.sum(voting_weights * support_probability)
achieved_influence_ratio = dp_value / total_possible_influence

print(f"6. Influence Coverage: Achieving {achieved_influence_ratio:.1%} of maximum possible influence")
print(f"7. Cost Efficiency: {dp_value/dp_cost:.5f} influence units per cost unit")

# Diplomatic Resource Allocation Optimization
# A case study in strategic resource distribution

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

# Set up the problem parameters
print("=== Diplomatic Resource Allocation Optimization ===\n")

# Define regions
regions = ['Asia-Pacific', 'Europe', 'Africa', 'Latin America']

# Benefit coefficients (diplomatic value per unit of resource)
# Higher values indicate greater strategic importance
benefit_coefficients = np.array([
    8.5,  # Asia-Pacific: High economic potential
    7.2,  # Europe: Stable partnerships
    6.8,  # Africa: Emerging opportunities  
    5.9   # Latin America: Regional stability
])

print("Benefit Coefficients (Diplomatic Value per Unit):")
for i, region in enumerate(regions):
    print(f"  {region}: {benefit_coefficients[i]}")

# Constraint parameters
total_budget = 100.0  # Total diplomatic resource units
min_allocation = np.array([5, 8, 10, 5])  # Minimum engagement requirements
max_allocation = np.array([40, 35, 30, 25])  # Maximum capacity constraints

print(f"\nTotal Available Budget: {total_budget} resource units")
print("\nMinimum Allocation Requirements:")
for i, region in enumerate(regions):
    print(f"  {region}: {min_allocation[i]} units")
    
print("\nMaximum Allocation Capacity:")
for i, region in enumerate(regions):
    print(f"  {region}: {max_allocation[i]} units")

# Set up the linear programming problem
# We need to convert maximization to minimization by negating coefficients
c = -benefit_coefficients  # Negate for maximization

# Inequality constraints (Ax <= b)
# Budget constraint: x1 + x2 + x3 + x4 <= total_budget
A_ub = np.array([[1, 1, 1, 1]])  # Sum of all allocations
b_ub = np.array([total_budget])

# Bounds for variables (min_allocation <= xi <= max_allocation)
bounds = [(min_allocation[i], max_allocation[i]) for i in range(4)]

print(f"\nBounds for each region: {bounds}")

# Solve the optimization problem
print("\n" + "="*50)
print("SOLVING OPTIMIZATION PROBLEM...")
print("="*50)

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

if result.success:
    optimal_allocation = result.x
    optimal_value = -result.fun  # Convert back to maximization value
    
    print(f"\nOptimization Status: SUCCESS")
    print(f"Maximum Diplomatic Benefit: {optimal_value:.2f}")
    print(f"Total Resources Used: {sum(optimal_allocation):.2f} / {total_budget}")
    
    print(f"\nOptimal Resource Allocation:")
    for i, region in enumerate(regions):
        percentage = (optimal_allocation[i] / sum(optimal_allocation)) * 100
        benefit = optimal_allocation[i] * benefit_coefficients[i]
        print(f"  {region}: {optimal_allocation[i]:.1f} units ({percentage:.1f}%) -> Benefit: {benefit:.1f}")
else:
    print("Optimization failed:", result.message)

# Create comprehensive visualizations
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('Diplomatic Resource Allocation Analysis', fontsize=16, fontweight='bold')

# 1. Optimal allocation pie chart
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4']
wedges, texts, autotexts = ax1.pie(optimal_allocation, labels=regions, autopct='%1.1f%%', 
                                   colors=colors, startangle=90, explode=(0.05, 0.05, 0.05, 0.05))
ax1.set_title('Optimal Resource Distribution', fontweight='bold')

# Enhance pie chart text
for autotext in autotexts:
    autotext.set_color('white')
    autotext.set_fontweight('bold')

# 2. Resource allocation vs constraints
x_pos = np.arange(len(regions))
width = 0.25

bars1 = ax2.bar(x_pos - width, min_allocation, width, label='Minimum Required', 
                color='lightcoral', alpha=0.7)
bars2 = ax2.bar(x_pos, optimal_allocation, width, label='Optimal Allocation', 
                color='steelblue', alpha=0.8)
bars3 = ax2.bar(x_pos + width, max_allocation, width, label='Maximum Capacity', 
                color='lightgreen', alpha=0.7)

ax2.set_xlabel('Regions')
ax2.set_ylabel('Resource Units')
ax2.set_title('Allocation Comparison: Optimal vs Constraints')
ax2.set_xticks(x_pos)
ax2.set_xticklabels(regions, rotation=45, ha='right')
ax2.legend()
ax2.grid(axis='y', alpha=0.3)

# Add value labels on bars
for bars in [bars1, bars2, bars3]:
    for bar in bars:
        height = bar.get_height()
        ax2.text(bar.get_x() + bar.get_width()/2., height + 0.5,
                f'{height:.1f}', ha='center', va='bottom', fontsize=9)

# 3. Benefit analysis
benefit_per_region = optimal_allocation * benefit_coefficients
bars = ax3.bar(regions, benefit_per_region, color=colors, alpha=0.8)
ax3.set_xlabel('Regions')
ax3.set_ylabel('Total Diplomatic Benefit')
ax3.set_title('Diplomatic Benefit by Region')
ax3.set_xticklabels(regions, rotation=45, ha='right')
ax3.grid(axis='y', alpha=0.3)

# Add value labels
for i, bar in enumerate(bars):
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height + 0.5,
            f'{height:.1f}', ha='center', va='bottom', fontweight='bold')

# 4. Efficiency analysis (Benefit per unit resource)
efficiency = benefit_coefficients
bars = ax4.bar(regions, efficiency, color='orange', alpha=0.8)
ax4.set_xlabel('Regions')
ax4.set_ylabel('Benefit per Resource Unit')
ax4.set_title('Resource Efficiency by Region')
ax4.set_xticklabels(regions, rotation=45, ha='right')
ax4.grid(axis='y', alpha=0.3)

# Add value labels
for i, bar in enumerate(bars):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height + 0.1,
            f'{height:.1f}', ha='center', va='bottom', fontweight='bold')

plt.tight_layout()
plt.show()

# Create summary table
print("\n" + "="*70)
print("DETAILED ANALYSIS SUMMARY")
print("="*70)

summary_data = {
    'Region': regions,
    'Min Required': min_allocation,
    'Optimal Allocation': optimal_allocation,
    'Max Capacity': max_allocation,
    'Benefit Coefficient': benefit_coefficients,
    'Total Benefit': benefit_per_region,
    'Utilization (%)': (optimal_allocation / max_allocation) * 100
}

df_summary = pd.DataFrame(summary_data)
df_summary = df_summary.round(2)

print(df_summary.to_string(index=False))

# Sensitivity analysis
print(f"\n" + "="*50)
print("SENSITIVITY ANALYSIS")
print("="*50)

print(f"Budget Utilization: {(sum(optimal_allocation)/total_budget)*100:.1f}%")
print(f"Average Benefit per Unit: {optimal_value/sum(optimal_allocation):.2f}")

# Calculate slack in constraints
constraint_slack = total_budget - sum(optimal_allocation)
print(f"Remaining Budget (Slack): {constraint_slack:.1f} units")

# Check which regions are at their limits
print(f"\nRegions at Maximum Capacity:")
for i, region in enumerate(regions):
    if abs(optimal_allocation[i] - max_allocation[i]) < 0.1:
        print(f"  {region}: {optimal_allocation[i]:.1f} / {max_allocation[i]} units")

print(f"\nRegions at Minimum Requirement:")
for i, region in enumerate(regions):
    if abs(optimal_allocation[i] - min_allocation[i]) < 0.1:
        print(f"  {region}: {optimal_allocation[i]:.1f} / {min_allocation[i]} units")

print(f"\n" + "="*50)
print("OPTIMIZATION COMPLETE")
print("="*50)

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

# Set up the problem parameters
class ResourceDiplomacyOptimizer:
    def __init__(self):
        # Production capacity (barrels per day)
        self.total_capacity = 2_000_000
        
        # Trading partner characteristics
        self.partners = {
            'Country_A': {'price': 85.0, 'strategic_weight': 0.8, 'min_commitment': 300_000},
            'Country_B': {'price': 82.0, 'strategic_weight': 0.6, 'min_commitment': 200_000},
            'Country_C': {'price': 78.0, 'strategic_weight': 0.9, 'min_commitment': 400_000}
        }
        
        # Weights for multi-objective optimization
        self.economic_weight = 0.7  # Weight for economic returns
        self.diplomatic_weight = 0.3  # Weight for diplomatic benefits
        
    def objective_function(self, x):
        """
        Multi-objective function combining economic and diplomatic benefits
        x = [allocation_A, allocation_B, allocation_C]
        """
        allocation_A, allocation_B, allocation_C = x
        
        # Economic benefit calculation
        economic_benefit = (
            allocation_A * self.partners['Country_A']['price'] +
            allocation_B * self.partners['Country_B']['price'] +
            allocation_C * self.partners['Country_C']['price']
        )
        
        # Diplomatic benefit calculation (normalized strategic importance)
        diplomatic_benefit = (
            allocation_A * self.partners['Country_A']['strategic_weight'] +
            allocation_B * self.partners['Country_B']['strategic_weight'] +
            allocation_C * self.partners['Country_C']['strategic_weight']
        )
        
        # Combined objective (negative because we minimize)
        total_benefit = (
            self.economic_weight * economic_benefit +
            self.diplomatic_weight * diplomatic_benefit * 100_000  # Scale diplomatic benefits
        )
        
        return -total_benefit  # Negative for minimization
    
    def constraints(self):
        """Define constraints for the optimization problem"""
        constraints = []
        
        # Total capacity constraint
        constraints.append({
            'type': 'eq',
            'fun': lambda x: self.total_capacity - sum(x)
        })
        
        # Minimum commitment constraints
        constraints.append({
            'type': 'ineq',
            'fun': lambda x: x[0] - self.partners['Country_A']['min_commitment']
        })
        constraints.append({
            'type': 'ineq',
            'fun': lambda x: x[1] - self.partners['Country_B']['min_commitment']
        })
        constraints.append({
            'type': 'ineq',
            'fun': lambda x: x[2] - self.partners['Country_C']['min_commitment']
        })
        
        return constraints
    
    def solve_optimization(self):
        """Solve the resource allocation optimization problem"""
        # Initial guess (equal distribution after minimum commitments)
        remaining_capacity = self.total_capacity - sum(
            partner['min_commitment'] for partner in self.partners.values()
        )
        
        x0 = [
            self.partners['Country_A']['min_commitment'] + remaining_capacity/3,
            self.partners['Country_B']['min_commitment'] + remaining_capacity/3,
            self.partners['Country_C']['min_commitment'] + remaining_capacity/3
        ]
        
        # Bounds (non-negative allocations, maximum is total capacity)
        bounds = [(0, self.total_capacity) for _ in range(3)]
        
        # Solve optimization
        result = minimize(
            self.objective_function,
            x0,
            method='SLSQP',
            bounds=bounds,
            constraints=self.constraints()
        )
        
        return result
    
    def analyze_sensitivity(self, price_variations=None):
        """Analyze sensitivity to price changes"""
        if price_variations is None:
            price_variations = np.linspace(0.8, 1.2, 10)  # ±20% price variation
        
        results = []
        base_prices = [self.partners[country]['price'] for country in self.partners.keys()]
        
        for variation in price_variations:
            # Temporarily modify prices
            for country in self.partners.keys():
                self.partners[country]['price'] *= variation
            
            # Solve optimization
            result = self.solve_optimization()
            
            if result.success:
                allocation = result.x
                total_revenue = -result.fun
                results.append({
                    'price_multiplier': variation,
                    'country_A_allocation': allocation[0],
                    'country_B_allocation': allocation[1],
                    'country_C_allocation': allocation[2],
                    'total_benefit': total_revenue
                })
            
            # Restore original prices
            for i, country in enumerate(self.partners.keys()):
                self.partners[country]['price'] = base_prices[i]
        
        return pd.DataFrame(results)

# Initialize optimizer
optimizer = ResourceDiplomacyOptimizer()

# Solve the optimization problem
print("=== Resource Diplomacy Optimization Results ===")
result = optimizer.solve_optimization()

if result.success:
    optimal_allocation = result.x
    print(f"\nOptimization Status: {result.message}")
    print(f"\nOptimal Resource Allocation:")
    print(f"Country A: {optimal_allocation[0]:,.0f} barrels/day ({optimal_allocation[0]/optimizer.total_capacity*100:.1f}%)")
    print(f"Country B: {optimal_allocation[1]:,.0f} barrels/day ({optimal_allocation[1]/optimizer.total_capacity*100:.1f}%)")
    print(f"Country C: {optimal_allocation[2]:,.0f} barrels/day ({optimal_allocation[2]/optimizer.total_capacity*100:.1f}%)")
    print(f"\nTotal Allocation: {sum(optimal_allocation):,.0f} barrels/day")
    
    # Calculate economic and diplomatic benefits
    economic_benefit = sum(optimal_allocation[i] * list(optimizer.partners.values())[i]['price'] 
                          for i in range(3))
    diplomatic_benefit = sum(optimal_allocation[i] * list(optimizer.partners.values())[i]['strategic_weight'] 
                           for i in range(3))
    
    print(f"\nDaily Economic Benefit: ${economic_benefit:,.0f}")
    print(f"Diplomatic Benefit Index: {diplomatic_benefit:,.0f}")
    
else:
    print(f"Optimization failed: {result.message}")

# Perform sensitivity analysis
print("\n=== Conducting Sensitivity Analysis ===")
sensitivity_data = optimizer.analyze_sensitivity()

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

# 1. Optimal Allocation Pie Chart
ax1 = plt.subplot(2, 3, 1)
countries = ['Country A', 'Country B', 'Country C']
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']
wedges, texts, autotexts = ax1.pie(optimal_allocation, labels=countries, autopct='%1.1f%%', 
                                   colors=colors, startangle=90, explode=(0.05, 0.05, 0.05))
ax1.set_title('Optimal Resource Allocation\n(Barrels per Day)', fontsize=14, fontweight='bold')

# 2. Economic vs Diplomatic Trade-off
ax2 = plt.subplot(2, 3, 2)
prices = [optimizer.partners[country]['price'] for country in optimizer.partners.keys()]
strategic_weights = [optimizer.partners[country]['strategic_weight'] for country in optimizer.partners.keys()]
scatter = ax2.scatter(prices, strategic_weights, s=[alloc/5000 for alloc in optimal_allocation], 
                     c=colors, alpha=0.7, edgecolors='black', linewidth=2)
ax2.set_xlabel('Oil Price ($/barrel)', fontsize=12)
ax2.set_ylabel('Strategic Importance', fontsize=12)
ax2.set_title('Price vs Strategic Importance\n(Bubble size = Allocation)', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

for i, country in enumerate(countries):
    ax2.annotate(country, (prices[i], strategic_weights[i]), 
                xytext=(5, 5), textcoords='offset points', fontsize=10)

# 3. Sensitivity Analysis - Allocation Changes
ax3 = plt.subplot(2, 3, 3)
ax3.plot(sensitivity_data['price_multiplier'], sensitivity_data['country_A_allocation']/1000, 
         'o-', label='Country A', color=colors[0], linewidth=2, markersize=6)
ax3.plot(sensitivity_data['price_multiplier'], sensitivity_data['country_B_allocation']/1000, 
         's-', label='Country B', color=colors[1], linewidth=2, markersize=6)
ax3.plot(sensitivity_data['price_multiplier'], sensitivity_data['country_C_allocation']/1000, 
         '^-', label='Country C', color=colors[2], linewidth=2, markersize=6)
ax3.set_xlabel('Price Multiplier', fontsize=12)
ax3.set_ylabel('Allocation (Thousands of barrels/day)', fontsize=12)
ax3.set_title('Sensitivity to Price Changes', fontsize=14, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Total Benefit vs Price Multiplier
ax4 = plt.subplot(2, 3, 4)
ax4.plot(sensitivity_data['price_multiplier'], sensitivity_data['total_benefit']/1_000_000, 
         'o-', color='green', linewidth=3, markersize=8)
ax4.set_xlabel('Price Multiplier', fontsize=12)
ax4.set_ylabel('Total Benefit (Millions $)', fontsize=12)
ax4.set_title('Total Economic Benefit\nvs Price Variation', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

# 5. Constraint Analysis
ax5 = plt.subplot(2, 3, 5)
min_commitments = [optimizer.partners[country]['min_commitment'] for country in optimizer.partners.keys()]
actual_allocations = optimal_allocation

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

bars1 = ax5.bar(x_pos - width/2, [x/1000 for x in min_commitments], width, 
                label='Minimum Commitments', color='lightcoral', alpha=0.7)
bars2 = ax5.bar(x_pos + width/2, [x/1000 for x in actual_allocations], width, 
                label='Optimal Allocation', color='lightblue', alpha=0.7)

ax5.set_xlabel('Countries', fontsize=12)
ax5.set_ylabel('Allocation (Thousands of barrels/day)', fontsize=12)
ax5.set_title('Minimum Commitments vs\nOptimal Allocation', fontsize=14, fontweight='bold')
ax5.set_xticks(x_pos)
ax5.set_xticklabels(countries)
ax5.legend()
ax5.grid(True, alpha=0.3, axis='y')

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

# 6. Revenue Breakdown
ax6 = plt.subplot(2, 3, 6)
revenues = [optimal_allocation[i] * prices[i] for i in range(3)]
colors_revenue = ['#FF9999', '#66B2FF', '#99FF99']

bars = ax6.bar(countries, [r/1_000_000 for r in revenues], color=colors_revenue, alpha=0.8, edgecolor='black')
ax6.set_ylabel('Daily Revenue (Millions $)', fontsize=12)
ax6.set_title('Daily Revenue by Trading Partner', fontsize=14, fontweight='bold')
ax6.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for i, bar in enumerate(bars):
    height = bar.get_height()
    ax6.text(bar.get_x() + bar.get_width()/2., height + 0.5,
             f'${height:.1f}M', ha='center', va='bottom', fontsize=11, fontweight='bold')

plt.tight_layout()
plt.show()

# Summary statistics
print("\n=== Summary Statistics ===")
total_revenue = sum(revenues)
print(f"Total Daily Revenue: ${total_revenue:,.0f}")
print(f"Annual Revenue Projection: ${total_revenue * 365:,.0f}")

utilization_rates = [optimal_allocation[i] / optimizer.partners[list(optimizer.partners.keys())[i]]['min_commitment'] 
                    for i in range(3)]
print(f"\nUtilization above minimum commitments:")
for i, country in enumerate(countries):
    excess = optimal_allocation[i] - min_commitments[i]
    print(f"{country}: +{excess:,.0f} barrels/day ({(utilization_rates[i]-1)*100:.1f}% above minimum)")

# Pipeline and Maritime Transport Route Optimization
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
import pandas as pd
from matplotlib.patches import FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')

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

class TransportOptimizer:
    def __init__(self, platforms, refineries):
        self.platforms = platforms
        self.refineries = refineries
        self.m = len(platforms)  # number of platforms
        self.n = len(refineries)  # number of refineries
        
        # Generate supply and demand data
        self.supply = np.array([p['capacity'] for p in platforms])
        self.demand = np.array([r['demand'] for r in refineries])
        
        # Generate cost matrices based on distances and transport modes
        self.pipeline_costs = self._generate_pipeline_costs()
        self.ship_costs = self._generate_ship_costs()
        
    def _calculate_distance(self, p1, p2):
        """Calculate Euclidean distance between two points"""
        return np.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)
    
    def _generate_pipeline_costs(self):
        """Generate pipeline transportation costs based on distance"""
        costs = np.zeros((self.m, self.n))
        base_pipeline_cost = 15  # $/barrel per unit distance
        
        for i, platform in enumerate(self.platforms):
            for j, refinery in enumerate(self.refineries):
                distance = self._calculate_distance(platform['location'], refinery['location'])
                # Pipeline cost increases with distance, but has economies of scale
                costs[i, j] = base_pipeline_cost * distance * 0.8
                
                # Some routes may not have pipeline infrastructure
                if distance > 300:  # Long distances make pipelines impractical
                    costs[i, j] = float('inf')
                    
        return costs
    
    def _generate_ship_costs(self):
        """Generate maritime transportation costs"""
        costs = np.zeros((self.m, self.n))
        base_ship_cost = 8  # $/barrel per unit distance
        
        for i, platform in enumerate(self.platforms):
            for j, refinery in enumerate(self.refineries):
                distance = self._calculate_distance(platform['location'], refinery['location'])
                # Ship costs are more linear with distance but have fixed loading costs
                fixed_cost = 500  # Fixed loading/unloading cost per route
                variable_cost = base_ship_cost * distance
                costs[i, j] = fixed_cost / 1000 + variable_cost  # Normalize fixed cost
                
        return costs
    
    def optimize_routes(self):
        """Solve the transportation optimization problem using linear programming"""
        # Create cost vector for the linear program
        # Variables: [x_11, x_12, ..., x_mn, y_11, y_12, ..., y_mn]
        pipeline_costs_flat = self.pipeline_costs.flatten()
        ship_costs_flat = self.ship_costs.flatten()
        c = np.concatenate([pipeline_costs_flat, ship_costs_flat])
        
        # Replace infinite costs with a large number
        c[c == float('inf')] = 1e6
        
        # Constraint matrices
        A_eq = []
        b_eq = []
        
        # Supply constraints: sum over j of (x_ij + y_ij) <= S_i
        for i in range(self.m):
            constraint = np.zeros(2 * self.m * self.n)
            for j in range(self.n):
                # Pipeline variable x_ij
                constraint[i * self.n + j] = 1
                # Ship variable y_ij
                constraint[self.m * self.n + i * self.n + j] = 1
            A_eq.append(constraint)
            b_eq.append(self.supply[i])
        
        # Demand constraints: sum over i of (x_ij + y_ij) >= D_j
        # Convert to equality by adding slack variables would be complex,
        # so we'll use inequality constraints with bounds
        
        A_ub = []
        b_ub = []
        
        # Convert supply constraints to inequality (≤)
        for i in range(self.m):
            constraint = np.zeros(2 * self.m * self.n)
            for j in range(self.n):
                constraint[i * self.n + j] = 1
                constraint[self.m * self.n + i * self.n + j] = 1
            A_ub.append(constraint)
            b_ub.append(self.supply[i])
        
        # Demand constraints: -sum over i of (x_ij + y_ij) ≤ -D_j
        for j in range(self.n):
            constraint = np.zeros(2 * self.m * self.n)
            for i in range(self.m):
                constraint[i * self.n + j] = -1
                constraint[self.m * self.n + i * self.n + j] = -1
            A_ub.append(constraint)
            b_ub.append(-self.demand[j])
        
        # Bounds: all variables >= 0
        bounds = [(0, None) for _ in range(2 * self.m * self.n)]
        
        # Solve the linear program
        result = linprog(c, A_ub=np.array(A_ub), b_ub=np.array(b_ub), 
                        bounds=bounds, method='highs')
        
        if result.success:
            # Extract solution
            solution = result.x
            pipeline_solution = solution[:self.m * self.n].reshape((self.m, self.n))
            ship_solution = solution[self.m * self.n:].reshape((self.m, self.n))
            
            return {
                'pipeline_flows': pipeline_solution,
                'ship_flows': ship_solution,
                'total_cost': result.fun,
                'success': True
            }
        else:
            return {'success': False, 'message': result.message}

def create_visualization(optimizer, solution):
    """Create comprehensive visualizations of the optimization results"""
    
    fig, axes = plt.subplots(2, 2, figsize=(16, 12))
    fig.suptitle('Pipeline and Maritime Transport Route Optimization Results', 
                 fontsize=16, fontweight='bold')
    
    # 1. Geographic representation of routes
    ax1 = axes[0, 0]
    
    # Plot platforms and refineries
    platform_locs = np.array([p['location'] for p in optimizer.platforms])
    refinery_locs = np.array([r['location'] for r in optimizer.refineries])
    
    ax1.scatter(platform_locs[:, 0], platform_locs[:, 1], 
               c='red', s=200, marker='o', label='Oil Platforms', alpha=0.8)
    ax1.scatter(refinery_locs[:, 0], refinery_locs[:, 1], 
               c='blue', s=200, marker='s', label='Refineries', alpha=0.8)
    
    # Draw routes with thickness proportional to flow
    pipeline_flows = solution['pipeline_flows']
    ship_flows = solution['ship_flows']
    
    max_flow = max(np.max(pipeline_flows), np.max(ship_flows))
    
    for i in range(optimizer.m):
        for j in range(optimizer.n):
            if pipeline_flows[i, j] > 0.1:  # Only show significant flows
                thickness = (pipeline_flows[i, j] / max_flow) * 5
                ax1.plot([platform_locs[i, 0], refinery_locs[j, 0]], 
                        [platform_locs[i, 1], refinery_locs[j, 1]], 
                        'g-', linewidth=thickness, alpha=0.7, label='Pipeline' if i==0 and j==0 else "")
            
            if ship_flows[i, j] > 0.1:
                thickness = (ship_flows[i, j] / max_flow) * 5
                ax1.plot([platform_locs[i, 0], refinery_locs[j, 0]], 
                        [platform_locs[i, 1], refinery_locs[j, 1]], 
                        'b--', linewidth=thickness, alpha=0.7, label='Maritime' if i==0 and j==1 else "")
    
    ax1.set_xlabel('X Coordinate (km)')
    ax1.set_ylabel('Y Coordinate (km)')
    ax1.set_title('Geographic Route Visualization')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. Cost comparison heatmap
    ax2 = axes[0, 1]
    cost_comparison = optimizer.pipeline_costs + optimizer.ship_costs
    # Replace inf values for visualization
    cost_comparison[cost_comparison > 1000] = np.nan
    
    im2 = ax2.imshow(cost_comparison, cmap='YlOrRd', aspect='auto')
    ax2.set_xlabel('Refinery Index')
    ax2.set_ylabel('Platform Index')
    ax2.set_title('Total Transportation Cost Matrix\n($/barrel)')
    
    # Add colorbar
    plt.colorbar(im2, ax=ax2)
    
    # Add text annotations
    for i in range(optimizer.m):
        for j in range(optimizer.n):
            if not np.isnan(cost_comparison[i, j]):
                ax2.text(j, i, f'{cost_comparison[i, j]:.0f}', 
                        ha='center', va='center', fontsize=8)
    
    # 3. Flow distribution
    ax3 = axes[1, 0]
    
    # Create stacked bar chart for each platform-refinery pair
    x_labels = [f'P{i+1}→R{j+1}' for i in range(optimizer.m) for j in range(optimizer.n) 
                if pipeline_flows[i,j] > 0.1 or ship_flows[i,j] > 0.1]
    pipeline_values = [pipeline_flows[i,j] for i in range(optimizer.m) for j in range(optimizer.n) 
                      if pipeline_flows[i,j] > 0.1 or ship_flows[i,j] > 0.1]
    ship_values = [ship_flows[i,j] for i in range(optimizer.m) for j in range(optimizer.n) 
                  if pipeline_flows[i,j] > 0.1 or ship_flows[i,j] > 0.1]
    
    if x_labels:  # Only plot if there are active routes
        x_pos = np.arange(len(x_labels))
        ax3.bar(x_pos, pipeline_values, label='Pipeline', alpha=0.8, color='green')
        ax3.bar(x_pos, ship_values, bottom=pipeline_values, label='Maritime', alpha=0.8, color='blue')
        
        ax3.set_xlabel('Route (Platform → Refinery)')
        ax3.set_ylabel('Flow (1000 barrels/day)')
        ax3.set_title('Transportation Flow by Route and Mode')
        ax3.set_xticks(x_pos)
        ax3.set_xticklabels(x_labels, rotation=45)
        ax3.legend()
    
    # 4. Supply and demand analysis
    ax4 = axes[1, 1]
    
    platform_names = [f'Platform {i+1}' for i in range(optimizer.m)]
    refinery_names = [f'Refinery {j+1}' for j in range(optimizer.n)]
    
    total_supply = optimizer.supply
    total_outflow = np.sum(pipeline_flows + ship_flows, axis=1)
    total_demand = optimizer.demand
    total_inflow = np.sum(pipeline_flows + ship_flows, axis=0)
    
    x_supply = np.arange(len(platform_names))
    x_demand = np.arange(len(refinery_names)) + len(platform_names) + 1
    
    ax4.bar(x_supply, total_supply, alpha=0.7, label='Available Supply', color='orange')
    ax4.bar(x_supply, total_outflow, alpha=0.7, label='Actual Outflow', color='red')
    ax4.bar(x_demand, total_demand, alpha=0.7, label='Required Demand', color='lightblue')
    ax4.bar(x_demand, total_inflow, alpha=0.7, label='Actual Inflow', color='navy')
    
    # Combine labels
    all_labels = platform_names + [''] + refinery_names
    all_x = np.concatenate([x_supply, [len(platform_names)], x_demand])
    
    ax4.set_xticks(all_x)
    ax4.set_xticklabels(all_labels, rotation=45)
    ax4.set_ylabel('Volume (1000 barrels/day)')
    ax4.set_title('Supply and Demand Balance Analysis')
    ax4.legend()
    
    plt.tight_layout()
    plt.show()
    
    return fig

def print_detailed_results(optimizer, solution):
    """Print detailed analysis of the optimization results"""
    
    print("="*80)
    print("TRANSPORTATION OPTIMIZATION RESULTS")
    print("="*80)
    
    pipeline_flows = solution['pipeline_flows']
    ship_flows = solution['ship_flows']
    total_flows = pipeline_flows + ship_flows
    
    print(f"\nTotal Minimum Cost: ${solution['total_cost']:.2f}")
    print(f"Average Cost per Barrel: ${solution['total_cost']/np.sum(total_flows):.2f}")
    
    print("\nFLOW ALLOCATION SUMMARY:")
    print("-" * 50)
    
    total_pipeline = np.sum(pipeline_flows)
    total_ship = np.sum(ship_flows)
    total_volume = total_pipeline + total_ship
    
    print(f"Total Pipeline Transport: {total_pipeline:.1f} (1000 barrels/day)")
    print(f"Total Maritime Transport: {total_ship:.1f} (1000 barrels/day)")
    print(f"Pipeline Share: {(total_pipeline/total_volume)*100:.1f}%")
    print(f"Maritime Share: {(total_ship/total_volume)*100:.1f}%")
    
    print("\nROUTE-SPECIFIC ANALYSIS:")
    print("-" * 50)
    
    for i in range(optimizer.m):
        for j in range(optimizer.n):
            if total_flows[i, j] > 0.1:
                pipeline_cost = optimizer.pipeline_costs[i, j] * pipeline_flows[i, j]
                ship_cost = optimizer.ship_costs[i, j] * ship_flows[i, j]
                total_cost = pipeline_cost + ship_cost
                
                print(f"\nPlatform {i+1} → Refinery {j+1}:")
                print(f"  Pipeline: {pipeline_flows[i, j]:.1f} units @ ${optimizer.pipeline_costs[i, j]:.2f}/unit = ${pipeline_cost:.2f}")
                print(f"  Maritime: {ship_flows[i, j]:.1f} units @ ${optimizer.ship_costs[i, j]:.2f}/unit = ${ship_cost:.2f}")
                print(f"  Total Route Cost: ${total_cost:.2f}")
    
    print("\nCAPACITY UTILIZATION:")
    print("-" * 50)
    
    for i in range(optimizer.m):
        utilization = np.sum(total_flows[i, :]) / optimizer.supply[i] * 100
        print(f"Platform {i+1}: {utilization:.1f}% ({np.sum(total_flows[i, :]):.1f}/{optimizer.supply[i]:.1f})")
    
    print("\nDEMAND SATISFACTION:")
    print("-" * 50)
    
    for j in range(optimizer.n):
        satisfaction = np.sum(total_flows[:, j]) / optimizer.demand[j] * 100
        print(f"Refinery {j+1}: {satisfaction:.1f}% ({np.sum(total_flows[:, j]):.1f}/{optimizer.demand[j]:.1f})")

# Main execution
if __name__ == "__main__":
    # Define oil platforms with locations and capacities
    platforms = [
        {'location': (50, 100), 'capacity': 1500},   # Platform 1
        {'location': (150, 200), 'capacity': 2000},  # Platform 2
        {'location': (250, 150), 'capacity': 1200},  # Platform 3
        {'location': (100, 300), 'capacity': 1800},  # Platform 4
    ]
    
    # Define refineries with locations and demand
    refineries = [
        {'location': (200, 50), 'demand': 1800},     # Refinery 1
        {'location': (300, 250), 'demand': 1500},    # Refinery 2
        {'location': (400, 100), 'demand': 1200},    # Refinery 3
        {'location': (350, 300), 'demand': 1000},    # Refinery 4
    ]
    
    # Create optimizer instance
    optimizer = TransportOptimizer(platforms, refineries)
    
    # Display problem setup
    print("TRANSPORTATION OPTIMIZATION PROBLEM SETUP")
    print("="*60)
    print(f"Number of Oil Platforms: {len(platforms)}")
    print(f"Number of Refineries: {len(refineries)}")
    print(f"Total Supply Capacity: {sum(p['capacity'] for p in platforms):,} (1000 barrels/day)")
    print(f"Total Demand: {sum(r['demand'] for r in refineries):,} (1000 barrels/day)")
    
    # Solve optimization problem
    print("\nSolving optimization problem...")
    solution = optimizer.optimize_routes()
    
    if solution['success']:
        print("✓ Optimization completed successfully!")
        
        # Print detailed results
        print_detailed_results(optimizer, solution)
        
        # Create visualizations
        fig = create_visualization(optimizer, solution)
        
    else:
        print("✗ Optimization failed:", solution['message'])