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

c = [0, 0, -1]  # Minimize -v3 (maximize v3)

capacity_values = np.linspace(
    base_capacities[i] * (1 - parameter_range),
    base_capacities[i] * (1 + parameter_range),
    n_points
)

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

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

class ForagingOptimizer:
    def __init__(self, R_A, h_A, R_B, h_B, T_travel):
        """
        Initialize the foraging model parameters
        
        Parameters:
        R_A, R_B: Maximum energy gain rates for patches A and B
        h_A, h_B: Handling time parameters for patches A and B
        T_travel: Travel time between patches
        """
        self.R_A = R_A
        self.h_A = h_A
        self.R_B = R_B
        self.h_B = h_B
        self.T_travel = T_travel
    
    def energy_gain_patch_A(self, t_A):
        """Calculate energy gain from patch A"""
        return (self.R_A * t_A) / (1 + self.h_A * t_A)
    
    def energy_gain_patch_B(self, t_B):
        """Calculate energy gain from patch B"""
        return (self.R_B * t_B) / (1 + self.h_B * t_B)
    
    def total_energy_rate(self, times):
        """Calculate total energy gain rate"""
        t_A, t_B = times
        if t_A < 0 or t_B < 0:
            return -np.inf
        
        total_energy = self.energy_gain_patch_A(t_A) + self.energy_gain_patch_B(t_B)
        total_time = t_A + t_B + self.T_travel
        
        return total_energy / total_time
    
    def negative_energy_rate(self, times):
        """Negative energy rate for minimization"""
        return -self.total_energy_rate(times)
    
    def marginal_gain_rate_A(self, t_A):
        """Calculate marginal gain rate for patch A"""
        return self.R_A / ((1 + self.h_A * t_A) ** 2)
    
    def marginal_gain_rate_B(self, t_B):
        """Calculate marginal gain rate for patch B"""
        return self.R_B / ((1 + self.h_B * t_B) ** 2)
    
    def optimize_foraging(self):
        """Find optimal foraging times using scipy optimization"""
        # Initial guess
        initial_guess = [5.0, 5.0]
        
        # Bounds to ensure positive times
        bounds = [(0.1, 50.0), (0.1, 50.0)]
        
        # Optimize
        result = minimize(self.negative_energy_rate, initial_guess, 
                         bounds=bounds, method='L-BFGS-B')
        
        optimal_times = result.x
        optimal_rate = -result.fun
        
        return optimal_times, optimal_rate

# Set up the problem parameters
# Patch A: High-quality food (berries)
R_A = 15.0  # Maximum energy gain rate (energy units/time)
h_A = 0.8   # Handling time parameter

# Patch B: Lower-quality food (seeds)
R_B = 8.0   # Maximum energy gain rate
h_B = 0.3   # Handling time parameter

# Travel time between patches
T_travel = 2.0

# Create the foraging optimizer
forager = ForagingOptimizer(R_A, h_A, R_B, h_B, T_travel)

# Find optimal solution
optimal_times, optimal_rate = forager.optimize_foraging()
t_A_opt, t_B_opt = optimal_times

print("=== FORAGING STRATEGY OPTIMIZATION RESULTS ===")
print(f"Optimal time in Patch A (berries): {t_A_opt:.2f} time units")
print(f"Optimal time in Patch B (seeds): {t_B_opt:.2f} time units")
print(f"Optimal energy gain rate: {optimal_rate:.3f} energy units/time")
print(f"Total foraging cycle time: {t_A_opt + t_B_opt + T_travel:.2f} time units")

# Calculate energy gains at optimal times
energy_A_opt = forager.energy_gain_patch_A(t_A_opt)
energy_B_opt = forager.energy_gain_patch_B(t_B_opt)

print(f"\nEnergy gained from Patch A: {energy_A_opt:.3f} energy units")
print(f"Energy gained from Patch B: {energy_B_opt:.3f} energy units")
print(f"Total energy gained: {energy_A_opt + energy_B_opt:.3f} energy units")

# Calculate marginal gain rates at optimal times
marginal_A = forager.marginal_gain_rate_A(t_A_opt)
marginal_B = forager.marginal_gain_rate_B(t_B_opt)

print(f"\nMarginal gain rate in Patch A: {marginal_A:.6f}")
print(f"Marginal gain rate in Patch B: {marginal_B:.6f}")
print(f"Difference (should be close to 0): {abs(marginal_A - marginal_B):.6f}")

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

# Plot 1: Energy gain vs time for both patches
t_range = np.linspace(0.1, 20, 1000)
energy_A = [forager.energy_gain_patch_A(t) for t in t_range]
energy_B = [forager.energy_gain_patch_B(t) for t in t_range]

ax1.plot(t_range, energy_A, 'b-', linewidth=2, label='Patch A (Berries)', alpha=0.8)
ax1.plot(t_range, energy_B, 'r-', linewidth=2, label='Patch B (Seeds)', alpha=0.8)
ax1.axvline(t_A_opt, color='blue', linestyle='--', alpha=0.7, label=f'Optimal t_A = {t_A_opt:.1f}')
ax1.axvline(t_B_opt, color='red', linestyle='--', alpha=0.7, label=f'Optimal t_B = {t_B_opt:.1f}')
ax1.set_xlabel('Time in Patch')
ax1.set_ylabel('Energy Gain')
ax1.set_title('Energy Gain vs Time in Each Patch')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Marginal gain rates
marginal_A_range = [forager.marginal_gain_rate_A(t) for t in t_range]
marginal_B_range = [forager.marginal_gain_rate_B(t) for t in t_range]

ax2.plot(t_range, marginal_A_range, 'b-', linewidth=2, label='Patch A Marginal Rate', alpha=0.8)
ax2.plot(t_range, marginal_B_range, 'r-', linewidth=2, label='Patch B Marginal Rate', alpha=0.8)
ax2.axhline(marginal_A, color='green', linestyle=':', linewidth=2, label=f'Optimal Marginal Rate = {marginal_A:.3f}')
ax2.axvline(t_A_opt, color='blue', linestyle='--', alpha=0.7)
ax2.axvline(t_B_opt, color='red', linestyle='--', alpha=0.7)
ax2.set_xlabel('Time in Patch')
ax2.set_ylabel('Marginal Gain Rate')
ax2.set_title('Marginal Gain Rates vs Time')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: 3D surface plot of total energy rate
t_A_range = np.linspace(0.1, 15, 50)
t_B_range = np.linspace(0.1, 15, 50)
T_A, T_B = np.meshgrid(t_A_range, t_B_range)

# Calculate total energy rate for all combinations
Z = np.zeros_like(T_A)
for i in range(len(t_A_range)):
    for j in range(len(t_B_range)):
        Z[j, i] = forager.total_energy_rate([T_A[j, i], T_B[j, i]])

contour = ax3.contour(T_A, T_B, Z, levels=15, alpha=0.7)
ax3.clabel(contour, inline=True, fontsize=8)
ax3.plot(t_A_opt, t_B_opt, 'go', markersize=10, label=f'Optimal Point ({t_A_opt:.1f}, {t_B_opt:.1f})')
ax3.set_xlabel('Time in Patch A')
ax3.set_ylabel('Time in Patch B')
ax3.set_title('Total Energy Rate Contours')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Sensitivity analysis
travel_times = np.linspace(0.5, 10, 20)
optimal_rates = []
optimal_t_A = []
optimal_t_B = []

for T_travel_test in travel_times:
    forager_test = ForagingOptimizer(R_A, h_A, R_B, h_B, T_travel_test)
    opt_times, opt_rate = forager_test.optimize_foraging()
    optimal_rates.append(opt_rate)
    optimal_t_A.append(opt_times[0])
    optimal_t_B.append(opt_times[1])

ax4.plot(travel_times, optimal_rates, 'g-', linewidth=2, label='Optimal Energy Rate')
ax4.axvline(T_travel, color='orange', linestyle='--', label=f'Current Travel Time = {T_travel}')
ax4.set_xlabel('Travel Time Between Patches')
ax4.set_ylabel('Optimal Energy Rate')
ax4.set_title('Sensitivity to Travel Time')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Compare strategies
print("\n=== STRATEGY COMPARISON ===")
strategies = {
    'Optimal Strategy': (t_A_opt, t_B_opt),
    'Equal Time Strategy': (5.0, 5.0),
    'Patch A Only': (10.0, 0.1),
    'Patch B Only': (0.1, 10.0)
}

print(f"{'Strategy':<20} {'t_A':<8} {'t_B':<8} {'Energy Rate':<12} {'Efficiency':<12}")
print("-" * 65)

for name, (t_A, t_B) in strategies.items():
    rate = forager.total_energy_rate([t_A, t_B])
    efficiency = (rate / optimal_rate) * 100
    print(f"{name:<20} {t_A:<8.1f} {t_B:<8.1f} {rate:<12.3f} {efficiency:<12.1f}%")

# Mathematical verification
print("\n=== MATHEMATICAL VERIFICATION ===")
print("According to Marginal Value Theorem, at optimum:")
print("Marginal gain rate in Patch A = Marginal gain rate in Patch B")
print(f"Left side: dE_A/dt_A = R_A / (1 + h_A * t_A)^2 = {marginal_A:.6f}")
print(f"Right side: dE_B/dt_B = R_B / (1 + h_B * t_B)^2 = {marginal_B:.6f}")
print(f"Verification: |Left - Right| = {abs(marginal_A - marginal_B):.8f} ≈ 0 ✓")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd
from scipy import constants

# Set up plotting parameters
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 12

# Physical constants
G = constants.G  # m³/kg/s²
AU = constants.au  # meters
R_sun = 6.96e8  # meters
M_sun = 1.989e30  # kg
R_earth = 6.371e6  # meters
R_jupiter = 7.1492e7  # meters

print("=== Exoplanet Parameter Estimation using Transit Photometry ===\n")

# STEP 1: Define our example system parameters (these would be "unknown" in reality)
print("1. TRUE SYSTEM PARAMETERS (what we're trying to estimate):")
print("-" * 60)

# Star properties
stellar_mass = 1.1 * M_sun  # kg
stellar_radius = 1.05 * R_sun  # meters
stellar_temp = 5800  # K

# Planet properties  
planet_radius = 1.2 * R_earth  # meters
orbital_period = 3.5 * 24 * 3600  # seconds (3.5 days)
impact_parameter = 0.3  # how close to center the planet passes

print(f"Stellar mass: {stellar_mass/M_sun:.2f} M☉")
print(f"Stellar radius: {stellar_radius/R_sun:.2f} R☉")
print(f"Planet radius: {planet_radius/R_earth:.2f} R⊕")
print(f"Orbital period: {orbital_period/(24*3600):.2f} days")
print(f"Impact parameter: {impact_parameter:.2f}")

# Calculate derived parameters
semi_major_axis = ((G * stellar_mass * orbital_period**2) / (4 * np.pi**2))**(1/3)
radius_ratio = planet_radius / stellar_radius
transit_depth = radius_ratio**2

print(f"Semi-major axis: {semi_major_axis/AU:.4f} AU")
print(f"Radius ratio (Rp/R*): {radius_ratio:.4f}")
print(f"Transit depth: {transit_depth:.6f} ({transit_depth*1e6:.1f} ppm)")

# STEP 2: Calculate transit duration
print(f"\n2. TRANSIT DURATION CALCULATION:")
print("-" * 40)

# Transit duration formula
def calculate_transit_duration(period, stellar_radius, semi_major_axis, radius_ratio, impact_param):
    """Calculate transit duration using the standard formula"""
    term = (stellar_radius / semi_major_axis) * np.sqrt((1 + radius_ratio)**2 - impact_param**2)
    duration = (period / np.pi) * np.arcsin(term)
    return duration

transit_duration = calculate_transit_duration(orbital_period, stellar_radius, 
                                            semi_major_axis, radius_ratio, impact_parameter)

print(f"Transit duration: {transit_duration/3600:.2f} hours")

# STEP 3: Generate synthetic transit light curve
print(f"\n3. GENERATING SYNTHETIC TRANSIT DATA:")
print("-" * 45)

def transit_model(time, period, t0, depth, duration, baseline=1.0):
    """
    Simple box-shaped transit model
    
    Parameters:
    - time: array of time values
    - period: orbital period
    - t0: transit center time
    - depth: transit depth (fractional)
    - duration: transit duration
    - baseline: out-of-transit flux level
    """
    # Calculate phase
    phase = ((time - t0) % period) / period
    # Center phase around transit (phase = 0.5 means transit center)
    phase = np.where(phase > 0.5, phase - 1, phase)
    
    # Simple box model
    flux = np.full_like(time, baseline)
    in_transit = np.abs(phase * period) < duration / 2
    flux[in_transit] = baseline * (1 - depth)
    
    return flux

# Generate time series (observe for 10 days)
observation_time = 10 * 24 * 3600  # 10 days in seconds
time_points = np.linspace(0, observation_time, 10000)

# Add some random noise to make it realistic
np.random.seed(42)  # for reproducibility
noise_level = 0.001  # 0.1% photometric precision
synthetic_flux = transit_model(time_points, orbital_period, orbital_period/2, 
                              transit_depth, transit_duration) + np.random.normal(0, noise_level, len(time_points))

print(f"Generated {len(time_points)} data points over {observation_time/(24*3600):.1f} days")
print(f"Photometric precision: {noise_level*100:.1f}%")
print(f"Expected number of transits: {observation_time/orbital_period:.1f}")

# STEP 4: Parameter estimation using curve fitting
print(f"\n4. PARAMETER ESTIMATION:")
print("-" * 30)

# Define fitting function
def transit_fit_function(time, period, t0, depth, duration):
    return transit_model(time, period, t0, depth, duration, baseline=1.0)

# Initial parameter guesses (what we might estimate from the data)
initial_guess = [
    3.0 * 24 * 3600,  # period guess (3 days)
    1.5 * 24 * 3600,  # t0 guess
    0.001,            # depth guess
    2.0 * 3600        # duration guess (2 hours)
]

print("Initial parameter guesses:")
print(f"Period: {initial_guess[0]/(24*3600):.2f} days")
print(f"Transit center: {initial_guess[1]/(24*3600):.2f} days")
print(f"Depth: {initial_guess[2]:.6f} ({initial_guess[2]*1e6:.1f} ppm)")
print(f"Duration: {initial_guess[3]/3600:.2f} hours")

# Perform curve fitting
try:
    fitted_params, param_covariance = curve_fit(
        transit_fit_function, time_points, synthetic_flux, 
        p0=initial_guess, maxfev=5000
    )
    
    # Calculate parameter uncertainties
    param_errors = np.sqrt(np.diag(param_covariance))
    
    print(f"\nFitted parameters:")
    print(f"Period: {fitted_params[0]/(24*3600):.3f} ± {param_errors[0]/(24*3600):.3f} days")
    print(f"Transit center: {fitted_params[1]/(24*3600):.3f} ± {param_errors[1]/(24*3600):.3f} days")
    print(f"Depth: {fitted_params[2]:.6f} ± {param_errors[2]:.6f} ({fitted_params[2]*1e6:.1f} ± {param_errors[2]*1e6:.1f} ppm)")
    print(f"Duration: {fitted_params[3]/3600:.3f} ± {param_errors[3]/3600:.3f} hours")
    
    # Calculate derived parameters
    fitted_radius_ratio = np.sqrt(fitted_params[2])
    fitted_planet_radius = fitted_radius_ratio * stellar_radius
    fitted_semi_major_axis = ((G * stellar_mass * fitted_params[0]**2) / (4 * np.pi**2))**(1/3)
    
    print(f"\nDerived parameters:")
    print(f"Radius ratio: {fitted_radius_ratio:.4f}")
    print(f"Planet radius: {fitted_planet_radius/R_earth:.2f} R⊕")
    print(f"Semi-major axis: {fitted_semi_major_axis/AU:.4f} AU")
    
    # Calculate accuracy
    period_accuracy = abs(fitted_params[0] - orbital_period) / orbital_period * 100
    depth_accuracy = abs(fitted_params[2] - transit_depth) / transit_depth * 100
    duration_accuracy = abs(fitted_params[3] - transit_duration) / transit_duration * 100
    
    print(f"\nFitting accuracy:")
    print(f"Period error: {period_accuracy:.2f}%")
    print(f"Depth error: {depth_accuracy:.2f}%")
    print(f"Duration error: {duration_accuracy:.2f}%")
    
    fitting_successful = True
    
except Exception as e:
    print(f"Fitting failed: {e}")
    fitting_successful = False

# STEP 5: Create comprehensive plots
print(f"\n5. CREATING VISUALIZATION PLOTS:")
print("-" * 40)

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

# Plot 1: Full light curve
ax1.plot(time_points/(24*3600), synthetic_flux, 'b.', alpha=0.3, markersize=1, label='Synthetic data')
if fitting_successful:
    model_flux = transit_fit_function(time_points, *fitted_params)
    ax1.plot(time_points/(24*3600), model_flux, 'r-', linewidth=2, label='Fitted model')
ax1.set_xlabel('Time (days)')
ax1.set_ylabel('Relative Flux')
ax1.set_title('Full Transit Light Curve')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Zoomed view of single transit
# Find the first transit
first_transit_center = orbital_period/2
plot_window = 4 * transit_duration
mask = np.abs(time_points - first_transit_center) < plot_window
time_zoom = time_points[mask]
flux_zoom = synthetic_flux[mask]

ax2.plot((time_zoom - first_transit_center)/3600, flux_zoom, 'b.', alpha=0.6, markersize=2, label='Data')
if fitting_successful:
    model_zoom = transit_fit_function(time_zoom, *fitted_params)
    ax2.plot((time_zoom - first_transit_center)/3600, model_zoom, 'r-', linewidth=2, label='Model')
ax2.set_xlabel('Time from transit center (hours)')
ax2.set_ylabel('Relative Flux')
ax2.set_title('Single Transit (Zoomed)')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Phase-folded light curve
phase = ((time_points - orbital_period/2) % orbital_period) / orbital_period
phase = np.where(phase > 0.5, phase - 1, phase)

ax3.plot(phase, synthetic_flux, 'b.', alpha=0.3, markersize=1, label='Data')
if fitting_successful:
    phase_model = ((time_points - fitted_params[1]) % fitted_params[0]) / fitted_params[0]
    phase_model = np.where(phase_model > 0.5, phase_model - 1, phase_model)
    sort_idx = np.argsort(phase_model)
    model_flux_all = transit_fit_function(time_points, *fitted_params)
    ax3.plot(phase_model[sort_idx], model_flux_all[sort_idx], 'r-', linewidth=2, label='Model')
ax3.set_xlabel('Orbital Phase')
ax3.set_ylabel('Relative Flux')
ax3.set_title('Phase-Folded Light Curve')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-0.1, 0.1)

# Plot 4: Parameter comparison
if fitting_successful:
    true_params = [orbital_period/(24*3600), transit_depth*1e6, transit_duration/3600, 
                   planet_radius/R_earth]
    fitted_params_plot = [fitted_params[0]/(24*3600), fitted_params[2]*1e6, 
                         fitted_params[3]/3600, fitted_planet_radius/R_earth]
    param_names = ['Period (days)', 'Depth (ppm)', 'Duration (hrs)', 'Planet Radius (R⊕)']
    
    x = np.arange(len(param_names))
    width = 0.35
    
    bars1 = ax4.bar(x - width/2, true_params, width, label='True values', alpha=0.7)
    bars2 = ax4.bar(x + width/2, fitted_params_plot, width, label='Fitted values', alpha=0.7)
    
    ax4.set_xlabel('Parameters')
    ax4.set_ylabel('Value')
    ax4.set_title('Parameter Comparison')
    ax4.set_xticks(x)
    ax4.set_xticklabels(param_names, rotation=45, ha='right')
    ax4.legend()
    ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# STEP 6: Statistical analysis
print(f"\n6. STATISTICAL ANALYSIS:")
print("-" * 30)

if fitting_successful:
    # Chi-squared goodness of fit
    model_flux_all = transit_fit_function(time_points, *fitted_params)
    residuals = synthetic_flux - model_flux_all
    chi_squared = np.sum((residuals / noise_level)**2)
    reduced_chi_squared = chi_squared / (len(time_points) - len(fitted_params))
    
    print(f"Chi-squared: {chi_squared:.2f}")
    print(f"Reduced chi-squared: {reduced_chi_squared:.2f}")
    print(f"RMS residuals: {np.std(residuals)*1e6:.1f} ppm")
    
    # Signal-to-noise ratio
    transit_snr = (transit_depth / noise_level) * np.sqrt(np.sum(np.abs(phase) < 0.02))
    print(f"Transit signal-to-noise ratio: {transit_snr:.1f}")

# STEP 7: Summary table
print(f"\n7. RESULTS SUMMARY:")
print("-" * 25)

if fitting_successful:
    summary_data = {
        'Parameter': ['Orbital Period', 'Transit Depth', 'Transit Duration', 'Planet Radius', 'Semi-major Axis'],
        'True Value': [f"{orbital_period/(24*3600):.3f} days",
                      f"{transit_depth*1e6:.1f} ppm",
                      f"{transit_duration/3600:.2f} hours",
                      f"{planet_radius/R_earth:.2f} R⊕",
                      f"{semi_major_axis/AU:.4f} AU"],
        'Fitted Value': [f"{fitted_params[0]/(24*3600):.3f} days",
                        f"{fitted_params[2]*1e6:.1f} ppm",
                        f"{fitted_params[3]/3600:.2f} hours",
                        f"{fitted_planet_radius/R_earth:.2f} R⊕",
                        f"{fitted_semi_major_axis/AU:.4f} AU"],
        'Uncertainty': [f"±{param_errors[0]/(24*3600):.3f} days",
                       f"±{param_errors[2]*1e6:.1f} ppm",
                       f"±{param_errors[3]/3600:.2f} hours",
                       f"±{fitted_radius_ratio*param_errors[2]/(2*fitted_params[2])*stellar_radius/R_earth:.2f} R⊕",
                       "N/A"]
    }
    
    summary_df = pd.DataFrame(summary_data)
    print(summary_df.to_string(index=False))

print(f"\n" + "="*70)
print("ANALYSIS COMPLETE!")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.stats import chi2
import warnings
warnings.filterwarnings('ignore')

# Set up matplotlib for better plots
plt.style.use('default')
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 12

class ExoplanetAnalysis:
    def __init__(self):
        """Initialize the exoplanet analysis class"""
        self.G = 6.67430e-11  # Gravitational constant (m^3 kg^-1 s^-2)
        self.M_sun = 1.989e30  # Solar mass (kg)
        self.R_sun = 6.96e8   # Solar radius (m)
        self.au = 1.496e11    # Astronomical unit (m)
        
    def generate_synthetic_data(self, seed=42):
        """Generate synthetic transit and RV data for testing"""
        np.random.seed(seed)
        
        # True parameters for our synthetic planet
        self.true_params = {
            'period': 3.52,      # days
            'transit_depth': 0.008,  # fractional depth
            'transit_duration': 0.12,  # days
            'rv_amplitude': 45.0,    # m/s
            'systemic_velocity': 15.0,  # m/s
            'phase_offset': 0.25    # radians
        }
        
        # Generate transit data
        self.transit_time = np.linspace(-0.3, 0.3, 200)  # days around transit
        self.transit_flux = self.transit_model(self.transit_time, 
                                              self.true_params['transit_depth'],
                                              self.true_params['transit_duration'])
        # Add realistic noise
        self.transit_flux += np.random.normal(0, 0.0005, len(self.transit_flux))
        self.transit_error = np.full_like(self.transit_flux, 0.0005)
        
        # Generate RV data
        self.rv_time = np.linspace(0, 14, 25)  # 4 orbital periods
        self.rv_velocity = self.rv_model(self.rv_time,
                                        self.true_params['period'],
                                        self.true_params['rv_amplitude'],
                                        self.true_params['systemic_velocity'],
                                        self.true_params['phase_offset'])
        # Add realistic noise
        self.rv_velocity += np.random.normal(0, 3.0, len(self.rv_velocity))
        self.rv_error = np.full_like(self.rv_velocity, 3.0)
        
        print("Synthetic data generated with true parameters:")
        for key, value in self.true_params.items():
            print(f"  {key}: {value}")
    
    def transit_model(self, t, depth, duration):
        """Simple box-shaped transit model"""
        model = np.ones_like(t)
        in_transit = np.abs(t) < duration/2
        model[in_transit] = 1 - depth
        return model
    
    def rv_model(self, t, period, amplitude, systemic, phase):
        """Sinusoidal radial velocity model"""
        return amplitude * np.sin(2*np.pi*t/period + phase) + systemic
    
    def chi_squared_transit(self, params):
        """Calculate chi-squared for transit data"""
        depth, duration = params
        model = self.transit_model(self.transit_time, depth, duration)
        chi2 = np.sum((self.transit_flux - model)**2 / self.transit_error**2)
        return chi2
    
    def chi_squared_rv(self, params):
        """Calculate chi-squared for RV data"""
        period, amplitude, systemic, phase = params
        model = self.rv_model(self.rv_time, period, amplitude, systemic, phase)
        chi2 = np.sum((self.rv_velocity - model)**2 / self.rv_error**2)
        return chi2
    
    def fit_transit_data(self):
        """Fit transit parameters using chi-squared minimization"""
        print("\n=== TRANSIT DATA FITTING ===")
        
        # Initial guess
        initial_guess = [0.01, 0.1]  # depth, duration
        
        # Bounds for parameters
        bounds = [(0.001, 0.05), (0.05, 0.3)]
        
        # Perform optimization
        result = minimize(self.chi_squared_transit, initial_guess, 
                         bounds=bounds, method='L-BFGS-B')
        
        self.transit_best_params = result.x
        self.transit_chi2 = result.fun
        
        # Calculate degrees of freedom and reduced chi-squared
        dof_transit = len(self.transit_time) - len(initial_guess)
        reduced_chi2_transit = self.transit_chi2 / dof_transit
        
        print(f"Best-fit transit parameters:")
        print(f"  Depth: {self.transit_best_params[0]:.6f} (true: {self.true_params['transit_depth']:.6f})")
        print(f"  Duration: {self.transit_best_params[1]:.6f} days (true: {self.true_params['transit_duration']:.6f})")
        print(f"  Chi-squared: {self.transit_chi2:.2f}")
        print(f"  Reduced chi-squared: {reduced_chi2_transit:.3f}")
        
        return result
    
    def fit_rv_data(self):
        """Fit RV parameters using chi-squared minimization"""
        print("\n=== RADIAL VELOCITY DATA FITTING ===")
        
        # Initial guess
        initial_guess = [3.5, 40.0, 10.0, 0.0]  # period, amplitude, systemic, phase
        
        # Bounds for parameters
        bounds = [(2.0, 5.0), (10.0, 100.0), (-50.0, 50.0), (-np.pi, np.pi)]
        
        # Perform optimization
        result = minimize(self.chi_squared_rv, initial_guess, 
                         bounds=bounds, method='L-BFGS-B')
        
        self.rv_best_params = result.x
        self.rv_chi2 = result.fun
        
        # Calculate degrees of freedom and reduced chi-squared
        dof_rv = len(self.rv_time) - len(initial_guess)
        reduced_chi2_rv = self.rv_chi2 / dof_rv
        
        print(f"Best-fit RV parameters:")
        print(f"  Period: {self.rv_best_params[0]:.6f} days (true: {self.true_params['period']:.6f})")
        print(f"  Amplitude: {self.rv_best_params[1]:.6f} m/s (true: {self.true_params['rv_amplitude']:.6f})")
        print(f"  Systemic velocity: {self.rv_best_params[2]:.6f} m/s (true: {self.true_params['systemic_velocity']:.6f})")
        print(f"  Phase offset: {self.rv_best_params[3]:.6f} rad (true: {self.true_params['phase_offset']:.6f})")
        print(f"  Chi-squared: {self.rv_chi2:.2f}")
        print(f"  Reduced chi-squared: {reduced_chi2_rv:.3f}")
        
        return result
    
    def calculate_physical_parameters(self):
        """Calculate physical parameters from best-fit values"""
        print("\n=== DERIVED PHYSICAL PARAMETERS ===")
        
        # Assume stellar parameters (typical for Sun-like star)
        M_star = 1.0 * self.M_sun  # kg
        R_star = 1.0 * self.R_sun  # m
        
        # Planet radius from transit depth
        R_planet = np.sqrt(self.transit_best_params[0]) * R_star
        R_planet_earth = R_planet / 6.371e6  # Convert to Earth radii
        
        # Semi-major axis from Kepler's third law
        P_seconds = self.rv_best_params[0] * 24 * 3600  # Convert to seconds
        a = ((self.G * M_star * P_seconds**2) / (4 * np.pi**2))**(1/3)
        a_au = a / self.au  # Convert to AU
        
        # Planet mass from RV amplitude (assuming i = 90°)
        K = self.rv_best_params[1]  # m/s
        M_planet = K * np.sqrt(1 - 0**2) * (M_star)**(2/3) * (P_seconds/(2*np.pi*self.G))**(1/3)
        M_planet_earth = M_planet / 5.972e24  # Convert to Earth masses
        
        # Planet density
        V_planet = (4/3) * np.pi * R_planet**3
        rho_planet = M_planet / V_planet
        rho_earth = 5515  # kg/m³
        rho_relative = rho_planet / rho_earth
        
        print(f"Planet radius: {R_planet_earth:.3f} Earth radii")
        print(f"Planet mass: {M_planet_earth:.3f} Earth masses")
        print(f"Semi-major axis: {a_au:.4f} AU")
        print(f"Planet density: {rho_relative:.3f} × Earth density")
        
        # Store results
        self.physical_params = {
            'radius_earth': R_planet_earth,
            'mass_earth': M_planet_earth,
            'semimajor_axis_au': a_au,
            'density_relative': rho_relative
        }
    
    def plot_results(self):
        """Create comprehensive plots of the fits"""
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        
        # Plot 1: Transit light curve
        ax1.errorbar(self.transit_time, self.transit_flux, yerr=self.transit_error,
                    fmt='o', color='blue', alpha=0.7, markersize=4, label='Data')
        
        # Best-fit model
        time_fine = np.linspace(self.transit_time.min(), self.transit_time.max(), 1000)
        model_fine = self.transit_model(time_fine, *self.transit_best_params)
        ax1.plot(time_fine, model_fine, 'r-', linewidth=2, label='Best fit')
        
        # True model
        true_model = self.transit_model(time_fine, 
                                       self.true_params['transit_depth'],
                                       self.true_params['transit_duration'])
        ax1.plot(time_fine, true_model, 'g--', linewidth=2, label='True model')
        
        ax1.set_xlabel('Time from transit center (days)')
        ax1.set_ylabel('Relative flux')
        ax1.set_title('Transit Light Curve Fit')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # Plot 2: Residuals for transit
        model_data = self.transit_model(self.transit_time, *self.transit_best_params)
        residuals = (self.transit_flux - model_data) / self.transit_error
        ax2.errorbar(self.transit_time, residuals, yerr=1.0,
                    fmt='o', color='red', alpha=0.7, markersize=4)
        ax2.axhline(y=0, color='black', linestyle='-', alpha=0.5)
        ax2.axhline(y=2, color='gray', linestyle='--', alpha=0.5)
        ax2.axhline(y=-2, color='gray', linestyle='--', alpha=0.5)
        ax2.set_xlabel('Time from transit center (days)')
        ax2.set_ylabel('Normalized residuals (σ)')
        ax2.set_title('Transit Fit Residuals')
        ax2.grid(True, alpha=0.3)
        
        # Plot 3: RV curve
        ax3.errorbar(self.rv_time, self.rv_velocity, yerr=self.rv_error,
                    fmt='o', color='blue', alpha=0.7, markersize=6, label='Data')
        
        # Best-fit model
        time_fine_rv = np.linspace(self.rv_time.min(), self.rv_time.max(), 1000)
        model_fine_rv = self.rv_model(time_fine_rv, *self.rv_best_params)
        ax3.plot(time_fine_rv, model_fine_rv, 'r-', linewidth=2, label='Best fit')
        
        # True model
        true_model_rv = self.rv_model(time_fine_rv, 
                                     self.true_params['period'],
                                     self.true_params['rv_amplitude'],
                                     self.true_params['systemic_velocity'],
                                     self.true_params['phase_offset'])
        ax3.plot(time_fine_rv, true_model_rv, 'g--', linewidth=2, label='True model')
        
        ax3.set_xlabel('Time (days)')
        ax3.set_ylabel('Radial velocity (m/s)')
        ax3.set_title('Radial Velocity Curve Fit')
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # Plot 4: RV residuals
        model_data_rv = self.rv_model(self.rv_time, *self.rv_best_params)
        residuals_rv = (self.rv_velocity - model_data_rv) / self.rv_error
        ax4.errorbar(self.rv_time, residuals_rv, yerr=1.0,
                    fmt='o', color='red', alpha=0.7, markersize=6)
        ax4.axhline(y=0, color='black', linestyle='-', alpha=0.5)
        ax4.axhline(y=2, color='gray', linestyle='--', alpha=0.5)
        ax4.axhline(y=-2, color='gray', linestyle='--', alpha=0.5)
        ax4.set_xlabel('Time (days)')
        ax4.set_ylabel('Normalized residuals (σ)')
        ax4.set_title('RV Fit Residuals')
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        # Additional plot: Parameter comparison
        self.plot_parameter_comparison()
    
    def plot_parameter_comparison(self):
        """Plot comparison between true and fitted parameters"""
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
        
        # Transit parameters
        transit_params = ['Transit Depth', 'Transit Duration (days)']
        true_transit = [self.true_params['transit_depth'], self.true_params['transit_duration']]
        fitted_transit = [self.transit_best_params[0], self.transit_best_params[1]]
        
        x_pos = np.arange(len(transit_params))
        width = 0.35
        
        ax1.bar(x_pos - width/2, true_transit, width, label='True', alpha=0.7, color='green')
        ax1.bar(x_pos + width/2, fitted_transit, width, label='Fitted', alpha=0.7, color='red')
        ax1.set_xlabel('Parameters')
        ax1.set_ylabel('Values')
        ax1.set_title('Transit Parameters Comparison')
        ax1.set_xticks(x_pos)
        ax1.set_xticklabels(transit_params)
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # RV parameters
        rv_params = ['Period (days)', 'Amplitude (m/s)', 'Systemic Vel (m/s)', 'Phase (rad)']
        true_rv = [self.true_params['period'], self.true_params['rv_amplitude'],
                  self.true_params['systemic_velocity'], self.true_params['phase_offset']]
        fitted_rv = list(self.rv_best_params)
        
        x_pos = np.arange(len(rv_params))
        
        ax2.bar(x_pos - width/2, true_rv, width, label='True', alpha=0.7, color='green')
        ax2.bar(x_pos + width/2, fitted_rv, width, label='Fitted', alpha=0.7, color='red')
        ax2.set_xlabel('Parameters')
        ax2.set_ylabel('Values')
        ax2.set_title('RV Parameters Comparison')
        ax2.set_xticks(x_pos)
        ax2.set_xticklabels(rv_params, rotation=45)
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
    
    def monte_carlo_uncertainty(self, n_trials=1000):
        """Estimate parameter uncertainties using Monte Carlo method"""
        print("\n=== MONTE CARLO UNCERTAINTY ESTIMATION ===")
        
        # Store results
        transit_samples = []
        rv_samples = []
        
        for i in range(n_trials):
            # Add noise to data
            noisy_transit = self.transit_flux + np.random.normal(0, self.transit_error)
            noisy_rv = self.rv_velocity + np.random.normal(0, self.rv_error)
            
            # Temporarily replace data
            original_transit = self.transit_flux.copy()
            original_rv = self.rv_velocity.copy()
            
            self.transit_flux = noisy_transit
            self.rv_velocity = noisy_rv
            
            # Fit parameters
            try:
                # Transit fit
                result_transit = minimize(self.chi_squared_transit, self.transit_best_params, 
                                        bounds=[(0.001, 0.05), (0.05, 0.3)], method='L-BFGS-B')
                if result_transit.success:
                    transit_samples.append(result_transit.x)
                
                # RV fit
                result_rv = minimize(self.chi_squared_rv, self.rv_best_params, 
                                   bounds=[(2.0, 5.0), (10.0, 100.0), (-50.0, 50.0), (-np.pi, np.pi)], 
                                   method='L-BFGS-B')
                if result_rv.success:
                    rv_samples.append(result_rv.x)
                    
            except:
                continue
            
            # Restore original data
            self.transit_flux = original_transit
            self.rv_velocity = original_rv
        
        # Calculate uncertainties
        transit_samples = np.array(transit_samples)
        rv_samples = np.array(rv_samples)
        
        if len(transit_samples) > 0:
            transit_std = np.std(transit_samples, axis=0)
            print(f"Transit parameter uncertainties (1σ):")
            print(f"  Depth: {self.transit_best_params[0]:.6f} ± {transit_std[0]:.6f}")
            print(f"  Duration: {self.transit_best_params[1]:.6f} ± {transit_std[1]:.6f}")
        
        if len(rv_samples) > 0:
            rv_std = np.std(rv_samples, axis=0)
            print(f"RV parameter uncertainties (1σ):")
            print(f"  Period: {self.rv_best_params[0]:.6f} ± {rv_std[0]:.6f}")
            print(f"  Amplitude: {self.rv_best_params[1]:.6f} ± {rv_std[1]:.6f}")
            print(f"  Systemic velocity: {self.rv_best_params[2]:.6f} ± {rv_std[2]:.6f}")
            print(f"  Phase offset: {self.rv_best_params[3]:.6f} ± {rv_std[3]:.6f}")

# Run the complete analysis
print("🚀 Starting Exoplanet Parameter Optimization Analysis")
print("=" * 60)

# Initialize analysis
analysis = ExoplanetAnalysis()

# Generate synthetic data
analysis.generate_synthetic_data()

# Fit both datasets
analysis.fit_transit_data()
analysis.fit_rv_data()

# Calculate physical parameters
analysis.calculate_physical_parameters()

# Create plots
analysis.plot_results()

# Estimate uncertainties
analysis.monte_carlo_uncertainty(n_trials=500)

print("\n" + "=" * 60)
print("🎯 Analysis Complete!")
print("The optimization successfully recovered the planet parameters from noisy observational data.")