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

# Set up plotting style
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

class QuantumEntanglementPurification:
    def __init__(self):
        """Initialize quantum entanglement purification optimizer"""
        # Define Pauli matrices
        self.I = np.array([[1, 0], [0, 1]], dtype=complex)
        self.X = np.array([[0, 1], [1, 0]], dtype=complex)
        self.Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
        self.Z = np.array([[1, 0], [0, -1]], dtype=complex)
        
        # Define Bell states
        self.phi_plus = np.array([1, 0, 0, 1], dtype=complex) / np.sqrt(2)
        self.phi_minus = np.array([1, 0, 0, -1], dtype=complex) / np.sqrt(2)
        self.psi_plus = np.array([0, 1, 1, 0], dtype=complex) / np.sqrt(2)
        self.psi_minus = np.array([0, 1, -1, 0], dtype=complex) / np.sqrt(2)
        
    def werner_state(self, p):
        """Create a Werner state with fidelity parameter p"""
        phi_plus_dm = np.outer(self.phi_plus, np.conj(self.phi_plus))
        identity_4 = np.eye(4, dtype=complex)
        return p * phi_plus_dm + (1 - p) * identity_4 / 4
    
    def rotation_operator(self, theta, phi=0):
        """Create rotation operator R(theta, phi) = exp(-i*theta*(cos(phi)*X + sin(phi)*Y)/2)"""
        return np.cos(theta/2) * self.I - 1j * np.sin(theta/2) * (np.cos(phi) * self.X + np.sin(phi) * self.Y)
    
    def bilateral_rotation(self, theta_a, theta_b, phi_a=0, phi_b=0):
        """Apply bilateral rotation to both qubits"""
        R_a = self.rotation_operator(theta_a, phi_a)
        R_b = self.rotation_operator(theta_b, phi_b)
        return np.kron(R_a, R_b)
    
    def cnot_gate(self):
        """CNOT gate matrix"""
        return np.array([[1, 0, 0, 0],
                        [0, 1, 0, 0],
                        [0, 0, 0, 1],
                        [0, 0, 1, 0]], dtype=complex)
    
    def bilateral_cnot(self):
        """Apply bilateral CNOT operations"""
        cnot = self.cnot_gate()
        return np.kron(cnot, cnot)
    
    def measurement_projector(self, outcome1, outcome2):
        """
        Create measurement projector for 4-qubit system
        Projects onto |outcome1⟩⊗|outcome2⟩ for qubits 2 and 4 (ancilla qubits)
        """
        # Create single-qubit projectors
        proj_0 = np.array([[1, 0], [0, 0]], dtype=complex)
        proj_1 = np.array([[0, 0], [0, 1]], dtype=complex)
        
        projectors = [proj_0, proj_1]
        
        # For 4-qubit system: qubit1 ⊗ qubit2 ⊗ qubit3 ⊗ qubit4
        # We measure qubits 2 and 4 (ancilla qubits)
        identity = np.eye(2, dtype=complex)
        
        # Create 4-qubit measurement projector
        # I ⊗ P_outcome1 ⊗ I ⊗ P_outcome2
        projector = np.kron(np.kron(np.kron(identity, projectors[outcome1]), 
                                   identity), projectors[outcome2])
        
        return projector
    
    def bbpssw_purification(self, p_initial, theta_a, theta_b, phi_a=0, phi_b=0):
        """
        Perform BBPSSW purification protocol
        
        Args:
            p_initial: Initial fidelity parameter
            theta_a, theta_b: Rotation angles for Alice and Bob
            phi_a, phi_b: Rotation phases for Alice and Bob
        
        Returns:
            success_probability: Probability of successful purification
            final_fidelity: Fidelity after purification
        """
        # Create two copies of Werner state
        rho_1 = self.werner_state(p_initial)
        rho_2 = self.werner_state(p_initial)
        
        # Initial 4-qubit state (tensor product of two Werner states)
        rho_initial = np.kron(rho_1, rho_2)
        
        # Apply bilateral rotations
        U_rot = np.kron(self.bilateral_rotation(theta_a, theta_b, phi_a, phi_b),
                       self.bilateral_rotation(theta_a, theta_b, phi_a, phi_b))
        rho_after_rot = U_rot @ rho_initial @ np.conj(U_rot).T
        
        # Apply bilateral CNOT
        U_cnot = self.bilateral_cnot()
        rho_after_cnot = U_cnot @ rho_after_rot @ np.conj(U_cnot).T
        
        # Calculate success probability and final state
        success_prob = 0
        final_state = np.zeros((16, 16), dtype=complex)
        
        # Sum over matching measurement outcomes (both ancilla qubits give same result)
        for outcome in [0, 1]:
            # Measurement projector for matching outcomes
            P_meas = self.measurement_projector(outcome, outcome)
            
            # Post-measurement state
            rho_post = P_meas @ rho_after_cnot @ P_meas
            prob = np.real(np.trace(rho_post))
            
            if prob > 1e-10:  # Avoid numerical issues
                success_prob += prob
                final_state += rho_post
        
        # Normalize final state
        if success_prob > 1e-10:
            final_state = final_state / success_prob
            
            # Trace out ancilla qubits to get the 2-qubit purified state
            # The purified state is in qubits 1 and 3
            purified_state = self.partial_trace_ancilla(final_state)
            
            # Calculate fidelity with maximally entangled state
            phi_plus_dm = np.outer(self.phi_plus, np.conj(self.phi_plus))
            final_fidelity = np.real(np.trace(purified_state @ phi_plus_dm))
        else:
            final_fidelity = 0
        
        return success_prob, final_fidelity
    
    def partial_trace_ancilla(self, rho_4qubit):
        """
        Trace out ancilla qubits (qubits 2 and 4) from 4-qubit state
        Returns 2-qubit state of qubits 1 and 3
        """
        # Reshape to group qubits properly
        rho_reshaped = rho_4qubit.reshape(2, 2, 2, 2, 2, 2, 2, 2)
        
        # Trace out qubits 2 and 4 (indices 1 and 3 in reshaped array)
        rho_2qubit = np.trace(rho_reshaped, axis1=1, axis2=5)  # Trace out qubit 2
        rho_2qubit = np.trace(rho_2qubit, axis1=1, axis2=3)    # Trace out qubit 4
        
        # Reshape back to 4x4 matrix
        return rho_2qubit.reshape(4, 4)
    
    def optimize_purification(self, p_initial, method='bounded'):
        """
        Optimize BBPSSW purification protocol parameters
        
        Args:
            p_initial: Initial fidelity parameter
            method: Optimization method ('bounded' or 'symmetric')
        
        Returns:
            optimal_params: Optimized parameters
            optimal_results: Results with optimal parameters
        """
        def objective(params):
            """Objective function to maximize: weighted sum of success probability and fidelity improvement"""
            if method == 'symmetric':
                theta_a = theta_b = params[0]
                phi_a = phi_b = 0
            else:
                theta_a, theta_b = params
                phi_a = phi_b = 0
            
            success_prob, final_fidelity = self.bbpssw_purification(p_initial, theta_a, theta_b, phi_a, phi_b)
            
            # Objective: maximize success probability × fidelity improvement
            fidelity_improvement = max(0, final_fidelity - p_initial)  # Ensure non-negative
            return -(success_prob * fidelity_improvement)  # Negative because we minimize
        
        # Set bounds and initial guess
        if method == 'symmetric':
            bounds = [(0, 2*np.pi)]
            x0 = [np.pi/2]
        else:
            bounds = [(0, 2*np.pi), (0, 2*np.pi)]
            x0 = [np.pi/2, np.pi/2]
        
        # Optimize
        result = minimize(objective, x0, method='L-BFGS-B', bounds=bounds)
        
        # Extract optimal parameters
        if method == 'symmetric':
            optimal_theta_a = optimal_theta_b = result.x[0]
        else:
            optimal_theta_a, optimal_theta_b = result.x
        
        # Calculate results with optimal parameters
        success_prob, final_fidelity = self.bbpssw_purification(
            p_initial, optimal_theta_a, optimal_theta_b
        )
        
        optimal_params = {
            'theta_a': optimal_theta_a,
            'theta_b': optimal_theta_b,
            'phi_a': 0,
            'phi_b': 0
        }
        
        optimal_results = {
            'success_probability': success_prob,
            'final_fidelity': final_fidelity,
            'fidelity_improvement': final_fidelity - p_initial,
            'optimization_success': result.success
        }
        
        return optimal_params, optimal_results

# Initialize the purification optimizer
purifier = QuantumEntanglementPurification()

# Example 1: Optimize for a specific initial fidelity
p_initial_example = 0.6
print(f"Optimizing purification for initial fidelity p = {p_initial_example}")
print("=" * 60)

# Optimize with symmetric constraint
optimal_params_sym, optimal_results_sym = purifier.optimize_purification(p_initial_example, method='symmetric')

print("Symmetric optimization results:")
print(f"Optimal θ_A = θ_B = {optimal_params_sym['theta_a']:.4f} rad ({optimal_params_sym['theta_a']*180/np.pi:.2f}°)")
print(f"Success probability: {optimal_results_sym['success_probability']:.4f}")
print(f"Final fidelity: {optimal_results_sym['final_fidelity']:.4f}")
print(f"Fidelity improvement: {optimal_results_sym['fidelity_improvement']:.4f}")
print()

# Optimize without symmetric constraint
optimal_params_full, optimal_results_full = purifier.optimize_purification(p_initial_example, method='bounded')

print("Full optimization results:")
print(f"Optimal θ_A = {optimal_params_full['theta_a']:.4f} rad ({optimal_params_full['theta_a']*180/np.pi:.2f}°)")
print(f"Optimal θ_B = {optimal_params_full['theta_b']:.4f} rad ({optimal_params_full['theta_b']*180/np.pi:.2f}°)")
print(f"Success probability: {optimal_results_full['success_probability']:.4f}")
print(f"Final fidelity: {optimal_results_full['final_fidelity']:.4f}")
print(f"Fidelity improvement: {optimal_results_full['fidelity_improvement']:.4f}")
print()

# Example 2: Analyze purification performance across different initial fidelities
print("Analyzing purification performance across different initial fidelities...")
print("=" * 60)

p_values = np.linspace(0.3, 0.9, 15)  # Reduced range for better convergence
results_analysis = {
    'p_initial': [],
    'success_prob': [],
    'final_fidelity': [],
    'fidelity_improvement': [],
    'optimal_theta_a': [],
    'optimal_theta_b': []
}

for p in p_values:
    try:
        optimal_params, optimal_results = purifier.optimize_purification(p, method='symmetric')
        
        results_analysis['p_initial'].append(p)
        results_analysis['success_prob'].append(optimal_results['success_probability'])
        results_analysis['final_fidelity'].append(optimal_results['final_fidelity'])
        results_analysis['fidelity_improvement'].append(optimal_results['fidelity_improvement'])
        results_analysis['optimal_theta_a'].append(optimal_params['theta_a'])
        results_analysis['optimal_theta_b'].append(optimal_params['theta_b'])
        
        print(f"p = {p:.3f}: Success = {optimal_results['success_probability']:.4f}, Final F = {optimal_results['final_fidelity']:.4f}")
        
    except Exception as e:
        print(f"Failed for p = {p:.3f}: {e}")

# Create comprehensive plots
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Quantum Entanglement Purification Optimization Analysis', fontsize=16, fontweight='bold')

# Plot 1: Success Probability vs Initial Fidelity
axes[0, 0].plot(results_analysis['p_initial'], results_analysis['success_prob'], 'o-', linewidth=2, markersize=6)
axes[0, 0].set_xlabel('Initial Fidelity (p)')
axes[0, 0].set_ylabel('Success Probability')
axes[0, 0].set_title('Success Probability vs Initial Fidelity')
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Final Fidelity vs Initial Fidelity
axes[0, 1].plot(results_analysis['p_initial'], results_analysis['final_fidelity'], 'o-', linewidth=2, markersize=6, color='orange')
axes[0, 1].plot(results_analysis['p_initial'], results_analysis['p_initial'], '--', color='gray', alpha=0.7, label='No improvement line')
axes[0, 1].set_xlabel('Initial Fidelity (p)')
axes[0, 1].set_ylabel('Final Fidelity')
axes[0, 1].set_title('Final Fidelity vs Initial Fidelity')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Fidelity Improvement vs Initial Fidelity
axes[0, 2].plot(results_analysis['p_initial'], results_analysis['fidelity_improvement'], 'o-', linewidth=2, markersize=6, color='green')
axes[0, 2].axhline(y=0, color='gray', linestyle='--', alpha=0.7)
axes[0, 2].set_xlabel('Initial Fidelity (p)')
axes[0, 2].set_ylabel('Fidelity Improvement')
axes[0, 2].set_title('Fidelity Improvement vs Initial Fidelity')
axes[0, 2].grid(True, alpha=0.3)

# Plot 4: Optimal Rotation Angles
axes[1, 0].plot(results_analysis['p_initial'], np.array(results_analysis['optimal_theta_a'])*180/np.pi, 'o-', linewidth=2, markersize=6, label='θ_A = θ_B', color='purple')
axes[1, 0].set_xlabel('Initial Fidelity (p)')
axes[1, 0].set_ylabel('Optimal Rotation Angle (degrees)')
axes[1, 0].set_title('Optimal Rotation Angles vs Initial Fidelity')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Plot 5: 2D Performance Heatmap for a specific initial fidelity
theta_range = np.linspace(0, np.pi, 20)  # Reduced resolution for faster computation
performance_map = np.zeros((len(theta_range), len(theta_range)))

p_heatmap = 0.6  # Use the example initial fidelity
for i, theta_a in enumerate(theta_range):
    for j, theta_b in enumerate(theta_range):
        try:
            success_prob, final_fidelity = purifier.bbpssw_purification(p_heatmap, theta_a, theta_b)
            performance_map[i, j] = success_prob * max(0, final_fidelity - p_heatmap)
        except:
            performance_map[i, j] = 0

im = axes[1, 1].imshow(performance_map, extent=[0, 180, 0, 180], origin='lower', cmap='viridis', aspect='auto')
axes[1, 1].set_xlabel('θ_B (degrees)')
axes[1, 1].set_ylabel('θ_A (degrees)')
axes[1, 1].set_title(f'Performance Heatmap (p = {p_heatmap})')
plt.colorbar(im, ax=axes[1, 1], label='Success Prob × Fidelity Improvement')

# Plot 6: Efficiency Analysis
if len(results_analysis['success_prob']) > 0:
    efficiency = np.array(results_analysis['success_prob']) * np.array(results_analysis['fidelity_improvement'])
    axes[1, 2].plot(results_analysis['p_initial'], efficiency, 'o-', linewidth=2, markersize=6, color='brown')
    axes[1, 2].set_xlabel('Initial Fidelity (p)')
    axes[1, 2].set_ylabel('Purification Efficiency')
    axes[1, 2].set_title('Purification Efficiency vs Initial Fidelity')
    axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print summary statistics
if len(results_analysis['fidelity_improvement']) > 0:
    print("\nSummary Statistics:")
    print("=" * 40)
    max_improvement_idx = np.argmax(results_analysis['fidelity_improvement'])
    efficiency = np.array(results_analysis['success_prob']) * np.array(results_analysis['fidelity_improvement'])
    max_efficiency_idx = np.argmax(efficiency)

    print(f"Maximum fidelity improvement: {results_analysis['fidelity_improvement'][max_improvement_idx]:.4f}")
    print(f"  at initial fidelity p = {results_analysis['p_initial'][max_improvement_idx]:.3f}")
    print(f"Maximum efficiency: {efficiency[max_efficiency_idx]:.4f}")
    print(f"  at initial fidelity p = {results_analysis['p_initial'][max_efficiency_idx]:.3f}")
    print(f"Average success probability: {np.mean(results_analysis['success_prob']):.4f}")
    print(f"Average fidelity improvement: {np.mean(results_analysis['fidelity_improvement']):.4f}")

# Detailed theoretical analysis
print("\nTheoretical Analysis:")
print("=" * 40)
print("The BBPSSW protocol works by:")
print("1. Creating quantum correlations between two noisy Bell pairs")
print("2. Using bilateral rotations to optimize interference patterns")
print("3. Performing bilateral CNOT to create entanglement between pairs")
print("4. Post-selecting on matching measurement outcomes")
print("5. The surviving state has enhanced fidelity due to quantum interference")
print()
print("Key insights from optimization:")
print("- Optimal angles depend on initial fidelity level")
print("- Maximum improvement occurs at intermediate fidelities (~0.5-0.7)")
print("- Success probability trades off with fidelity improvement")
print("- Protocol efficiency peaks at specific operating points")

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

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

class HabitatOptimizer:
    def __init__(self, alpha=10, beta=0.5, gamma=5, delta=0.1):
        """
        Initialize the habitat selection model
        
        Parameters:
        alpha: food acquisition efficiency
        beta: territory defense cost per unit area
        gamma: competition pressure coefficient
        delta: competition decay rate
        """
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma
        self.delta = delta
        
    def fitness_function(self, params):
        """
        Calculate fitness based on territory area and distance from competitors
        
        Parameters:
        params: [area, distance] - territory area and distance from competitors
        
        Returns:
        fitness: negative fitness (for minimization)
        """
        area, distance = params
        
        # Food benefit (diminishing returns with square root)
        food_benefit = self.alpha * np.sqrt(area)
        
        # Territory defense cost (linear with area)
        defense_cost = self.beta * area
        
        # Competition pressure (exponential decay with distance)
        competition_cost = self.gamma * np.exp(-self.delta * distance)
        
        fitness = food_benefit - defense_cost - competition_cost
        
        return -fitness  # Negative because we're minimizing
    
    def optimize_habitat(self, area_bounds=(1, 50), distance_bounds=(0, 20)):
        """
        Find optimal territory area and distance from competitors
        
        Parameters:
        area_bounds: (min, max) territory area
        distance_bounds: (min, max) distance from competitors
        
        Returns:
        result: optimization result
        """
        bounds = [area_bounds, distance_bounds]
        
        # Use differential evolution for global optimization
        result = differential_evolution(
            self.fitness_function,
            bounds,
            seed=42,
            maxiter=1000
        )
        
        return result
    
    def plot_fitness_landscape(self, area_range=(1, 50), distance_range=(0, 20)):
        """
        Create a 3D plot of the fitness landscape
        """
        # Create meshgrid for plotting
        area_vals = np.linspace(area_range[0], area_range[1], 50)
        distance_vals = np.linspace(distance_range[0], distance_range[1], 50)
        A, D = np.meshgrid(area_vals, distance_vals)
        
        # Calculate fitness for each point
        fitness_vals = np.zeros_like(A)
        for i in range(A.shape[0]):
            for j in range(A.shape[1]):
                fitness_vals[i, j] = -self.fitness_function([A[i, j], D[i, j]])
        
        # Create 3D plot
        fig = plt.figure(figsize=(12, 8))
        ax = fig.add_subplot(111, projection='3d')
        
        surface = ax.plot_surface(A, D, fitness_vals, cmap='viridis', alpha=0.7)
        ax.set_xlabel('Territory Area')
        ax.set_ylabel('Distance from Competitors')
        ax.set_zlabel('Fitness')
        ax.set_title('Fitness Landscape for Habitat Selection')
        
        plt.colorbar(surface, shrink=0.5)
        plt.tight_layout()
        plt.show()
        
        return A, D, fitness_vals
    
    def plot_fitness_components(self, optimal_area, optimal_distance):
        """
        Plot individual components of the fitness function
        """
        area_vals = np.linspace(1, 50, 100)
        distance_vals = np.linspace(0, 20, 100)
        
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 10))
        
        # Food benefit vs area
        food_benefit = self.alpha * np.sqrt(area_vals)
        ax1.plot(area_vals, food_benefit, 'g-', linewidth=2, label='Food Benefit')
        ax1.axvline(optimal_area, color='red', linestyle='--', alpha=0.7, label=f'Optimal Area = {optimal_area:.1f}')
        ax1.set_xlabel('Territory Area')
        ax1.set_ylabel('Food Benefit')
        ax1.set_title('Food Benefit vs Territory Area')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # Defense cost vs area
        defense_cost = self.beta * area_vals
        ax2.plot(area_vals, defense_cost, 'r-', linewidth=2, label='Defense Cost')
        ax2.axvline(optimal_area, color='red', linestyle='--', alpha=0.7, label=f'Optimal Area = {optimal_area:.1f}')
        ax2.set_xlabel('Territory Area')
        ax2.set_ylabel('Defense Cost')
        ax2.set_title('Defense Cost vs Territory Area')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        # Competition pressure vs distance
        competition_cost = self.gamma * np.exp(-self.delta * distance_vals)
        ax3.plot(distance_vals, competition_cost, 'orange', linewidth=2, label='Competition Pressure')
        ax3.axvline(optimal_distance, color='red', linestyle='--', alpha=0.7, label=f'Optimal Distance = {optimal_distance:.1f}')
        ax3.set_xlabel('Distance from Competitors')
        ax3.set_ylabel('Competition Pressure')
        ax3.set_title('Competition Pressure vs Distance')
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # Combined fitness function
        area_grid, distance_grid = np.meshgrid(area_vals, distance_vals)
        fitness_grid = np.zeros_like(area_grid)
        for i in range(area_grid.shape[0]):
            for j in range(area_grid.shape[1]):
                fitness_grid[i, j] = -self.fitness_function([area_grid[i, j], distance_grid[i, j]])
        
        contour = ax4.contour(area_grid, distance_grid, fitness_grid, levels=20, cmap='viridis')
        ax4.clabel(contour, inline=True, fontsize=8)
        ax4.plot(optimal_area, optimal_distance, 'r*', markersize=15, label=f'Optimal Point ({optimal_area:.1f}, {optimal_distance:.1f})')
        ax4.set_xlabel('Territory Area')
        ax4.set_ylabel('Distance from Competitors')
        ax4.set_title('Fitness Contour Plot')
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
    
    def sensitivity_analysis(self, optimal_area, optimal_distance):
        """
        Analyze sensitivity to parameter changes
        """
        # Parameter variations
        param_variations = {
            'alpha': np.linspace(5, 15, 20),
            'beta': np.linspace(0.1, 1.0, 20),
            'gamma': np.linspace(1, 10, 20),
            'delta': np.linspace(0.05, 0.2, 20)
        }
        
        fig, axes = plt.subplots(2, 2, figsize=(15, 10))
        axes = axes.flatten()
        
        for idx, (param_name, param_vals) in enumerate(param_variations.items()):
            optimal_areas = []
            optimal_distances = []
            optimal_fitness = []
            
            for param_val in param_vals:
                # Create temporary optimizer with modified parameter
                temp_params = {
                    'alpha': self.alpha,
                    'beta': self.beta,
                    'gamma': self.gamma,
                    'delta': self.delta
                }
                temp_params[param_name] = param_val
                
                temp_optimizer = HabitatOptimizer(**temp_params)
                result = temp_optimizer.optimize_habitat()
                
                optimal_areas.append(result.x[0])
                optimal_distances.append(result.x[1])
                optimal_fitness.append(-result.fun)
            
            # Plot results
            ax = axes[idx]
            ax2 = ax.twinx()
            
            line1 = ax.plot(param_vals, optimal_areas, 'b-', linewidth=2, label='Optimal Area')
            line2 = ax2.plot(param_vals, optimal_distances, 'r-', linewidth=2, label='Optimal Distance')
            
            ax.set_xlabel(f'{param_name}')
            ax.set_ylabel('Optimal Area', color='b')
            ax2.set_ylabel('Optimal Distance', color='r')
            ax.set_title(f'Sensitivity to {param_name}')
            ax.grid(True, alpha=0.3)
            
            # Add legend
            lines = line1 + line2
            labels = [l.get_label() for l in lines]
            ax.legend(lines, labels, loc='upper left')
        
        plt.tight_layout()
        plt.show()

# Initialize the habitat optimizer
optimizer = HabitatOptimizer(alpha=10, beta=0.5, gamma=5, delta=0.1)

# Find optimal habitat selection
print("Optimizing habitat selection...")
result = optimizer.optimize_habitat()

print(f"\nOptimization Results:")
print(f"Optimal territory area: {result.x[0]:.2f} units")
print(f"Optimal distance from competitors: {result.x[1]:.2f} units")
print(f"Maximum fitness: {-result.fun:.2f}")
print(f"Success: {result.success}")

# Calculate fitness components at optimal point
optimal_area, optimal_distance = result.x
food_benefit = optimizer.alpha * np.sqrt(optimal_area)
defense_cost = optimizer.beta * optimal_area
competition_cost = optimizer.gamma * np.exp(-optimizer.delta * optimal_distance)

print(f"\nFitness Components at Optimal Point:")
print(f"Food benefit: {food_benefit:.2f}")
print(f"Defense cost: {defense_cost:.2f}")
print(f"Competition cost: {competition_cost:.2f}")
print(f"Net fitness: {food_benefit - defense_cost - competition_cost:.2f}")

# Create visualizations
print("\nGenerating visualizations...")

# 1. 3D fitness landscape
print("Creating 3D fitness landscape...")
A, D, fitness_vals = optimizer.plot_fitness_landscape()

# 2. Fitness components analysis
print("Analyzing fitness components...")
optimizer.plot_fitness_components(optimal_area, optimal_distance)

# 3. Sensitivity analysis
print("Performing sensitivity analysis...")
optimizer.sensitivity_analysis(optimal_area, optimal_distance)

# Additional analysis: Compare different scenarios
print("\nComparing different environmental scenarios...")

scenarios = {
    'Low Competition': {'alpha': 10, 'beta': 0.5, 'gamma': 2, 'delta': 0.1},
    'High Competition': {'alpha': 10, 'beta': 0.5, 'gamma': 8, 'delta': 0.1},
    'Abundant Food': {'alpha': 15, 'beta': 0.5, 'gamma': 5, 'delta': 0.1},
    'Scarce Food': {'alpha': 5, 'beta': 0.5, 'gamma': 5, 'delta': 0.1},
    'High Defense Cost': {'alpha': 10, 'beta': 1.0, 'gamma': 5, 'delta': 0.1}
}

scenario_results = {}
for scenario_name, params in scenarios.items():
    scenario_optimizer = HabitatOptimizer(**params)
    scenario_result = scenario_optimizer.optimize_habitat()
    scenario_results[scenario_name] = {
        'area': scenario_result.x[0],
        'distance': scenario_result.x[1],
        'fitness': -scenario_result.fun
    }

# Plot scenario comparison
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))

scenarios_list = list(scenario_results.keys())
areas = [scenario_results[s]['area'] for s in scenarios_list]
distances = [scenario_results[s]['distance'] for s in scenarios_list]
fitness_values = [scenario_results[s]['fitness'] for s in scenarios_list]

ax1.bar(scenarios_list, areas, color='skyblue', alpha=0.7)
ax1.set_ylabel('Optimal Territory Area')
ax1.set_title('Optimal Territory Area by Scenario')
ax1.tick_params(axis='x', rotation=45)

ax2.bar(scenarios_list, distances, color='lightcoral', alpha=0.7)
ax2.set_ylabel('Optimal Distance from Competitors')
ax2.set_title('Optimal Distance by Scenario')
ax2.tick_params(axis='x', rotation=45)

ax3.bar(scenarios_list, fitness_values, color='lightgreen', alpha=0.7)
ax3.set_ylabel('Maximum Fitness')
ax3.set_title('Maximum Fitness by Scenario')
ax3.tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

print("\nScenario Analysis Results:")
for scenario, results in scenario_results.items():
    print(f"{scenario}:")
    print(f"  Optimal area: {results['area']:.2f}")
    print(f"  Optimal distance: {results['distance']:.2f}")
    print(f"  Maximum fitness: {results['fitness']:.2f}")
    print()

print("Analysis complete!")

def fitness_function(self, params):
    area, distance = params
    food_benefit = self.alpha * np.sqrt(area)
    defense_cost = self.beta * area
    competition_cost = self.gamma * np.exp(-self.delta * distance)
    fitness = food_benefit - defense_cost - competition_cost
    return -fitness  # Negative for minimization

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

# Set up the plotting style
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

class BirdWingOptimizer:
    def __init__(self):
        # Physical constants and parameters
        self.air_density = 1.225  # kg/m³ at sea level
        self.bird_mass = 0.5  # kg (typical small bird)
        self.flight_speed = 10.0  # m/s
        self.gravity = 9.81  # m/s²
        
        # Wing design constraints
        self.min_aspect_ratio = 3.0
        self.max_aspect_ratio = 12.0
        self.min_wingspan = 0.3  # meters
        self.max_wingspan = 1.5  # meters
        
    def calculate_wing_properties(self, wingspan, aspect_ratio, chord_coeffs):
        """
        Calculate wing properties based on design parameters
        
        Parameters:
        - wingspan: total wingspan in meters
        - aspect_ratio: wingspan²/wing_area
        - chord_coeffs: coefficients for chord distribution polynomial
        """
        # Calculate wing area from aspect ratio
        wing_area = wingspan**2 / aspect_ratio
        
        # Mean chord length
        mean_chord = wing_area / wingspan
        
        # Chord distribution along span (using polynomial)
        span_positions = np.linspace(-wingspan/2, wingspan/2, 50)
        normalized_positions = span_positions / (wingspan/2)
        
        # Chord distribution: c(y) = c_root * (1 + a₁*η + a₂*η²)
        # where η is normalized span position
        chord_distribution = mean_chord * (1 + 
                                         chord_coeffs[0] * normalized_positions + 
                                         chord_coeffs[1] * normalized_positions**2)
        
        # Ensure non-negative chord lengths
        chord_distribution = np.maximum(chord_distribution, 0.01)
        
        return wing_area, mean_chord, span_positions, chord_distribution
    
    def calculate_oswald_efficiency(self, aspect_ratio, chord_coeffs):
        """
        Calculate Oswald efficiency factor based on wing geometry
        """
        # Simplified model: efficiency depends on aspect ratio and taper
        base_efficiency = 0.8
        ar_factor = aspect_ratio / (aspect_ratio + 2)
        taper_penalty = 0.1 * abs(chord_coeffs[0])  # Penalty for extreme taper
        
        return base_efficiency * ar_factor - taper_penalty
    
    def calculate_lift_coefficient(self, wing_area):
        """
        Calculate required lift coefficient for steady flight
        """
        # Lift must equal weight: L = mg = 0.5 * ρ * V² * S * C_L
        required_lift = self.bird_mass * self.gravity
        dynamic_pressure = 0.5 * self.air_density * self.flight_speed**2
        
        return required_lift / (dynamic_pressure * wing_area)
    
    def calculate_drag_coefficient(self, lift_coeff, aspect_ratio, oswald_eff):
        """
        Calculate total drag coefficient
        """
        # Profile drag (simplified)
        profile_drag = 0.008  # Typical value for bird wings
        
        # Induced drag
        induced_drag = lift_coeff**2 / (np.pi * aspect_ratio * oswald_eff)
        
        return profile_drag + induced_drag
    
    def objective_function(self, params):
        """
        Objective function to minimize (negative lift-to-drag ratio)
        """
        wingspan, aspect_ratio, chord_coeff1, chord_coeff2 = params
        chord_coeffs = [chord_coeff1, chord_coeff2]
        
        # Calculate wing properties
        wing_area, mean_chord, _, _ = self.calculate_wing_properties(
            wingspan, aspect_ratio, chord_coeffs)
        
        # Calculate aerodynamic coefficients
        lift_coeff = self.calculate_lift_coefficient(wing_area)
        oswald_eff = self.calculate_oswald_efficiency(aspect_ratio, chord_coeffs)
        drag_coeff = self.calculate_drag_coefficient(lift_coeff, aspect_ratio, oswald_eff)
        
        # Lift-to-drag ratio (we want to maximize this, so minimize its negative)
        ld_ratio = lift_coeff / drag_coeff
        
        # Add penalties for extreme designs
        penalty = 0
        if aspect_ratio > 10:
            penalty += (aspect_ratio - 10)**2
        if abs(chord_coeff1) > 0.5:
            penalty += (abs(chord_coeff1) - 0.5)**2
            
        return -ld_ratio + penalty
    
    def optimize_wing(self):
        """
        Perform wing optimization
        """
        # Initial guess: [wingspan, aspect_ratio, chord_coeff1, chord_coeff2]
        initial_guess = [0.8, 7.0, -0.2, 0.1]
        
        # Bounds for parameters
        bounds = [
            (self.min_wingspan, self.max_wingspan),  # wingspan
            (self.min_aspect_ratio, self.max_aspect_ratio),  # aspect ratio
            (-0.5, 0.5),  # chord coefficient 1
            (-0.3, 0.3)   # chord coefficient 2
        ]
        
        # Optimize
        result = minimize(self.objective_function, initial_guess, 
                         bounds=bounds, method='L-BFGS-B')
        
        return result

# Create optimizer and run optimization
optimizer = BirdWingOptimizer()
optimization_result = optimizer.optimize_wing()

# Extract optimal parameters
optimal_wingspan, optimal_aspect_ratio, optimal_chord1, optimal_chord2 = optimization_result.x
optimal_chord_coeffs = [optimal_chord1, optimal_chord2]

print("=== Wing Optimization Results ===")
print(f"Optimal Wingspan: {optimal_wingspan:.3f} m")
print(f"Optimal Aspect Ratio: {optimal_aspect_ratio:.3f}")
print(f"Chord Coefficients: [{optimal_chord1:.3f}, {optimal_chord2:.3f}]")
print(f"Optimization Success: {optimization_result.success}")
print(f"Final Objective Value: {optimization_result.fun:.6f}")

# Calculate properties of optimal wing
wing_area, mean_chord, span_positions, chord_distribution = optimizer.calculate_wing_properties(
    optimal_wingspan, optimal_aspect_ratio, optimal_chord_coeffs)

lift_coeff = optimizer.calculate_lift_coefficient(wing_area)
oswald_eff = optimizer.calculate_oswald_efficiency(optimal_aspect_ratio, optimal_chord_coeffs)
drag_coeff = optimizer.calculate_drag_coefficient(lift_coeff, optimal_aspect_ratio, oswald_eff)
ld_ratio = lift_coeff / drag_coeff

print(f"\n=== Aerodynamic Properties ===")
print(f"Wing Area: {wing_area:.4f} m²")
print(f"Mean Chord: {mean_chord:.4f} m")
print(f"Lift Coefficient: {lift_coeff:.4f}")
print(f"Drag Coefficient: {drag_coeff:.4f}")
print(f"Oswald Efficiency: {oswald_eff:.4f}")
print(f"Lift-to-Drag Ratio: {ld_ratio:.2f}")

# Compare with different bird species (reference data)
bird_species = {
    'Sparrow': {'wingspan': 0.22, 'aspect_ratio': 4.5},
    'Crow': {'wingspan': 0.85, 'aspect_ratio': 6.2},
    'Seagull': {'wingspan': 1.2, 'aspect_ratio': 8.1},
    'Albatross': {'wingspan': 3.1, 'aspect_ratio': 15.2}
}

print(f"\n=== Comparison with Real Birds ===")
for species, data in bird_species.items():
    print(f"{species}: Wingspan={data['wingspan']:.2f}m, AR={data['aspect_ratio']:.1f}")

# Visualization
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# 1. Wing planform view
ax1.fill_between(span_positions, chord_distribution/2, -chord_distribution/2, 
                 alpha=0.7, color='skyblue', label='Optimal Wing')
ax1.plot(span_positions, chord_distribution/2, 'b-', linewidth=2)
ax1.plot(span_positions, -chord_distribution/2, 'b-', linewidth=2)
ax1.set_xlabel('Span Position (m)')
ax1.set_ylabel('Chord Length (m)')
ax1.set_title('Optimal Wing Planform')
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.set_aspect('equal')

# 2. Chord distribution
ax2.plot(span_positions, chord_distribution, 'r-', linewidth=2, label='Chord Distribution')
ax2.axhline(y=mean_chord, color='g', linestyle='--', alpha=0.7, label=f'Mean Chord ({mean_chord:.3f}m)')
ax2.set_xlabel('Span Position (m)')
ax2.set_ylabel('Chord Length (m)')
ax2.set_title('Chord Distribution Along Wingspan')
ax2.grid(True, alpha=0.3)
ax2.legend()

# 3. Aspect ratio vs L/D ratio analysis
aspect_ratios = np.linspace(3, 12, 50)
ld_ratios = []

for ar in aspect_ratios:
    # Keep other parameters constant
    test_area = optimal_wingspan**2 / ar
    test_lift = optimizer.calculate_lift_coefficient(test_area)
    test_oswald = optimizer.calculate_oswald_efficiency(ar, optimal_chord_coeffs)
    test_drag = optimizer.calculate_drag_coefficient(test_lift, ar, test_oswald)
    ld_ratios.append(test_lift / test_drag)

ax3.plot(aspect_ratios, ld_ratios, 'g-', linewidth=2, label='L/D vs Aspect Ratio')
ax3.axvline(x=optimal_aspect_ratio, color='r', linestyle='--', 
            label=f'Optimal AR ({optimal_aspect_ratio:.2f})')
ax3.axhline(y=ld_ratio, color='r', linestyle='--', alpha=0.7,
            label=f'Optimal L/D ({ld_ratio:.2f})')
ax3.set_xlabel('Aspect Ratio')
ax3.set_ylabel('Lift-to-Drag Ratio')
ax3.set_title('Aerodynamic Efficiency vs Aspect Ratio')
ax3.grid(True, alpha=0.3)
ax3.legend()

# 4. Comparison with real bird species
species_names = list(bird_species.keys())
species_wingspans = [bird_species[sp]['wingspan'] for sp in species_names]
species_aspect_ratios = [bird_species[sp]['aspect_ratio'] for sp in species_names]

ax4.scatter(species_wingspans, species_aspect_ratios, 
           s=100, c='orange', alpha=0.7, label='Real Bird Species')
ax4.scatter(optimal_wingspan, optimal_aspect_ratio, 
           s=150, c='red', marker='*', label='Optimized Design')

for i, species in enumerate(species_names):
    ax4.annotate(species, (species_wingspans[i], species_aspect_ratios[i]),
                xytext=(5, 5), textcoords='offset points', fontsize=9)

ax4.set_xlabel('Wingspan (m)')
ax4.set_ylabel('Aspect Ratio')
ax4.set_title('Design Comparison: Optimized vs Real Birds')
ax4.grid(True, alpha=0.3)
ax4.legend()

plt.tight_layout()
plt.show()

# Performance analysis across different flight speeds
speeds = np.linspace(5, 20, 30)
ld_ratios_speed = []

for speed in speeds:
    optimizer.flight_speed = speed
    test_lift = optimizer.calculate_lift_coefficient(wing_area)
    test_drag = optimizer.calculate_drag_coefficient(test_lift, optimal_aspect_ratio, oswald_eff)
    ld_ratios_speed.append(test_lift / test_drag)

# Reset to original speed
optimizer.flight_speed = 10.0

plt.figure(figsize=(10, 6))
plt.plot(speeds, ld_ratios_speed, 'b-', linewidth=2, label='L/D Ratio vs Speed')
plt.axvline(x=10.0, color='r', linestyle='--', label='Design Speed (10 m/s)')
plt.xlabel('Flight Speed (m/s)')
plt.ylabel('Lift-to-Drag Ratio')
plt.title('Aerodynamic Performance vs Flight Speed')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

print("\n=== Analysis Complete ===")
print("The optimization has found the ideal wing geometry for maximum aerodynamic efficiency!")

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

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

class MetabolicPathwayOptimizer:
    """
    A class to optimize metabolic pathways using linear and nonlinear programming
    """
    
    def __init__(self, enzyme_capacities, km_values, vmax_values):
        """
        Initialize the optimizer with enzyme parameters
        
        Parameters:
        - enzyme_capacities: Maximum flux through each enzyme
        - km_values: Michaelis-Menten constants for each enzyme
        - vmax_values: Maximum reaction rates for each enzyme
        """
        self.enzyme_capacities = enzyme_capacities
        self.km_values = km_values
        self.vmax_values = vmax_values
        self.results = {}
    
    def linear_optimization(self):
        """
        Solve the linear programming version of the metabolic optimization problem
        """
        # Objective: maximize v3 (ATP production)
        # In scipy.optimize.linprog, we minimize, so we use -v3
        c = [0, 0, -1]  # Coefficients for v1, v2, v3
        
        # Inequality constraints: vi <= vi_max
        A_ub = [[1, 0, 0],    # v1 <= v1_max
                [0, 1, 0],    # v2 <= v2_max
                [0, 0, 1]]    # v3 <= v3_max
        b_ub = self.enzyme_capacities
        
        # Equality constraints: mass balance
        A_eq = [[1, -1, 0],   # v1 - v2 = 0
                [0, 1, -1]]   # v2 - v3 = 0
        b_eq = [0, 0]
        
        # Bounds: all fluxes >= 0
        bounds = [(0, None), (0, None), (0, None)]
        
        # Solve linear program
        result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, 
                        bounds=bounds, method='highs')
        
        self.results['linear'] = {
            'fluxes': result.x,
            'objective': -result.fun,
            'success': result.success
        }
        
        return result.x, -result.fun
    
    def michaelis_menten_rate(self, substrate_conc, vmax, km):
        """
        Calculate reaction rate using Michaelis-Menten kinetics
        """
        return vmax * substrate_conc / (km + substrate_conc)
    
    def nonlinear_optimization(self, initial_concentrations):
        """
        Solve nonlinear optimization with Michaelis-Menten kinetics
        """
        def objective(x):
            # x = [v1, v2, v3, S1, S2, S3] where Si are substrate concentrations
            return -x[2]  # Maximize v3 (ATP production)
        
        def constraint_mm_kinetics(x):
            """
            Michaelis-Menten kinetics constraints
            """
            v1, v2, v3, s1, s2, s3 = x
            
            # Calculate theoretical rates based on substrate concentrations
            v1_theoretical = self.michaelis_menten_rate(s1, self.vmax_values[0], self.km_values[0])
            v2_theoretical = self.michaelis_menten_rate(s2, self.vmax_values[1], self.km_values[1])
            v3_theoretical = self.michaelis_menten_rate(s3, self.vmax_values[2], self.km_values[2])
            
            return np.array([
                v1 - v1_theoretical,
                v2 - v2_theoretical,
                v3 - v3_theoretical
            ])
        
        def constraint_mass_balance(x):
            """
            Mass balance constraints
            """
            v1, v2, v3, s1, s2, s3 = x
            return np.array([
                v1 - v2,  # v1 = v2
                v2 - v3   # v2 = v3
            ])
        
        # Initial guess: [v1, v2, v3, s1, s2, s3]
        x0 = np.array([1.0, 1.0, 1.0] + initial_concentrations)
        
        # Bounds
        bounds = [(0, self.enzyme_capacities[0]),  # v1
                  (0, self.enzyme_capacities[1]),  # v2
                  (0, self.enzyme_capacities[2]),  # v3
                  (0.1, 10),  # s1 (substrate concentrations)
                  (0.1, 10),  # s2
                  (0.1, 10)]  # s3
        
        # Constraints
        constraints = [
            NonlinearConstraint(constraint_mm_kinetics, 0, 0),
            NonlinearConstraint(constraint_mass_balance, 0, 0)
        ]
        
        # Solve nonlinear program
        result = minimize(objective, x0, bounds=bounds, constraints=constraints, 
                         method='SLSQP')
        
        if result.success:
            fluxes = result.x[:3]
            concentrations = result.x[3:]
            self.results['nonlinear'] = {
                'fluxes': fluxes,
                'concentrations': concentrations,
                'objective': -result.fun,
                'success': result.success
            }
        
        return result
    
    def sensitivity_analysis(self, parameter_range=0.2, n_points=20):
        """
        Perform sensitivity analysis on enzyme capacities
        """
        base_capacities = np.array(self.enzyme_capacities)
        sensitivity_results = []
        
        for i in range(len(base_capacities)):
            # Vary the i-th enzyme capacity
            capacity_values = np.linspace(
                base_capacities[i] * (1 - parameter_range),
                base_capacities[i] * (1 + parameter_range),
                n_points
            )
            
            objectives = []
            for cap_val in capacity_values:
                # Create temporary capacities
                temp_capacities = base_capacities.copy()
                temp_capacities[i] = cap_val
                
                # Solve optimization with modified capacity
                temp_optimizer = MetabolicPathwayOptimizer(
                    temp_capacities, self.km_values, self.vmax_values
                )
                _, obj_val = temp_optimizer.linear_optimization()
                objectives.append(obj_val)
            
            sensitivity_results.append({
                'enzyme_index': i,
                'capacity_values': capacity_values,
                'objectives': objectives
            })
        
        self.results['sensitivity'] = sensitivity_results
        return sensitivity_results
    
    def plot_results(self):
        """
        Create comprehensive visualization of optimization results
        """
        fig = plt.figure(figsize=(16, 12))
        
        # Plot 1: Flux distribution comparison
        ax1 = fig.add_subplot(2, 3, 1)
        if 'linear' in self.results and 'nonlinear' in self.results:
            linear_fluxes = self.results['linear']['fluxes']
            nonlinear_fluxes = self.results['nonlinear']['fluxes']
            
            x_pos = np.arange(3)
            width = 0.35
            
            ax1.bar(x_pos - width/2, linear_fluxes, width, 
                   label='Linear Optimization', alpha=0.8)
            ax1.bar(x_pos + width/2, nonlinear_fluxes, width, 
                   label='Nonlinear Optimization', alpha=0.8)
            
            ax1.set_xlabel('Reaction Step')
            ax1.set_ylabel('Flux (mmol/h)')
            ax1.set_title('Flux Distribution Comparison')
            ax1.set_xticks(x_pos)
            ax1.set_xticklabels(['v₁', 'v₂', 'v₃'])
            ax1.legend()
            ax1.grid(True, alpha=0.3)
        
        # Plot 2: Objective function values
        ax2 = fig.add_subplot(2, 3, 2)
        if 'linear' in self.results and 'nonlinear' in self.results:
            methods = ['Linear', 'Nonlinear']
            objectives = [self.results['linear']['objective'], 
                         self.results['nonlinear']['objective']]
            
            bars = ax2.bar(methods, objectives, color=['skyblue', 'lightcoral'])
            ax2.set_ylabel('ATP Production Rate (mmol/h)')
            ax2.set_title('Optimization Results Comparison')
            ax2.grid(True, alpha=0.3)
            
            # Add value labels on bars
            for bar, obj in zip(bars, objectives):
                height = bar.get_height()
                ax2.text(bar.get_x() + bar.get_width()/2., height,
                        f'{obj:.3f}', ha='center', va='bottom')
        
        # Plot 3: Sensitivity analysis
        if 'sensitivity' in self.results:
            ax3 = fig.add_subplot(2, 3, 3)
            sensitivity_data = self.results['sensitivity']
            
            for i, data in enumerate(sensitivity_data):
                ax3.plot(data['capacity_values'], data['objectives'], 
                        marker='o', linewidth=2, markersize=4,
                        label=f'Enzyme {i+1}')
            
            ax3.set_xlabel('Enzyme Capacity (mmol/h)')
            ax3.set_ylabel('ATP Production Rate (mmol/h)')
            ax3.set_title('Sensitivity Analysis')
            ax3.legend()
            ax3.grid(True, alpha=0.3)
        
        # Plot 4: Metabolite concentrations (if available)
        if 'nonlinear' in self.results:
            ax4 = fig.add_subplot(2, 3, 4)
            concentrations = self.results['nonlinear']['concentrations']
            metabolites = ['Glucose', 'G6P', 'Pyruvate']
            
            bars = ax4.bar(metabolites, concentrations, color='lightgreen')
            ax4.set_ylabel('Concentration (mM)')
            ax4.set_title('Optimal Metabolite Concentrations')
            ax4.grid(True, alpha=0.3)
            
            # Add value labels
            for bar, conc in zip(bars, concentrations):
                height = bar.get_height()
                ax4.text(bar.get_x() + bar.get_width()/2., height,
                        f'{conc:.3f}', ha='center', va='bottom')
        
        # Plot 5: Michaelis-Menten curves
        ax5 = fig.add_subplot(2, 3, 5)
        substrate_range = np.linspace(0.1, 10, 100)
        
        for i in range(3):
            rates = [self.michaelis_menten_rate(s, self.vmax_values[i], self.km_values[i]) 
                    for s in substrate_range]
            ax5.plot(substrate_range, rates, linewidth=2, 
                    label=f'Enzyme {i+1}')
        
        ax5.set_xlabel('Substrate Concentration (mM)')
        ax5.set_ylabel('Reaction Rate (mmol/h)')
        ax5.set_title('Michaelis-Menten Kinetics')
        ax5.legend()
        ax5.grid(True, alpha=0.3)
        
        # Plot 6: 3D optimization landscape (simplified)
        ax6 = fig.add_subplot(2, 3, 6, projection='3d')
        
        # Create a simplified 3D surface for visualization
        v1_range = np.linspace(0, min(self.enzyme_capacities), 20)
        v2_range = np.linspace(0, min(self.enzyme_capacities), 20)
        V1, V2 = np.meshgrid(v1_range, v2_range)
        
        # For mass balance, v3 = v2 = v1, so objective = min(v1, v2, capacity)
        Z = np.minimum(np.minimum(V1, V2), min(self.enzyme_capacities))
        
        surf = ax6.plot_surface(V1, V2, Z, alpha=0.7, cmap='viridis')
        ax6.set_xlabel('v₁ (mmol/h)')
        ax6.set_ylabel('v₂ (mmol/h)')
        ax6.set_zlabel('ATP Production (mmol/h)')
        ax6.set_title('Optimization Landscape')
        
        plt.tight_layout()
        plt.show()
        
        return fig

# Example usage and demonstration
def main():
    """
    Main function to demonstrate metabolic pathway optimization
    """
    print("=== Metabolic Pathway Optimization Demo ===\n")
    
    # Define enzyme parameters
    enzyme_capacities = [5.0, 3.0, 4.0]  # mmol/h
    km_values = [0.5, 1.0, 0.8]  # mM
    vmax_values = [6.0, 4.0, 5.0]  # mmol/h
    initial_concentrations = [2.0, 1.5, 1.0]  # mM
    
    # Create optimizer
    optimizer = MetabolicPathwayOptimizer(enzyme_capacities, km_values, vmax_values)
    
    # Perform linear optimization
    print("1. Linear Optimization Results:")
    linear_fluxes, linear_objective = optimizer.linear_optimization()
    print(f"   Optimal fluxes: v₁={linear_fluxes[0]:.3f}, v₂={linear_fluxes[1]:.3f}, v₃={linear_fluxes[2]:.3f}")
    print(f"   Maximum ATP production: {linear_objective:.3f} mmol/h")
    print(f"   Success: {optimizer.results['linear']['success']}\n")
    
    # Perform nonlinear optimization
    print("2. Nonlinear Optimization Results:")
    nonlinear_result = optimizer.nonlinear_optimization(initial_concentrations)
    
    if nonlinear_result.success:
        nl_fluxes = optimizer.results['nonlinear']['fluxes']
        nl_concentrations = optimizer.results['nonlinear']['concentrations']
        nl_objective = optimizer.results['nonlinear']['objective']
        
        print(f"   Optimal fluxes: v₁={nl_fluxes[0]:.3f}, v₂={nl_fluxes[1]:.3f}, v₃={nl_fluxes[2]:.3f}")
        print(f"   Optimal concentrations: [Glucose]={nl_concentrations[0]:.3f}, [G6P]={nl_concentrations[1]:.3f}, [Pyruvate]={nl_concentrations[2]:.3f}")
        print(f"   Maximum ATP production: {nl_objective:.3f} mmol/h")
        print(f"   Success: {optimizer.results['nonlinear']['success']}\n")
    else:
        print("   Nonlinear optimization failed to converge\n")
    
    # Perform sensitivity analysis
    print("3. Sensitivity Analysis:")
    sensitivity_results = optimizer.sensitivity_analysis()
    
    for i, result in enumerate(sensitivity_results):
        max_obj = max(result['objectives'])
        min_obj = min(result['objectives'])
        sensitivity = (max_obj - min_obj) / min_obj * 100
        print(f"   Enzyme {i+1} sensitivity: {sensitivity:.1f}% change in objective")
    
    print("\n4. Generating comprehensive visualization...")
    
    # Generate plots
    fig = optimizer.plot_results()
    
    # Print summary statistics
    print("\n=== Summary Statistics ===")
    print(f"Linear optimization ATP yield: {linear_objective:.3f} mmol/h")
    if 'nonlinear' in optimizer.results:
        print(f"Nonlinear optimization ATP yield: {optimizer.results['nonlinear']['objective']:.3f} mmol/h")
    
    print(f"Bottleneck enzyme capacity: {min(enzyme_capacities):.1f} mmol/h")
    print(f"Theoretical maximum flux: {min(enzyme_capacities):.1f} mmol/h")
    
    return optimizer

# Run the demonstration
if __name__ == "__main__":
    optimizer = main()