import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
import seaborn as sns

# Set style for better-looking plots
plt.style.use('default')
sns.set_palette("husl")

# Define parameters for our example
# Alpha: cost coefficient for parental survival (higher = more costly)
# Beta: benefit coefficient for offspring survival (higher = more beneficial)
alpha = 0.3  # Moderate cost to parent survival
beta = 0.5   # Moderate benefit to offspring survival

# Define the investment range
investment = np.linspace(0, 10, 1000)

def parent_survival(I, alpha):
    """
    Parent survival probability as a function of reproductive investment.
    Uses exponential decay model: higher investment reduces parent survival.
    
    Parameters:
    I (float or array): Reproductive investment level
    alpha (float): Cost coefficient for parental survival
    
    Returns:
    float or array: Parent survival probability
    """
    return np.exp(-alpha * I)

def offspring_survival(I, beta):
    """
    Offspring survival probability as a function of reproductive investment.
    Uses saturating exponential model: higher investment increases offspring survival
    with diminishing returns.
    
    Parameters:
    I (float or array): Reproductive investment level
    beta (float): Benefit coefficient for offspring survival
    
    Returns:
    float or array: Offspring survival probability
    """
    return 1 - np.exp(-beta * I)

def fitness_function(I, alpha, beta):
    """
    Overall fitness function combining parent and offspring survival.
    Fitness is the product of parent survival and offspring survival probabilities.
    
    Parameters:
    I (float or array): Reproductive investment level
    alpha (float): Cost coefficient for parental survival
    beta (float): Benefit coefficient for offspring survival
    
    Returns:
    float or array: Overall fitness value
    """
    return parent_survival(I, alpha) * offspring_survival(I, beta)

# Calculate survival probabilities and fitness
parent_surv = parent_survival(investment, alpha)
offspring_surv = offspring_survival(investment, beta)
fitness = fitness_function(investment, alpha, beta)

# Find optimal investment level
def negative_fitness(I):
    """Helper function for optimization (minimize negative fitness = maximize fitness)"""
    return -fitness_function(I, alpha, beta)

# Find the optimal investment level
result = minimize_scalar(negative_fitness, bounds=(0, 10), method='bounded')
optimal_investment = result.x
optimal_fitness = -result.fun

print(f"Optimal reproductive investment: {optimal_investment:.3f}")
print(f"Maximum fitness: {optimal_fitness:.3f}")
print(f"Parent survival at optimum: {parent_survival(optimal_investment, alpha):.3f}")
print(f"Offspring survival at optimum: {offspring_survival(optimal_investment, beta):.3f}")

# Create comprehensive visualization
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('Reproductive Investment Optimization Analysis', fontsize=16, fontweight='bold')

# Plot 1: Parent survival probability
ax1.plot(investment, parent_surv, 'r-', linewidth=3, label='Parent Survival')
ax1.axvline(optimal_investment, color='black', linestyle='--', alpha=0.7, 
            label=f'Optimal Investment = {optimal_investment:.2f}')
ax1.set_xlabel('Reproductive Investment (I)')
ax1.set_ylabel('Parent Survival Probability')
ax1.set_title('Parent Survival: $S_p(I) = e^{-\\alpha I}$')
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.set_ylim(0, 1)

# Plot 2: Offspring survival probability
ax2.plot(investment, offspring_surv, 'b-', linewidth=3, label='Offspring Survival')
ax2.axvline(optimal_investment, color='black', linestyle='--', alpha=0.7,
            label=f'Optimal Investment = {optimal_investment:.2f}')
ax2.set_xlabel('Reproductive Investment (I)')
ax2.set_ylabel('Offspring Survival Probability')
ax2.set_title('Offspring Survival: $S_o(I) = 1 - e^{-\\beta I}$')
ax2.grid(True, alpha=0.3)
ax2.legend()
ax2.set_ylim(0, 1)

# Plot 3: Combined fitness function
ax3.plot(investment, fitness, 'g-', linewidth=3, label='Fitness Function')
ax3.plot(optimal_investment, optimal_fitness, 'ro', markersize=10, 
         label=f'Maximum at I = {optimal_investment:.2f}')
ax3.axvline(optimal_investment, color='black', linestyle='--', alpha=0.7)
ax3.set_xlabel('Reproductive Investment (I)')
ax3.set_ylabel('Fitness: $S_p(I) \\times S_o(I)$')
ax3.set_title('Overall Fitness Function')
ax3.grid(True, alpha=0.3)
ax3.legend()

# Plot 4: Trade-off visualization
ax4.plot(investment, parent_surv, 'r-', linewidth=2, label='Parent Survival', alpha=0.8)
ax4.plot(investment, offspring_surv, 'b-', linewidth=2, label='Offspring Survival', alpha=0.8)
ax4.plot(investment, fitness, 'g-', linewidth=3, label='Combined Fitness')
ax4.axvline(optimal_investment, color='black', linestyle='--', alpha=0.7,
            label=f'Optimal Point')
ax4.set_xlabel('Reproductive Investment (I)')
ax4.set_ylabel('Probability')
ax4.set_title('Trade-off Between Parent and Offspring Survival')
ax4.grid(True, alpha=0.3)
ax4.legend()
ax4.set_ylim(0, 1)

plt.tight_layout()
plt.show()

# Additional analysis: Parameter sensitivity
print("\n" + "="*60)
print("PARAMETER SENSITIVITY ANALYSIS")
print("="*60)

# Test different alpha and beta values
alpha_values = [0.1, 0.3, 0.5, 0.7]
beta_values = [0.2, 0.4, 0.6, 0.8]

fig2, ((ax5, ax6), (ax7, ax8)) = plt.subplots(2, 2, figsize=(15, 12))
fig2.suptitle('Parameter Sensitivity Analysis', fontsize=16, fontweight='bold')

# Vary alpha (parent cost coefficient)
colors_alpha = ['lightcoral', 'red', 'darkred', 'maroon']
for i, a in enumerate(alpha_values):
    fitness_alpha = fitness_function(investment, a, beta)
    ax5.plot(investment, fitness_alpha, color=colors_alpha[i], linewidth=2,
             label=f'α = {a}')
    
    # Find optimum for this alpha
    result_alpha = minimize_scalar(lambda x: -fitness_function(x, a, beta), 
                                   bounds=(0, 10), method='bounded')
    ax5.plot(result_alpha.x, -result_alpha.fun, 'o', color=colors_alpha[i], markersize=6)

ax5.set_xlabel('Reproductive Investment (I)')
ax5.set_ylabel('Fitness')
ax5.set_title('Effect of Parent Cost Coefficient (α)')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Vary beta (offspring benefit coefficient)
colors_beta = ['lightblue', 'blue', 'darkblue', 'navy']
for i, b in enumerate(beta_values):
    fitness_beta = fitness_function(investment, alpha, b)
    ax6.plot(investment, fitness_beta, color=colors_beta[i], linewidth=2,
             label=f'β = {b}')
    
    # Find optimum for this beta
    result_beta = minimize_scalar(lambda x: -fitness_function(x, alpha, b), 
                                  bounds=(0, 10), method='bounded')
    ax6.plot(result_beta.x, -result_beta.fun, 'o', color=colors_beta[i], markersize=6)

ax6.set_xlabel('Reproductive Investment (I)')
ax6.set_ylabel('Fitness')
ax6.set_title('Effect of Offspring Benefit Coefficient (β)')
ax6.legend()
ax6.grid(True, alpha=0.3)

# Optimal investment heatmap
alpha_range = np.linspace(0.1, 1.0, 20)
beta_range = np.linspace(0.1, 1.0, 20)
optimal_investments = np.zeros((len(alpha_range), len(beta_range)))

for i, a in enumerate(alpha_range):
    for j, b in enumerate(beta_range):
        result_ij = minimize_scalar(lambda x: -fitness_function(x, a, b), 
                                    bounds=(0, 10), method='bounded')
        optimal_investments[i, j] = result_ij.x

im = ax7.imshow(optimal_investments, extent=[beta_range[0], beta_range[-1], 
                                             alpha_range[0], alpha_range[-1]], 
                aspect='auto', origin='lower', cmap='viridis')
ax7.set_xlabel('Offspring Benefit Coefficient (β)')
ax7.set_ylabel('Parent Cost Coefficient (α)')
ax7.set_title('Optimal Investment Landscape')
plt.colorbar(im, ax=ax7, label='Optimal Investment Level')

# Fitness landscape
max_fitness_values = np.zeros((len(alpha_range), len(beta_range)))
for i, a in enumerate(alpha_range):
    for j, b in enumerate(beta_range):
        result_ij = minimize_scalar(lambda x: -fitness_function(x, a, b), 
                                    bounds=(0, 10), method='bounded')
        max_fitness_values[i, j] = -result_ij.fun

im2 = ax8.imshow(max_fitness_values, extent=[beta_range[0], beta_range[-1], 
                                             alpha_range[0], alpha_range[-1]], 
                 aspect='auto', origin='lower', cmap='plasma')
ax8.set_xlabel('Offspring Benefit Coefficient (β)')
ax8.set_ylabel('Parent Cost Coefficient (α)')
ax8.set_title('Maximum Fitness Landscape')
plt.colorbar(im2, ax=ax8, label='Maximum Fitness')

plt.tight_layout()
plt.show()

# Summary statistics
print(f"\nSUMMARY FOR BASE PARAMETERS (α={alpha}, β={beta}):")
print(f"Optimal investment level: {optimal_investment:.3f}")
print(f"Maximum fitness achieved: {optimal_fitness:.3f}")
print(f"Parent survival at optimum: {parent_survival(optimal_investment, alpha):.3f}")
print(f"Offspring survival at optimum: {offspring_survival(optimal_investment, beta):.3f}")
print(f"Investment efficiency: {optimal_fitness/optimal_investment:.3f} fitness per unit investment")

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import pandas as pd
from typing import List, Tuple
import math

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

class MateSelectionOptimizer:
    """
    A class to simulate and optimize mate selection strategies using
    the Secretary Problem framework applied to evolutionary biology.
    """
    
    def __init__(self, population_size: int = 100):
        """
        Initialize the optimizer with a given population size.
        
        Parameters:
        -----------
        population_size : int
            Total number of potential mates in the population
        """
        self.population_size = population_size
        self.optimal_threshold = int(population_size / math.e)
        
    def generate_population(self, distribution: str = 'normal') -> np.ndarray:
        """
        Generate a population with genetic fitness scores.
        
        Parameters:
        -----------
        distribution : str
            Distribution type for fitness scores ('normal', 'exponential', 'uniform')
            
        Returns:
        --------
        np.ndarray : Array of fitness scores
        """
        np.random.seed(42)  # For reproducibility
        
        if distribution == 'normal':
            # Normal distribution: μ=50, σ=15 (IQ-like scale)
            fitness_scores = np.random.normal(50, 15, self.population_size)
        elif distribution == 'exponential':
            # Exponential distribution (right-skewed)
            fitness_scores = np.random.exponential(30, self.population_size)
        elif distribution == 'uniform':
            # Uniform distribution
            fitness_scores = np.random.uniform(0, 100, self.population_size)
        else:
            raise ValueError("Distribution must be 'normal', 'exponential', or 'uniform'")
            
        return np.maximum(fitness_scores, 0)  # Ensure no negative fitness
    
    def simulate_random_strategy(self, fitness_scores: np.ndarray, 
                                num_simulations: int = 1000) -> List[float]:
        """
        Simulate random mate selection strategy.
        
        Parameters:
        -----------
        fitness_scores : np.ndarray
            Population fitness scores
        num_simulations : int
            Number of simulation runs
            
        Returns:
        --------
        List[float] : Selected fitness scores from each simulation
        """
        selected_scores = []
        
        for _ in range(num_simulations):
            # Randomly shuffle the population
            shuffled_scores = np.random.permutation(fitness_scores)
            # Randomly select a position to stop
            stop_position = np.random.randint(1, len(shuffled_scores) + 1)
            selected_scores.append(shuffled_scores[stop_position - 1])
            
        return selected_scores
    
    def simulate_optimal_strategy(self, fitness_scores: np.ndarray, 
                                 num_simulations: int = 1000) -> List[float]:
        """
        Simulate optimal mate selection strategy (Secretary Problem solution).
        
        Parameters:
        -----------
        fitness_scores : np.ndarray
            Population fitness scores
        num_simulations : int
            Number of simulation runs
            
        Returns:
        --------
        List[float] : Selected fitness scores from each simulation
        """
        selected_scores = []
        
        for _ in range(num_simulations):
            # Randomly shuffle the population
            shuffled_scores = np.random.permutation(fitness_scores)
            
            # Observe first n/e candidates
            observation_phase = shuffled_scores[:self.optimal_threshold]
            if len(observation_phase) == 0:
                # If threshold is 0, select the first candidate
                selected_scores.append(shuffled_scores[0])
                continue
                
            # Find the maximum in observation phase
            observation_max = np.max(observation_phase)
            
            # Selection phase: choose first candidate better than observation max
            selected = None
            for i in range(self.optimal_threshold, len(shuffled_scores)):
                if shuffled_scores[i] > observation_max:
                    selected = shuffled_scores[i]
                    break
            
            # If no candidate is better, select the last one
            if selected is None:
                selected = shuffled_scores[-1]
                
            selected_scores.append(selected)
            
        return selected_scores
    
    def simulate_threshold_strategy(self, fitness_scores: np.ndarray, 
                                   threshold_percentile: float = 80,
                                   num_simulations: int = 1000) -> List[float]:
        """
        Simulate threshold-based strategy (accept first candidate above percentile).
        
        Parameters:
        -----------
        fitness_scores : np.ndarray
            Population fitness scores
        threshold_percentile : float
            Percentile threshold (0-100)
        num_simulations : int
            Number of simulation runs
            
        Returns:
        --------
        List[float] : Selected fitness scores from each simulation
        """
        selected_scores = []
        threshold_value = np.percentile(fitness_scores, threshold_percentile)
        
        for _ in range(num_simulations):
            # Randomly shuffle the population
            shuffled_scores = np.random.permutation(fitness_scores)
            
            # Select first candidate above threshold
            selected = None
            for score in shuffled_scores:
                if score >= threshold_value:
                    selected = score
                    break
            
            # If no candidate meets threshold, select the last one
            if selected is None:
                selected = shuffled_scores[-1]
                
            selected_scores.append(selected)
            
        return selected_scores
    
    def calculate_success_probability(self, selected_scores: List[float], 
                                    true_maximum: float) -> float:
        """
        Calculate the probability of selecting the true best mate.
        
        Parameters:
        -----------
        selected_scores : List[float]
            Scores selected by the strategy
        true_maximum : float
            True maximum fitness in the population
            
        Returns:
        --------
        float : Success probability
        """
        successes = sum(1 for score in selected_scores 
                       if abs(score - true_maximum) < 1e-10)
        return successes / len(selected_scores)

# Initialize the optimizer
optimizer = MateSelectionOptimizer(population_size=100)

# Generate population with normal distribution
population_fitness = optimizer.generate_population('normal')
true_best_fitness = np.max(population_fitness)

print(f"Population size: {len(population_fitness)}")
print(f"Optimal observation threshold: {optimizer.optimal_threshold}")
print(f"True best fitness score: {true_best_fitness:.2f}")
print(f"Population mean fitness: {np.mean(population_fitness):.2f}")
print(f"Population std fitness: {np.std(population_fitness):.2f}")

# Run simulations with different strategies
num_sims = 10000

print("\n" + "="*50)
print("Running simulations...")
print("="*50)

# Random strategy
random_results = optimizer.simulate_random_strategy(population_fitness, num_sims)
random_success_prob = optimizer.calculate_success_probability(random_results, true_best_fitness)

# Optimal strategy (Secretary Problem)
optimal_results = optimizer.simulate_optimal_strategy(population_fitness, num_sims)
optimal_success_prob = optimizer.calculate_success_probability(optimal_results, true_best_fitness)

# Threshold strategies with different percentiles
threshold_70_results = optimizer.simulate_threshold_strategy(population_fitness, 70, num_sims)
threshold_70_success = optimizer.calculate_success_probability(threshold_70_results, true_best_fitness)

threshold_80_results = optimizer.simulate_threshold_strategy(population_fitness, 80, num_sims)
threshold_80_success = optimizer.calculate_success_probability(threshold_80_results, true_best_fitness)

threshold_90_results = optimizer.simulate_threshold_strategy(population_fitness, 90, num_sims)
threshold_90_success = optimizer.calculate_success_probability(threshold_90_results, true_best_fitness)

# Print results
print(f"\nRandom Strategy:")
print(f"  Average selected fitness: {np.mean(random_results):.2f}")
print(f"  Success probability: {random_success_prob:.4f}")

print(f"\nOptimal Strategy (Secretary Problem):")
print(f"  Average selected fitness: {np.mean(optimal_results):.2f}")
print(f"  Success probability: {optimal_success_prob:.4f}")
print(f"  Theoretical success probability: {1/math.e:.4f}")

print(f"\nThreshold Strategy (70th percentile):")
print(f"  Average selected fitness: {np.mean(threshold_70_results):.2f}")
print(f"  Success probability: {threshold_70_success:.4f}")

print(f"\nThreshold Strategy (80th percentile):")
print(f"  Average selected fitness: {np.mean(threshold_80_results):.2f}")
print(f"  Success probability: {threshold_80_success:.4f}")

print(f"\nThreshold Strategy (90th percentile):")
print(f"  Average selected fitness: {np.mean(threshold_90_results):.2f}")
print(f"  Success probability: {threshold_90_success:.4f}")

# Create comprehensive visualizations
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Mate Selection Optimization: Comparison of Strategies', fontsize=16, fontweight='bold')

# 1. Population distribution
axes[0, 0].hist(population_fitness, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
axes[0, 0].axvline(true_best_fitness, color='red', linestyle='--', linewidth=2, label=f'True Best: {true_best_fitness:.1f}')
axes[0, 0].axvline(np.mean(population_fitness), color='green', linestyle='--', linewidth=2, label=f'Mean: {np.mean(population_fitness):.1f}')
axes[0, 0].set_title('Population Fitness Distribution')
axes[0, 0].set_xlabel('Fitness Score')
axes[0, 0].set_ylabel('Frequency')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. Strategy comparison - average scores
strategies = ['Random', 'Optimal\n(Secretary)', 'Threshold\n(70%)', 'Threshold\n(80%)', 'Threshold\n(90%)']
avg_scores = [np.mean(random_results), np.mean(optimal_results), 
              np.mean(threshold_70_results), np.mean(threshold_80_results), np.mean(threshold_90_results)]
colors = ['gray', 'gold', 'lightcoral', 'lightgreen', 'lightblue']

bars = axes[0, 1].bar(strategies, avg_scores, color=colors, edgecolor='black', linewidth=1)
axes[0, 1].axhline(true_best_fitness, color='red', linestyle='--', linewidth=2, label=f'True Best: {true_best_fitness:.1f}')
axes[0, 1].axhline(np.mean(population_fitness), color='green', linestyle='--', linewidth=2, label=f'Population Mean: {np.mean(population_fitness):.1f}')
axes[0, 1].set_title('Average Selected Fitness by Strategy')
axes[0, 1].set_ylabel('Average Fitness Score')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Add value labels on bars
for i, (bar, score) in enumerate(zip(bars, avg_scores)):
    axes[0, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
                    f'{score:.1f}', ha='center', va='bottom', fontweight='bold')

# 3. Success probability comparison
success_probs = [random_success_prob, optimal_success_prob, threshold_70_success, 
                threshold_80_success, threshold_90_success]

bars2 = axes[0, 2].bar(strategies, success_probs, color=colors, edgecolor='black', linewidth=1)
axes[0, 2].axhline(1/math.e, color='red', linestyle='--', linewidth=2, label=f'Theoretical Optimal: {1/math.e:.3f}')
axes[0, 2].set_title('Probability of Selecting the Best Mate')
axes[0, 2].set_ylabel('Success Probability')
axes[0, 2].set_ylim(0, max(success_probs) * 1.2)
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)

# Add value labels on bars
for i, (bar, prob) in enumerate(zip(bars2, success_probs)):
    axes[0, 2].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005, 
                    f'{prob:.3f}', ha='center', va='bottom', fontweight='bold')

# 4. Distribution comparison of selected scores
data_for_plot = [random_results, optimal_results, threshold_80_results]
labels_for_plot = ['Random', 'Optimal (Secretary)', 'Threshold (80%)']
colors_for_plot = ['gray', 'gold', 'lightgreen']

for i, (data, label, color) in enumerate(zip(data_for_plot, labels_for_plot, colors_for_plot)):
    axes[1, 0].hist(data, bins=30, alpha=0.6, label=label, color=color, density=True)

axes[1, 0].axvline(true_best_fitness, color='red', linestyle='--', linewidth=2, label=f'True Best: {true_best_fitness:.1f}')
axes[1, 0].set_title('Distribution of Selected Fitness Scores')
axes[1, 0].set_xlabel('Selected Fitness Score')
axes[1, 0].set_ylabel('Density')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 5. Box plot comparison
box_data = [random_results, optimal_results, threshold_70_results, threshold_80_results, threshold_90_results]
bp = axes[1, 1].boxplot(box_data, labels=strategies, patch_artist=True)

for patch, color in zip(bp['boxes'], colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.7)

axes[1, 1].axhline(true_best_fitness, color='red', linestyle='--', linewidth=2, label=f'True Best: {true_best_fitness:.1f}')
axes[1, 1].set_title('Distribution Comparison (Box Plots)')
axes[1, 1].set_ylabel('Selected Fitness Score')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# 6. Performance vs Population Size Analysis
population_sizes = range(10, 201, 10)
optimal_performance = []
random_performance = []

for pop_size in population_sizes:
    temp_optimizer = MateSelectionOptimizer(pop_size)
    temp_population = temp_optimizer.generate_population('normal')
    temp_true_best = np.max(temp_population)
    
    # Run smaller simulations for speed
    temp_optimal = temp_optimizer.simulate_optimal_strategy(temp_population, 1000)
    temp_random = temp_optimizer.simulate_random_strategy(temp_population, 1000)
    
    optimal_performance.append(np.mean(temp_optimal))
    random_performance.append(np.mean(temp_random))

axes[1, 2].plot(population_sizes, optimal_performance, 'o-', color='gold', linewidth=2, 
                markersize=4, label='Optimal Strategy')
axes[1, 2].plot(population_sizes, random_performance, 's-', color='gray', linewidth=2, 
                markersize=4, label='Random Strategy')
axes[1, 2].set_title('Strategy Performance vs Population Size')
axes[1, 2].set_xlabel('Population Size')
axes[1, 2].set_ylabel('Average Selected Fitness')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Effect of different observation thresholds
print("\n" + "="*60)
print("ANALYSIS: Effect of Different Observation Thresholds")
print("="*60)

thresholds = range(1, min(50, optimizer.population_size), 2)
threshold_performance = []

for threshold in thresholds:
    # Temporary modification of the optimizer threshold
    original_threshold = optimizer.optimal_threshold
    optimizer.optimal_threshold = threshold
    
    temp_results = optimizer.simulate_optimal_strategy(population_fitness, 2000)
    temp_success = optimizer.calculate_success_probability(temp_results, true_best_fitness)
    threshold_performance.append(temp_success)
    
    # Restore original threshold
    optimizer.optimal_threshold = original_threshold

# Plot threshold analysis
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(thresholds, threshold_performance, 'o-', color='blue', linewidth=2, markersize=6)
plt.axvline(optimizer.population_size / math.e, color='red', linestyle='--', linewidth=2, 
            label=f'Theoretical Optimal: {optimizer.population_size / math.e:.1f}')
plt.axhline(1/math.e, color='red', linestyle=':', linewidth=2, 
            label=f'Theoretical Success Rate: {1/math.e:.3f}')
plt.title('Success Rate vs Observation Threshold')
plt.xlabel('Observation Threshold (Number of Candidates)')
plt.ylabel('Success Probability')
plt.legend()
plt.grid(True, alpha=0.3)

# Mathematical formula visualization
plt.subplot(1, 2, 2)
x = np.linspace(0.1, 0.9, 100)
y = -x * np.log(x)  # This is the success probability function p(r) = -r*ln(r)

plt.plot(x, y, 'purple', linewidth=3)
plt.axvline(1/math.e, color='red', linestyle='--', linewidth=2, 
            label=f'Optimal r = 1/e ≈ {1/math.e:.3f}')
plt.axhline(1/math.e, color='red', linestyle=':', linewidth=2, 
            label=f'Maximum Success Rate ≈ {1/math.e:.3f}')
plt.title('Theoretical Success Probability Function\n' + r'$P(r) = -r \ln(r)$')
plt.xlabel('Observation Ratio (r = threshold/population)')
plt.ylabel('Success Probability')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Summary statistics
print(f"\nSUMMARY STATISTICS:")
print(f"="*40)
print(f"Population Statistics:")
print(f"  Size: {len(population_fitness)}")
print(f"  Mean: {np.mean(population_fitness):.2f}")
print(f"  Std Dev: {np.std(population_fitness):.2f}")
print(f"  Best Individual: {true_best_fitness:.2f}")
print(f"  95th Percentile: {np.percentile(population_fitness, 95):.2f}")

print(f"\nStrategy Effectiveness Ranking:")
strategy_results = [
    ("Random", np.mean(random_results), random_success_prob),
    ("Optimal (Secretary)", np.mean(optimal_results), optimal_success_prob),
    ("Threshold (70%)", np.mean(threshold_70_results), threshold_70_success),
    ("Threshold (80%)", np.mean(threshold_80_results), threshold_80_success),
    ("Threshold (90%)", np.mean(threshold_90_results), threshold_90_success)
]

# Sort by success probability
strategy_results.sort(key=lambda x: x[2], reverse=True)

for i, (strategy, avg_score, success_rate) in enumerate(strategy_results, 1):
    print(f"  {i}. {strategy:20s} | Avg Score: {avg_score:6.2f} | Success Rate: {success_rate:.4f}")

print(f"\nKey Insights:")
print(f"  • The optimal strategy achieves {optimal_success_prob:.1%} success rate")
print(f"  • This is {optimal_success_prob/random_success_prob:.1f}x better than random selection")
print(f"  • High threshold strategies risk missing good candidates entirely")
print(f"  • The theoretical maximum success rate is {1/math.e:.1%}")

class MateSelectionOptimizer:
    def __init__(self, population_size: int = 100):
        self.population_size = population_size
        self.optimal_threshold = int(population_size / math.e)

def generate_population(self, distribution: str = 'normal') -> np.ndarray:

def simulate_random_strategy(self, fitness_scores, num_simulations=1000):

def simulate_optimal_strategy(self, fitness_scores, num_simulations=1000):

def simulate_threshold_strategy(self, fitness_scores, threshold_percentile=80, num_simulations=1000):

def calculate_success_probability(self, selected_scores, true_maximum):
    successes = sum(1 for score in selected_scores if abs(score - true_maximum) < 1e-10)
    return successes / len(selected_scores)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
import seaborn as sns
from matplotlib.patches import FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')

# Set up plotting style
plt.style.use('default')
sns.set_palette("husl")

class PredatorPreyModel:
    """
    A model for predator-prey coevolution focusing on hunting vs escape strategies.
    
    This class implements the mathematical framework for studying how predators
    and prey coevolve their strategies in response to each other.
    """
    
    def __init__(self, alpha_p=10.0, beta_p=5.0, alpha_r=8.0, beta_r=4.0, c=0.1):
        """
        Initialize the predator-prey model with fitness parameters.
        
        Parameters:
        - alpha_p: Benefit coefficient for predator success
        - beta_p: Cost coefficient for predator strategy investment
        - alpha_r: Benefit coefficient for prey survival
        - beta_r: Cost coefficient for prey strategy investment
        - c: Environmental constant affecting success probability
        """
        self.alpha_p = alpha_p  # Predator benefit from successful hunting
        self.beta_p = beta_p    # Cost of predator strategy investment
        self.alpha_r = alpha_r  # Prey benefit from escape
        self.beta_r = beta_r    # Cost of prey strategy investment
        self.c = c              # Environmental baseline factor
        
    def success_probability(self, s_p, s_r):
        """
        Calculate the probability of predator success given strategies.
        
        Formula: P(s_p, s_r) = s_p / (s_p + s_r + c)
        
        This represents how predator investment (s_p) competes against
        prey investment (s_r) with environmental factors (c).
        """
        return s_p / (s_p + s_r + self.c)
    
    def predator_fitness(self, s_p, s_r):
        """
        Calculate predator fitness: benefits from hunting success minus costs.
        
        Formula: F_p = α_p * P(s_p, s_r) - β_p * s_p²
        """
        prob = self.success_probability(s_p, s_r)
        return self.alpha_p * prob - self.beta_p * (s_p ** 2)
    
    def prey_fitness(self, s_p, s_r):
        """
        Calculate prey fitness: benefits from escape minus costs.
        
        Formula: F_r = α_r * (1 - P(s_p, s_r)) - β_r * s_r²
        """
        prob = self.success_probability(s_p, s_r)
        return self.alpha_r * (1 - prob) - self.beta_r * (s_r ** 2)
    
    def find_best_response(self, opponent_strategy, is_predator=True, bounds=(0, 2)):
        """
        Find the optimal strategy given the opponent's fixed strategy.
        
        This represents evolutionary adaptation: each species optimizes
        its strategy assuming the other species' strategy is fixed.
        """
        if is_predator:
            # Predator optimizing against fixed prey strategy
            result = minimize_scalar(
                lambda s_p: -self.predator_fitness(s_p, opponent_strategy),
                bounds=bounds,
                method='bounded'
            )
        else:
            # Prey optimizing against fixed predator strategy
            result = minimize_scalar(
                lambda s_r: -self.prey_fitness(opponent_strategy, s_r),
                bounds=bounds,
                method='bounded'
            )
        
        return result.x if result.success else bounds[0]
    
    def coevolutionary_dynamics(self, initial_strategies, generations=100, learning_rate=0.1):
        """
        Simulate the coevolutionary process over multiple generations.
        
        This models how strategies evolve through iterative adaptation,
        where each species responds to the other's current strategy.
        """
        s_p, s_r = initial_strategies
        history = {'predator': [s_p], 'prey': [s_r], 
                  'pred_fitness': [], 'prey_fitness': [], 'success_prob': []}
        
        for generation in range(generations):
            # Calculate current fitness values
            pred_fit = self.predator_fitness(s_p, s_r)
            prey_fit = self.prey_fitness(s_p, s_r)
            success_prob = self.success_probability(s_p, s_r)
            
            history['pred_fitness'].append(pred_fit)
            history['prey_fitness'].append(prey_fit)
            history['success_prob'].append(success_prob)
            
            # Each species adapts toward optimal response
            optimal_s_p = self.find_best_response(s_r, is_predator=True)
            optimal_s_r = self.find_best_response(s_p, is_predator=False)
            
            # Update strategies with learning rate (gradual adaptation)
            s_p += learning_rate * (optimal_s_p - s_p)
            s_r += learning_rate * (optimal_s_r - s_r)
            
            history['predator'].append(s_p)
            history['prey'].append(s_r)
        
        return history
    
    def create_fitness_landscape(self, resolution=50):
        """
        Create a 2D fitness landscape showing how fitness varies with strategies.
        
        This visualization helps understand the strategic interactions
        and identify equilibrium points.
        """
        s_range = np.linspace(0, 2, resolution)
        s_p_grid, s_r_grid = np.meshgrid(s_range, s_range)
        
        pred_fitness_grid = np.zeros_like(s_p_grid)
        prey_fitness_grid = np.zeros_like(s_r_grid)
        success_prob_grid = np.zeros_like(s_p_grid)
        
        for i in range(resolution):
            for j in range(resolution):
                s_p, s_r = s_p_grid[i, j], s_r_grid[i, j]
                pred_fitness_grid[i, j] = self.predator_fitness(s_p, s_r)
                prey_fitness_grid[i, j] = self.prey_fitness(s_p, s_r)
                success_prob_grid[i, j] = self.success_probability(s_p, s_r)
        
        return {
            's_p_grid': s_p_grid, 's_r_grid': s_r_grid,
            'predator_fitness': pred_fitness_grid,
            'prey_fitness': prey_fitness_grid,
            'success_probability': success_prob_grid
        }

# Initialize the model with realistic parameters
model = PredatorPreyModel(alpha_p=10.0, beta_p=5.0, alpha_r=8.0, beta_r=4.0, c=0.1)

# Run coevolutionary simulation
print("🔬 Running Predator-Prey Coevolution Simulation...")
print("=" * 60)

initial_strategies = (0.5, 0.5)  # Both species start with moderate investment
history = model.coevolutionary_dynamics(initial_strategies, generations=100, learning_rate=0.1)

print(f"Initial strategies: Predator = {initial_strategies[0]:.3f}, Prey = {initial_strategies[1]:.3f}")
print(f"Final strategies: Predator = {history['predator'][-1]:.3f}, Prey = {history['prey'][-1]:.3f}")
print(f"Final success probability: {history['success_prob'][-1]:.3f}")

# Generate fitness landscape data
landscape = model.create_fitness_landscape(resolution=40)

# Create comprehensive visualization
fig = plt.figure(figsize=(18, 12))
fig.suptitle('Predator-Prey Coevolutionary Arms Race: Strategic Optimization Analysis', 
             fontsize=16, fontweight='bold', y=0.98)

# 1. Strategy Evolution Over Time
ax1 = plt.subplot(2, 3, 1)
generations = range(len(history['predator']))
plt.plot(generations, history['predator'], 'r-', linewidth=2, label='Predator Strategy', alpha=0.8)
plt.plot(generations, history['prey'], 'b-', linewidth=2, label='Prey Strategy', alpha=0.8)
plt.xlabel('Generation')
plt.ylabel('Strategy Investment')
plt.title('Strategy Evolution Over Time')
plt.legend()
plt.grid(True, alpha=0.3)

# 2. Fitness Evolution
ax2 = plt.subplot(2, 3, 2)
plt.plot(history['pred_fitness'], 'r-', linewidth=2, label='Predator Fitness', alpha=0.8)
plt.plot(history['prey_fitness'], 'b-', linewidth=2, label='Prey Fitness', alpha=0.8)
plt.xlabel('Generation')
plt.ylabel('Fitness')
plt.title('Fitness Evolution')
plt.legend()
plt.grid(True, alpha=0.3)

# 3. Success Probability Over Time
ax3 = plt.subplot(2, 3, 3)
plt.plot(history['success_prob'], 'purple', linewidth=2, alpha=0.8)
plt.xlabel('Generation')
plt.ylabel('Predator Success Probability')
plt.title('Hunting Success Probability')
plt.grid(True, alpha=0.3)

# 4. Predator Fitness Landscape
ax4 = plt.subplot(2, 3, 4)
contour_pred = plt.contourf(landscape['s_p_grid'], landscape['s_r_grid'], 
                           landscape['predator_fitness'], levels=20, cmap='Reds', alpha=0.7)
plt.colorbar(contour_pred, ax=ax4, label='Predator Fitness')
plt.plot(history['predator'], history['prey'], 'k-', linewidth=2, alpha=0.8, label='Evolution Path')
plt.plot(history['predator'][0], history['prey'][0], 'go', markersize=8, label='Start')
plt.plot(history['predator'][-1], history['prey'][-1], 'ro', markersize=8, label='End')
plt.xlabel('Predator Strategy')
plt.ylabel('Prey Strategy')
plt.title('Predator Fitness Landscape')
plt.legend()

# 5. Prey Fitness Landscape
ax5 = plt.subplot(2, 3, 5)
contour_prey = plt.contourf(landscape['s_p_grid'], landscape['s_r_grid'], 
                           landscape['prey_fitness'], levels=20, cmap='Blues', alpha=0.7)
plt.colorbar(contour_prey, ax=ax5, label='Prey Fitness')
plt.plot(history['predator'], history['prey'], 'k-', linewidth=2, alpha=0.8, label='Evolution Path')
plt.plot(history['predator'][0], history['prey'][0], 'go', markersize=8, label='Start')
plt.plot(history['predator'][-1], history['prey'][-1], 'ro', markersize=8, label='End')
plt.xlabel('Predator Strategy')
plt.ylabel('Prey Strategy')
plt.title('Prey Fitness Landscape')
plt.legend()

# 6. Strategy Phase Portrait
ax6 = plt.subplot(2, 3, 6)
plt.plot(history['predator'], history['prey'], 'purple', linewidth=3, alpha=0.8, label='Coevolution Path')
plt.plot(history['predator'][0], history['prey'][0], 'go', markersize=10, label='Initial State')
plt.plot(history['predator'][-1], history['prey'][-1], 'ro', markersize=10, label='Final Equilibrium')

# Add arrows to show direction
n_arrows = 10
indices = np.linspace(0, len(history['predator'])-2, n_arrows, dtype=int)
for i in indices:
    dx = history['predator'][i+1] - history['predator'][i]
    dy = history['prey'][i+1] - history['prey'][i]
    plt.arrow(history['predator'][i], history['prey'][i], dx, dy, 
              head_width=0.02, head_length=0.02, fc='purple', ec='purple', alpha=0.6)

plt.xlabel('Predator Strategy')
plt.ylabel('Prey Strategy')
plt.title('Strategy Phase Portrait')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Detailed Analysis Output
print("\n📊 DETAILED ANALYSIS RESULTS")
print("=" * 60)

# Calculate equilibrium properties
final_pred_strategy = history['predator'][-1]
final_prey_strategy = history['prey'][-1]
final_success_prob = history['success_prob'][-1]
final_pred_fitness = history['pred_fitness'][-1]
final_prey_fitness = history['prey_fitness'][-1]

print(f"🎯 Equilibrium Strategies:")
print(f"   Predator investment: {final_pred_strategy:.4f}")
print(f"   Prey investment: {final_prey_strategy:.4f}")
print(f"   Strategy ratio (P/R): {final_pred_strategy/final_prey_strategy:.4f}")

print(f"\n🏹 Hunting Success:")
print(f"   Success probability: {final_success_prob:.4f} ({final_success_prob*100:.1f}%)")
print(f"   Escape probability: {1-final_success_prob:.4f} ({(1-final_success_prob)*100:.1f}%)")

print(f"\n💪 Final Fitness Values:")
print(f"   Predator fitness: {final_pred_fitness:.4f}")
print(f"   Prey fitness: {final_prey_fitness:.4f}")

# Strategy efficiency analysis
efficiency_pred = final_pred_fitness / final_pred_strategy if final_pred_strategy > 0 else 0
efficiency_prey = final_prey_fitness / final_prey_strategy if final_prey_strategy > 0 else 0

print(f"\n⚡ Strategy Efficiency (Fitness/Investment):")
print(f"   Predator efficiency: {efficiency_pred:.4f}")
print(f"   Prey efficiency: {efficiency_prey:.4f}")

print(f"\n🔄 Coevolutionary Dynamics:")
strategy_change_pred = abs(history['predator'][-1] - history['predator'][0])
strategy_change_prey = abs(history['prey'][-1] - history['prey'][0])
print(f"   Predator strategy change: {strategy_change_pred:.4f}")
print(f"   Prey strategy change: {strategy_change_prey:.4f}")
print(f"   Total evolutionary change: {strategy_change_pred + strategy_change_prey:.4f}")

print("\n" + "=" * 60)
print("🧬 Simulation completed successfully!")
print("The coevolutionary arms race has reached equilibrium.")

def success_probability(self, s_p, s_r):
    return s_p / (s_p + s_r + self.c)

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

# Set up the competition model parameters
class CompetitionModel:
    def __init__(self, R=100, c=0.5):
        """
        Initialize the competition model
        
        Parameters:
        R (float): Total resource availability
        c (float): Cost coefficient for competition investment
        """
        self.R = R  # Total resource pool
        self.c = c  # Cost of competition per unit squared
    
    def payoff_species1(self, x1, x2):
        """
        Calculate payoff for species 1
        
        π₁(x₁, x₂) = (x₁/(x₁ + x₂)) * R - c * x₁²
        """
        if x1 + x2 == 0:
            return -self.c * x1**2
        return (x1 / (x1 + x2)) * self.R - self.c * x1**2
    
    def payoff_species2(self, x1, x2):
        """
        Calculate payoff for species 2
        
        π₂(x₁, x₂) = (x₂/(x₁ + x₂)) * R - c * x₂²
        """
        if x1 + x2 == 0:
            return -self.c * x2**2
        return (x2 / (x1 + x2)) * self.R - self.c * x2**2
    
    def find_nash_equilibrium(self):
        """
        Find Nash equilibrium using optimization
        At Nash equilibrium, each species maximizes its payoff given the other's strategy
        """
        def objective(vars):
            x1, x2 = vars
            # We minimize the negative sum of squared deviations from best response
            # Best response for species 1: ∂π₁/∂x₁ = 0
            # Best response for species 2: ∂π₂/∂x₂ = 0
            
            # Partial derivatives
            if x1 + x2 > 0:
                dpi1_dx1 = (x2 * self.R) / ((x1 + x2)**2) - 2 * self.c * x1
                dpi2_dx2 = (x1 * self.R) / ((x1 + x2)**2) - 2 * self.c * x2
            else:
                dpi1_dx1 = -2 * self.c * x1
                dpi2_dx2 = -2 * self.c * x2
            
            return dpi1_dx1**2 + dpi2_dx2**2
        
        # Initial guess
        x0 = [5, 5]
        
        # Constraints: x1, x2 >= 0
        bounds = [(0, None), (0, None)]
        
        result = minimize(objective, x0, bounds=bounds, method='L-BFGS-B')
        return result.x
    
    def analytical_nash_equilibrium(self):
        """
        Calculate Nash equilibrium analytically
        
        From first-order conditions:
        ∂π₁/∂x₁ = x₂R/(x₁+x₂)² - 2cx₁ = 0
        ∂π₂/∂x₂ = x₁R/(x₁+x₂)² - 2cx₂ = 0
        
        By symmetry in this model, x₁* = x₂* at equilibrium
        Solving: R/(4x*) - 2cx* = 0
        Therefore: x* = √(R/(8c))
        """
        x_star = np.sqrt(self.R / (8 * self.c))
        return x_star, x_star

# Initialize the model
model = CompetitionModel(R=100, c=0.5)

# Find Nash equilibrium
nash_numerical = model.find_nash_equilibrium()
nash_analytical = model.analytical_nash_equilibrium()

print("Competition Strategy Optimization Results")
print("=" * 50)
print(f"Model Parameters:")
print(f"  Resource Pool (R): {model.R}")
print(f"  Competition Cost (c): {model.c}")
print()
print(f"Nash Equilibrium (Numerical): x₁* = {nash_numerical[0]:.3f}, x₂* = {nash_numerical[1]:.3f}")
print(f"Nash Equilibrium (Analytical): x₁* = {nash_analytical[0]:.3f}, x₂* = {nash_analytical[1]:.3f}")
print()

# Calculate payoffs at equilibrium
payoff1_nash = model.payoff_species1(nash_analytical[0], nash_analytical[1])
payoff2_nash = model.payoff_species2(nash_analytical[0], nash_analytical[1])
print(f"Equilibrium Payoffs:")
print(f"  Species 1: π₁* = {payoff1_nash:.3f}")
print(f"  Species 2: π₂* = {payoff2_nash:.3f}")
print(f"  Total Payoff: {payoff1_nash + payoff2_nash:.3f}")

# Create strategy space for visualization
x1_range = np.linspace(0, 20, 100)
x2_range = np.linspace(0, 20, 100)
X1, X2 = np.meshgrid(x1_range, x2_range)

# Calculate payoff surfaces
Payoff1 = np.zeros_like(X1)
Payoff2 = np.zeros_like(X2)

for i in range(len(x1_range)):
    for j in range(len(x2_range)):
        Payoff1[j, i] = model.payoff_species1(X1[j, i], X2[j, i])
        Payoff2[j, i] = model.payoff_species2(X1[j, i], X2[j, i])

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

# 1. 3D Payoff Surface for Species 1
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
surf1 = ax1.plot_surface(X1, X2, Payoff1, cmap='viridis', alpha=0.8)
ax1.scatter([nash_analytical[0]], [nash_analytical[1]], [payoff1_nash], 
           color='red', s=100, label='Nash Equilibrium')
ax1.set_xlabel('Species 1 Investment (x₁)')
ax1.set_ylabel('Species 2 Investment (x₂)')
ax1.set_zlabel('Species 1 Payoff (π₁)')
ax1.set_title('Species 1 Payoff Surface')
plt.colorbar(surf1, ax=ax1, shrink=0.5)

# 2. 3D Payoff Surface for Species 2
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
surf2 = ax2.plot_surface(X1, X2, Payoff2, cmap='plasma', alpha=0.8)
ax2.scatter([nash_analytical[0]], [nash_analytical[1]], [payoff2_nash], 
           color='red', s=100, label='Nash Equilibrium')
ax2.set_xlabel('Species 1 Investment (x₁)')
ax2.set_ylabel('Species 2 Investment (x₂)')
ax2.set_zlabel('Species 2 Payoff (π₂)')
ax2.set_title('Species 2 Payoff Surface')
plt.colorbar(surf2, ax=ax2, shrink=0.5)

# 3. Contour plot with Nash equilibrium
ax3 = fig.add_subplot(2, 3, 3)
contour1 = ax3.contour(X1, X2, Payoff1, levels=20, colors='blue', alpha=0.5)
contour2 = ax3.contour(X1, X2, Payoff2, levels=20, colors='red', alpha=0.5)
ax3.clabel(contour1, inline=True, fontsize=8, fmt='%.1f')
ax3.clabel(contour2, inline=True, fontsize=8, fmt='%.1f')
ax3.plot(nash_analytical[0], nash_analytical[1], 'ko', markersize=10, label='Nash Equilibrium')
ax3.set_xlabel('Species 1 Investment (x₁)')
ax3.set_ylabel('Species 2 Investment (x₂)')
ax3.set_title('Payoff Contours (Blue: Species 1, Red: Species 2)')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Best Response Functions
def best_response_1(x2):
    """Best response function for species 1 given species 2's strategy"""
    if x2 == 0:
        return 0
    # From ∂π₁/∂x₁ = 0: x₂R/(x₁+x₂)² = 2cx₁
    # This gives us a cubic equation that we solve numerically
    def objective(x1):
        if x1 + x2 <= 0:
            return float('inf')
        return abs((x2 * model.R) / ((x1 + x2)**2) - 2 * model.c * x1)
    
    result = minimize(objective, [1], bounds=[(0, None)], method='L-BFGS-B')
    return result.x[0] if result.success else 0

def best_response_2(x1):
    """Best response function for species 2 given species 1's strategy"""
    if x1 == 0:
        return 0
    def objective(x2):
        if x1 + x2 <= 0:
            return float('inf')
        return abs((x1 * model.R) / ((x1 + x2)**2) - 2 * model.c * x2)
    
    result = minimize(objective, [1], bounds=[(0, None)], method='L-BFGS-B')
    return result.x[0] if result.success else 0

ax4 = fig.add_subplot(2, 3, 4)
x_br = np.linspace(0.1, 15, 50)
br1 = [best_response_1(x) for x in x_br]
br2 = [best_response_2(x) for x in x_br]

ax4.plot(x_br, br1, 'b-', linewidth=2, label='BR₁(x₂)')
ax4.plot(br2, x_br, 'r-', linewidth=2, label='BR₂(x₁)')
ax4.plot(nash_analytical[0], nash_analytical[1], 'ko', markersize=10, label='Nash Equilibrium')
ax4.plot([0, 15], [0, 15], 'k--', alpha=0.5, label='45° line')
ax4.set_xlabel('Investment Level')
ax4.set_ylabel('Best Response')
ax4.set_title('Best Response Functions')
ax4.legend()
ax4.grid(True, alpha=0.3)

# 5. Payoff comparison across different strategies
strategies = np.linspace(1, 15, 15)
payoffs_coop = []  # Both species use same strategy
payoffs_defect = []  # One species increases investment

for s in strategies:
    # Cooperative strategy (both use s)
    p1_coop = model.payoff_species1(s, s)
    p2_coop = model.payoff_species2(s, s)
    payoffs_coop.append(p1_coop + p2_coop)
    
    # Defection strategy (one uses s, other uses Nash)
    p1_defect = model.payoff_species1(s, nash_analytical[1])
    p2_defect = model.payoff_species2(nash_analytical[0], s)
    payoffs_defect.append(p1_defect + p2_defect)

ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(strategies, payoffs_coop, 'g-', linewidth=2, label='Cooperative (same strategy)')
ax5.plot(strategies, payoffs_defect, 'r-', linewidth=2, label='Mixed strategies')
ax5.axhline(y=payoff1_nash + payoff2_nash, color='k', linestyle='--', 
           label=f'Nash Equilibrium Total = {payoff1_nash + payoff2_nash:.2f}')
ax5.axvline(x=nash_analytical[0], color='k', linestyle=':', alpha=0.5)
ax5.set_xlabel('Investment Level')
ax5.set_ylabel('Total Payoff')
ax5.set_title('Strategy Comparison')
ax5.legend()
ax5.grid(True, alpha=0.3)

# 6. Sensitivity Analysis
costs = np.linspace(0.1, 1.0, 20)
nash_investments = []
total_payoffs = []

for c in costs:
    temp_model = CompetitionModel(R=100, c=c)
    x_nash = np.sqrt(temp_model.R / (8 * c))
    nash_investments.append(x_nash)
    payoff = temp_model.payoff_species1(x_nash, x_nash) + temp_model.payoff_species2(x_nash, x_nash)
    total_payoffs.append(payoff)

ax6 = fig.add_subplot(2, 3, 6)
ax6_twin = ax6.twinx()

line1 = ax6.plot(costs, nash_investments, 'b-', linewidth=2, label='Nash Investment')
line2 = ax6_twin.plot(costs, total_payoffs, 'r-', linewidth=2, label='Total Payoff')

ax6.set_xlabel('Competition Cost (c)')
ax6.set_ylabel('Nash Investment Level', color='b')
ax6_twin.set_ylabel('Total Payoff', color='r')
ax6.set_title('Sensitivity to Competition Cost')

# Combine legends
lines1, labels1 = ax6.get_legend_handles_labels()
lines2, labels2 = ax6_twin.get_legend_handles_labels()
ax6.legend(lines1 + lines2, labels1 + labels2, loc='center right')

ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print additional analysis
print("\nDetailed Analysis:")
print("=" * 50)
print(f"At Nash equilibrium, each species invests {nash_analytical[0]:.3f} units")
print(f"This results in each species getting {model.R/2:.1f} units of resource")
print(f"Total competition investment: {2*nash_analytical[0]:.3f} units")
print(f"Total cost incurred: {2 * model.c * nash_analytical[0]**2:.3f} units")
print(f"Efficiency loss compared to no competition: {2 * model.c * nash_analytical[0]**2:.3f} units")