import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, quad
from scipy.optimize import fsolve
from mpl_toolkits.mplot3d import Axes3D

# Physical constants
g = 9.8  # gravity (m/s^2)

# Boundary conditions
x0, y0 = 0.0, 0.0  # Start point A
x1, y1 = 2.0, -1.0  # End point B

# 1. Brachistochrone (Cycloid) solution
def find_cycloid_parameter(x_target, y_target):
    """Find the parameter r for cycloid passing through (x_target, y_target)"""
    def equations(params):
        r, theta_end = params
        x = r * (theta_end - np.sin(theta_end))
        y = -r * (1 - np.cos(theta_end))
        return [x - x_target, y - y_target]
    
    # Initial guess
    r_guess = np.sqrt(x_target**2 + y_target**2) / 2
    theta_guess = np.pi
    solution = fsolve(equations, [r_guess, theta_guess])
    return solution[0], solution[1]

# Find cycloid parameters
r_cycloid, theta_max = find_cycloid_parameter(x1, y1)

# Generate cycloid curve
theta_array = np.linspace(0, theta_max, 500)
x_cycloid = r_cycloid * (theta_array - np.sin(theta_array))
y_cycloid = -r_cycloid * (1 - np.cos(theta_array))

# 2. Straight line path
x_line = np.linspace(x0, x1, 500)
y_line = y0 + (y1 - y0) / (x1 - x0) * (x_line - x0)

# 3. Parabolic path
x_parabola = np.linspace(x0, x1, 500)
# Parabola passing through (0,0) and (x1,y1) with vertex control
a_para = y1 / (x1**2)
y_parabola = a_para * x_parabola**2

# Function to calculate travel time along a curve
def calculate_travel_time(x_path, y_path):
    """Calculate time to travel along a path under gravity"""
    time = 0.0
    for i in range(1, len(x_path)):
        dx = x_path[i] - x_path[i-1]
        dy = y_path[i] - y_path[i-1]
        y_avg = (y_path[i] + y_path[i-1]) / 2
        
        if y_avg >= 0:  # No motion if above starting point
            continue
            
        ds = np.sqrt(dx**2 + dy**2)
        v_avg = np.sqrt(2 * g * abs(y_avg))
        
        if v_avg > 0:
            time += ds / v_avg
    
    return time

# Calculate times for each path
time_cycloid = calculate_travel_time(x_cycloid, y_cycloid)
time_line = calculate_travel_time(x_line, y_line)
time_parabola = calculate_travel_time(x_parabola, y_parabola)

# Calculate velocities along each path
def calculate_velocity_profile(x_path, y_path):
    """Calculate velocity at each point along the path"""
    velocities = np.zeros(len(y_path))
    for i in range(len(y_path)):
        if y_path[i] < 0:
            velocities[i] = np.sqrt(2 * g * abs(y_path[i]))
    return velocities

v_cycloid = calculate_velocity_profile(x_cycloid, y_cycloid)
v_line = calculate_velocity_profile(x_line, y_line)
v_parabola = calculate_velocity_profile(x_parabola, y_parabola)

# Visualization
fig = plt.figure(figsize=(18, 12))

# Plot 1: Path comparison
ax1 = plt.subplot(2, 3, 1)
ax1.plot(x_cycloid, y_cycloid, 'b-', linewidth=2.5, label=f'Brachistochrone (t={time_cycloid:.4f}s)')
ax1.plot(x_line, y_line, 'r--', linewidth=2, label=f'Straight line (t={time_line:.4f}s)')
ax1.plot(x_parabola, y_parabola, 'g:', linewidth=2, label=f'Parabola (t={time_parabola:.4f}s)')
ax1.plot([x0, x1], [y0, y1], 'ko', markersize=10)
ax1.text(x0-0.1, y0+0.1, 'A', fontsize=14, fontweight='bold')
ax1.text(x1+0.1, y1, 'B', fontsize=14, fontweight='bold')
ax1.set_xlabel('x (m)', fontsize=12)
ax1.set_ylabel('y (m)', fontsize=12)
ax1.set_title('Path Comparison', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.axis('equal')

# Plot 2: Velocity profiles
ax2 = plt.subplot(2, 3, 2)
ax2.plot(x_cycloid, v_cycloid, 'b-', linewidth=2.5, label='Brachistochrone')
ax2.plot(x_line, v_line, 'r--', linewidth=2, label='Straight line')
ax2.plot(x_parabola, v_parabola, 'g:', linewidth=2, label='Parabola')
ax2.set_xlabel('x (m)', fontsize=12)
ax2.set_ylabel('Velocity (m/s)', fontsize=12)
ax2.set_title('Velocity Along Path', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# Plot 3: Time comparison bar chart
ax3 = plt.subplot(2, 3, 3)
paths = ['Brachistochrone', 'Straight Line', 'Parabola']
times = [time_cycloid, time_line, time_parabola]
colors = ['blue', 'red', 'green']
bars = ax3.bar(paths, times, color=colors, alpha=0.7, edgecolor='black', linewidth=2)
ax3.set_ylabel('Travel Time (s)', fontsize=12)
ax3.set_title('Time Comparison', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3, axis='y')
for i, (bar, time) in enumerate(zip(bars, times)):
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height,
            f'{time:.4f}s\n({(time/time_cycloid-1)*100:+.1f}%)',
            ha='center', va='bottom', fontsize=10, fontweight='bold')

# Plot 4: Cycloid generation (rolling circle)
ax4 = plt.subplot(2, 3, 4)
n_circles = 6
theta_samples = np.linspace(0, theta_max, n_circles)
for theta in theta_samples:
    # Circle position
    cx = r_cycloid * theta
    cy = -r_cycloid
    circle = plt.Circle((cx, cy), r_cycloid, fill=False, color='gray', 
                        linestyle='--', alpha=0.5, linewidth=1)
    ax4.add_patch(circle)
    # Point on circle
    px = r_cycloid * (theta - np.sin(theta))
    py = -r_cycloid * (1 - np.cos(theta))
    ax4.plot(px, py, 'ro', markersize=6, alpha=0.7)

ax4.plot(x_cycloid, y_cycloid, 'b-', linewidth=2.5, label='Cycloid curve')
ax4.set_xlabel('x (m)', fontsize=12)
ax4.set_ylabel('y (m)', fontsize=12)
ax4.set_title('Cycloid Generation (Rolling Circle)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.axis('equal')

# Plot 5: Path length comparison
def calculate_path_length(x_path, y_path):
    length = 0.0
    for i in range(1, len(x_path)):
        dx = x_path[i] - x_path[i-1]
        dy = y_path[i] - y_path[i-1]
        length += np.sqrt(dx**2 + dy**2)
    return length

length_cycloid = calculate_path_length(x_cycloid, y_cycloid)
length_line = calculate_path_length(x_line, y_line)
length_parabola = calculate_path_length(x_parabola, y_parabola)

ax5 = plt.subplot(2, 3, 5)
lengths = [length_cycloid, length_line, length_parabola]
bars2 = ax5.bar(paths, lengths, color=colors, alpha=0.7, edgecolor='black', linewidth=2)
ax5.set_ylabel('Path Length (m)', fontsize=12)
ax5.set_title('Path Length Comparison', fontsize=14, fontweight='bold')
ax5.grid(True, alpha=0.3, axis='y')
for bar, length in zip(bars2, lengths):
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height,
            f'{length:.4f}m',
            ha='center', va='bottom', fontsize=10, fontweight='bold')

# Plot 6: 3D trajectory with time dimension
ax6 = fig.add_subplot(2, 3, 6, projection='3d')

# Calculate time arrays for each path
def calculate_time_array(x_path, y_path):
    time_array = np.zeros(len(x_path))
    for i in range(1, len(x_path)):
        dx = x_path[i] - x_path[i-1]
        dy = y_path[i] - y_path[i-1]
        y_avg = (y_path[i] + y_path[i-1]) / 2
        
        if y_avg >= 0:
            continue
            
        ds = np.sqrt(dx**2 + dy**2)
        v_avg = np.sqrt(2 * g * abs(y_avg))
        
        if v_avg > 0:
            time_array[i] = time_array[i-1] + ds / v_avg
    return time_array

t_cycloid_array = calculate_time_array(x_cycloid, y_cycloid)
t_line_array = calculate_time_array(x_line, y_line)
t_parabola_array = calculate_time_array(x_parabola, y_parabola)

ax6.plot(x_cycloid, y_cycloid, t_cycloid_array, 'b-', linewidth=2.5, label='Brachistochrone')
ax6.plot(x_line, y_line, t_line_array, 'r--', linewidth=2, label='Straight line')
ax6.plot(x_parabola, y_parabola, t_parabola_array, 'g:', linewidth=2, label='Parabola')
ax6.set_xlabel('x (m)', fontsize=10)
ax6.set_ylabel('y (m)', fontsize=10)
ax6.set_zlabel('Time (s)', fontsize=10)
ax6.set_title('3D: Position vs Time', fontsize=14, fontweight='bold')
ax6.legend(fontsize=9)
ax6.view_init(elev=20, azim=45)

plt.tight_layout()
plt.show()

# Print detailed results
print("="*70)
print("BRACHISTOCHRONE PROBLEM SOLUTION")
print("="*70)
print(f"\nBoundary Conditions:")
print(f"  Start point A: ({x0}, {y0})")
print(f"  End point B: ({x1}, {y1})")
print(f"  Gravity: g = {g} m/s²")
print(f"\nCycloid Parameters:")
print(f"  Radius r = {r_cycloid:.6f} m")
print(f"  Maximum angle θ_max = {theta_max:.6f} rad = {np.degrees(theta_max):.2f}°")
print(f"\n{'Path Type':<20} {'Time (s)':<15} {'Length (m)':<15} {'Time Diff (%)':<15}")
print("-"*70)
print(f"{'Brachistochrone':<20} {time_cycloid:<15.6f} {length_cycloid:<15.6f} {'0.00% (optimal)':<15}")
print(f"{'Straight Line':<20} {time_line:<15.6f} {length_line:<15.6f} {f'+{(time_line/time_cycloid-1)*100:.2f}%':<15}")
print(f"{'Parabola':<20} {time_parabola:<15.6f} {length_parabola:<15.6f} {f'+{(time_parabola/time_cycloid-1)*100:.2f}%':<15}")
print("="*70)
print(f"\nKey Insight: The brachistochrone is {(time_line/time_cycloid-1)*100:.2f}% faster than")
print(f"the straight line, even though it's {(length_cycloid/length_line-1)*100:.2f}% longer!")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import differential_evolution
import warnings
warnings.filterwarnings('ignore')

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

class MarsExplorationMission:
    """
    Mars Exploration Mission Design Problem
    Decision variables:
    - x[0]: Number of cameras (0-5)
    - x[1]: Number of spectrometers (0-3)
    - x[2]: Number of drilling instruments (0-2)
    - x[3]: Mission duration in days (180-720)
    - x[4]: Solar panel area in m^2 (10-50)
    - x[5]: Communication bandwidth in Mbps (1-20)
    """
    
    def __init__(self):
        # Instrument specifications
        self.camera_cost = 15  # Million USD
        self.spectro_cost = 25  # Million USD
        self.drill_cost = 40   # Million USD
        
        self.camera_mass = 5   # kg
        self.spectro_mass = 12 # kg
        self.drill_mass = 20   # kg
        
        self.camera_power = 50  # Watts
        self.spectro_power = 80 # Watts
        self.drill_power = 150  # Watts
        
        # Base costs
        self.base_spacecraft_cost = 200  # Million USD
        self.daily_operation_cost = 0.5  # Million USD per day
        
    def calculate_cost(self, x):
        """Calculate total mission cost in Million USD"""
        n_cameras = int(x[0])
        n_spectros = int(x[1])
        n_drills = int(x[2])
        duration = x[3]
        solar_area = x[4]
        bandwidth = x[5]
        
        # Instrument costs
        instrument_cost = (n_cameras * self.camera_cost + 
                          n_spectros * self.spectro_cost + 
                          n_drills * self.drill_cost)
        
        # Spacecraft cost (increases with mass and power requirements)
        total_mass = (n_cameras * self.camera_mass + 
                     n_spectros * self.spectro_mass + 
                     n_drills * self.drill_mass)
        
        spacecraft_cost = self.base_spacecraft_cost + total_mass * 2
        
        # Power system cost
        power_cost = solar_area * 3  # $3M per m^2
        
        # Communication system cost
        comm_cost = bandwidth * 5  # $5M per Mbps
        
        # Operations cost
        ops_cost = duration * self.daily_operation_cost
        
        total_cost = (instrument_cost + spacecraft_cost + 
                     power_cost + comm_cost + ops_cost)
        
        return total_cost
    
    def calculate_risk(self, x):
        """Calculate mission risk score (0-100, higher is worse)"""
        n_cameras = int(x[0])
        n_spectros = int(x[1])
        n_drills = int(x[2])
        duration = x[3]
        solar_area = x[4]
        bandwidth = x[5]
        
        # Complexity risk (more instruments = more risk)
        n_instruments = n_cameras + n_spectros + n_drills
        complexity_risk = n_instruments * 5
        
        # Duration risk (longer mission = more exposure to failures)
        duration_risk = (duration / 720) * 30
        
        # Power margin risk
        required_power = (n_cameras * self.camera_power + 
                         n_spectros * self.spectro_power + 
                         n_drills * self.drill_power)
        available_power = solar_area * 100  # 100W per m^2
        power_margin = available_power - required_power
        
        if power_margin < 0:
            power_risk = 50  # Severe risk if insufficient power
        elif power_margin < required_power * 0.2:
            power_risk = 30  # High risk if < 20% margin
        elif power_margin < required_power * 0.5:
            power_risk = 15  # Medium risk if < 50% margin
        else:
            power_risk = 5   # Low risk with good margin
        
        # Communication risk
        data_rate = n_instruments * 10  # MB per day per instrument
        comm_capability = bandwidth * 60  # MB per day
        
        if comm_capability < data_rate:
            comm_risk = 40
        elif comm_capability < data_rate * 1.5:
            comm_risk = 20
        else:
            comm_risk = 5
        
        total_risk = complexity_risk + duration_risk + power_risk + comm_risk
        
        return min(total_risk, 100)
    
    def calculate_science_return(self, x):
        """Calculate scientific return score (0-100, higher is better)"""
        n_cameras = int(x[0])
        n_spectros = int(x[1])
        n_drills = int(x[2])
        duration = x[3]
        solar_area = x[4]
        bandwidth = x[5]
        
        # Base science from instruments
        camera_science = n_cameras * 10
        spectro_science = n_spectros * 20
        drill_science = n_drills * 30
        
        # Synergy bonus (instruments work better together)
        if n_cameras > 0 and n_spectros > 0:
            camera_science *= 1.3
            spectro_science *= 1.2
        
        if n_spectros > 0 and n_drills > 0:
            spectro_science *= 1.3
            drill_science *= 1.2
        
        # Duration factor (longer mission = more data, but diminishing returns)
        duration_factor = np.log(duration / 180 + 1) / np.log(5)
        
        # Communication limitation
        data_generated = (n_cameras + n_spectros + n_drills) * duration * 10
        data_transmitted = bandwidth * duration * 60
        transmission_efficiency = min(data_transmitted / (data_generated + 1), 1.0)
        
        total_science = (camera_science + spectro_science + drill_science) * \
                       duration_factor * transmission_efficiency
        
        return min(total_science, 100)
    
    def evaluate(self, x):
        """Evaluate all three objectives"""
        cost = self.calculate_cost(x)
        risk = self.calculate_risk(x)
        science = self.calculate_science_return(x)
        
        # Return as minimization problem (negate science)
        return cost, risk, -science


class NSGA2Optimizer:
    """NSGA-II Multi-objective Genetic Algorithm"""
    
    def __init__(self, problem, pop_size=100, n_gen=100):
        self.problem = problem
        self.pop_size = pop_size
        self.n_gen = n_gen
        
        # Bounds for decision variables
        self.bounds = [
            (0, 5),      # cameras
            (0, 3),      # spectrometers
            (0, 2),      # drills
            (180, 720),  # duration
            (10, 50),    # solar area
            (1, 20)      # bandwidth
        ]
    
    def initialize_population(self):
        """Create initial random population"""
        pop = np.random.rand(self.pop_size, 6)
        for i in range(6):
            pop[:, i] = pop[:, i] * (self.bounds[i][1] - self.bounds[i][0]) + self.bounds[i][0]
        return pop
    
    def fast_non_dominated_sort(self, objectives):
        """Fast non-dominated sorting"""
        n = len(objectives)
        domination_count = np.zeros(n)
        dominated_solutions = [[] for _ in range(n)]
        fronts = [[]]
        
        for i in range(n):
            for j in range(n):
                if i == j:
                    continue
                if self.dominates(objectives[i], objectives[j]):
                    dominated_solutions[i].append(j)
                elif self.dominates(objectives[j], objectives[i]):
                    domination_count[i] += 1
            
            if domination_count[i] == 0:
                fronts[0].append(i)
        
        k = 0
        while len(fronts[k]) > 0:
            next_front = []
            for i in fronts[k]:
                for j in dominated_solutions[i]:
                    domination_count[j] -= 1
                    if domination_count[j] == 0:
                        next_front.append(j)
            k += 1
            fronts.append(next_front)
        
        return fronts[:-1]
    
    def dominates(self, obj1, obj2):
        """Check if obj1 dominates obj2 (minimization)"""
        return all(obj1 <= obj2) and any(obj1 < obj2)
    
    def crowding_distance(self, objectives, front):
        """Calculate crowding distance for solutions in a front"""
        n = len(front)
        if n <= 2:
            return np.full(n, np.inf)
        
        distances = np.zeros(n)
        
        for m in range(3):  # 3 objectives
            sorted_indices = np.argsort([objectives[i][m] for i in front])
            distances[sorted_indices[0]] = np.inf
            distances[sorted_indices[-1]] = np.inf
            
            obj_range = objectives[front[sorted_indices[-1]]][m] - \
                       objectives[front[sorted_indices[0]]][m]
            
            if obj_range == 0:
                continue
            
            for i in range(1, n - 1):
                distances[sorted_indices[i]] += \
                    (objectives[front[sorted_indices[i + 1]]][m] - \
                     objectives[front[sorted_indices[i - 1]]][m]) / obj_range
        
        return distances
    
    def tournament_selection(self, population, objectives, fronts):
        """Binary tournament selection"""
        idx1, idx2 = np.random.choice(len(population), 2, replace=False)
        
        # Find fronts for both
        front1, front2 = None, None
        for i, front in enumerate(fronts):
            if idx1 in front:
                front1 = i
            if idx2 in front:
                front2 = i
        
        if front1 < front2:
            return population[idx1].copy()
        elif front2 < front1:
            return population[idx2].copy()
        else:
            # Same front, use crowding distance
            front = fronts[front1]
            distances = self.crowding_distance(objectives, front)
            pos1 = front.index(idx1)
            pos2 = front.index(idx2)
            
            if distances[pos1] > distances[pos2]:
                return population[idx1].copy()
            else:
                return population[idx2].copy()
    
    def crossover(self, parent1, parent2):
        """Simulated Binary Crossover (SBX)"""
        eta = 20
        child1 = parent1.copy()
        child2 = parent2.copy()
        
        for i in range(len(parent1)):
            if np.random.rand() < 0.5:
                if abs(parent1[i] - parent2[i]) > 1e-6:
                    u = np.random.rand()
                    if u <= 0.5:
                        beta = (2 * u) ** (1 / (eta + 1))
                    else:
                        beta = (1 / (2 * (1 - u))) ** (1 / (eta + 1))
                    
                    child1[i] = 0.5 * ((1 + beta) * parent1[i] + (1 - beta) * parent2[i])
                    child2[i] = 0.5 * ((1 - beta) * parent1[i] + (1 + beta) * parent2[i])
        
        return child1, child2
    
    def mutate(self, individual):
        """Polynomial mutation"""
        eta = 20
        mutated = individual.copy()
        
        for i in range(len(individual)):
            if np.random.rand() < 1.0 / len(individual):
                u = np.random.rand()
                lb, ub = self.bounds[i]
                
                delta1 = (individual[i] - lb) / (ub - lb)
                delta2 = (ub - individual[i]) / (ub - lb)
                
                if u < 0.5:
                    deltaq = (2 * u + (1 - 2 * u) * (1 - delta1) ** (eta + 1)) ** (1 / (eta + 1)) - 1
                else:
                    deltaq = 1 - (2 * (1 - u) + 2 * (u - 0.5) * (1 - delta2) ** (eta + 1)) ** (1 / (eta + 1))
                
                mutated[i] = individual[i] + deltaq * (ub - lb)
                mutated[i] = np.clip(mutated[i], lb, ub)
        
        return mutated
    
    def optimize(self):
        """Run NSGA-II optimization"""
        population = self.initialize_population()
        
        for gen in range(self.n_gen):
            # Evaluate population
            objectives = np.array([self.problem.evaluate(ind) for ind in population])
            
            # Non-dominated sorting
            fronts = self.fast_non_dominated_sort(objectives)
            
            # Generate offspring
            offspring = []
            while len(offspring) < self.pop_size:
                parent1 = self.tournament_selection(population, objectives, fronts)
                parent2 = self.tournament_selection(population, objectives, fronts)
                
                child1, child2 = self.crossover(parent1, parent2)
                child1 = self.mutate(child1)
                child2 = self.mutate(child2)
                
                offspring.append(child1)
                if len(offspring) < self.pop_size:
                    offspring.append(child2)
            
            offspring = np.array(offspring)
            
            # Combine parent and offspring populations
            combined = np.vstack([population, offspring])
            combined_obj = np.array([self.problem.evaluate(ind) for ind in combined])
            
            # Select next generation
            fronts = self.fast_non_dominated_sort(combined_obj)
            new_population = []
            
            for front in fronts:
                if len(new_population) + len(front) <= self.pop_size:
                    new_population.extend([combined[i] for i in front])
                else:
                    # Sort by crowding distance
                    distances = self.crowding_distance(combined_obj, front)
                    sorted_indices = np.argsort(distances)[::-1]
                    remaining = self.pop_size - len(new_population)
                    new_population.extend([combined[front[i]] for i in sorted_indices[:remaining]])
                    break
            
            population = np.array(new_population)
            
            if (gen + 1) % 20 == 0:
                print(f"Generation {gen + 1}/{self.n_gen} completed")
        
        # Final evaluation
        objectives = np.array([self.problem.evaluate(ind) for ind in population])
        fronts = self.fast_non_dominated_sort(objectives)
        pareto_front_indices = fronts[0]
        
        pareto_solutions = population[pareto_front_indices]
        pareto_objectives = objectives[pareto_front_indices]
        
        return pareto_solutions, pareto_objectives


# Initialize problem
print("=" * 60)
print("Mars Exploration Mission Multi-Objective Optimization")
print("=" * 60)
print("\nObjectives:")
print("  1. Minimize Cost (Million USD)")
print("  2. Minimize Risk (0-100)")
print("  3. Maximize Scientific Return (0-100)")
print("\nRunning NSGA-II optimization...")
print("=" * 60)

mission = MarsExplorationMission()
optimizer = NSGA2Optimizer(mission, pop_size=100, n_gen=100)

# Run optimization
pareto_solutions, pareto_objectives = optimizer.optimize()

print("\n" + "=" * 60)
print(f"Optimization complete! Found {len(pareto_solutions)} Pareto-optimal solutions")
print("=" * 60)

# Convert objectives for display (negate science back to positive)
pareto_objectives_display = pareto_objectives.copy()
pareto_objectives_display[:, 2] = -pareto_objectives_display[:, 2]

# Create visualizations
fig = plt.figure(figsize=(20, 12))

# 3D Pareto Front
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
scatter = ax1.scatter(pareto_objectives_display[:, 0], 
                     pareto_objectives_display[:, 1], 
                     pareto_objectives_display[:, 2],
                     c=pareto_objectives_display[:, 2], 
                     cmap='viridis', 
                     s=100, 
                     alpha=0.6,
                     edgecolors='black')
ax1.set_xlabel('Cost (Million USD)', fontsize=10, fontweight='bold')
ax1.set_ylabel('Risk Score', fontsize=10, fontweight='bold')
ax1.set_zlabel('Scientific Return', fontsize=10, fontweight='bold')
ax1.set_title('3D Pareto Front', fontsize=12, fontweight='bold')
plt.colorbar(scatter, ax=ax1, label='Scientific Return', shrink=0.5)

# Cost vs Risk
ax2 = fig.add_subplot(2, 3, 2)
scatter2 = ax2.scatter(pareto_objectives_display[:, 0], 
                       pareto_objectives_display[:, 1],
                       c=pareto_objectives_display[:, 2], 
                       cmap='viridis', 
                       s=100, 
                       alpha=0.6,
                       edgecolors='black')
ax2.set_xlabel('Cost (Million USD)', fontsize=10, fontweight='bold')
ax2.set_ylabel('Risk Score', fontsize=10, fontweight='bold')
ax2.set_title('Cost vs Risk Trade-off', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
plt.colorbar(scatter2, ax=ax2, label='Scientific Return')

# Cost vs Science
ax3 = fig.add_subplot(2, 3, 3)
scatter3 = ax3.scatter(pareto_objectives_display[:, 0], 
                       pareto_objectives_display[:, 2],
                       c=pareto_objectives_display[:, 1], 
                       cmap='coolwarm_r', 
                       s=100, 
                       alpha=0.6,
                       edgecolors='black')
ax3.set_xlabel('Cost (Million USD)', fontsize=10, fontweight='bold')
ax3.set_ylabel('Scientific Return', fontsize=10, fontweight='bold')
ax3.set_title('Cost vs Scientific Return', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
plt.colorbar(scatter3, ax=ax3, label='Risk Score')

# Risk vs Science
ax4 = fig.add_subplot(2, 3, 4)
scatter4 = ax4.scatter(pareto_objectives_display[:, 1], 
                       pareto_objectives_display[:, 2],
                       c=pareto_objectives_display[:, 0], 
                       cmap='plasma', 
                       s=100, 
                       alpha=0.6,
                       edgecolors='black')
ax4.set_xlabel('Risk Score', fontsize=10, fontweight='bold')
ax4.set_ylabel('Scientific Return', fontsize=10, fontweight='bold')
ax4.set_title('Risk vs Scientific Return', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
plt.colorbar(scatter4, ax=ax4, label='Cost (M USD)')

# Decision variable distributions
ax5 = fig.add_subplot(2, 3, 5)
variables = ['Cameras', 'Spectros', 'Drills', 'Duration', 'Solar', 'Bandwidth']
for i in range(3):  # Plot only instrument counts
    ax5.hist(pareto_solutions[:, i], bins=10, alpha=0.5, label=variables[i])
ax5.set_xlabel('Count', fontsize=10, fontweight='bold')
ax5.set_ylabel('Frequency', fontsize=10, fontweight='bold')
ax5.set_title('Instrument Distribution in Pareto Solutions', fontsize=12, fontweight='bold')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Mission configuration scatter
ax6 = fig.add_subplot(2, 3, 6)
total_instruments = (pareto_solutions[:, 0] + 
                    pareto_solutions[:, 1] + 
                    pareto_solutions[:, 2])
scatter6 = ax6.scatter(pareto_solutions[:, 3], 
                       total_instruments,
                       c=pareto_objectives_display[:, 2], 
                       cmap='viridis', 
                       s=100, 
                       alpha=0.6,
                       edgecolors='black')
ax6.set_xlabel('Mission Duration (days)', fontsize=10, fontweight='bold')
ax6.set_ylabel('Total Instruments', fontsize=10, fontweight='bold')
ax6.set_title('Mission Duration vs Instrument Count', fontsize=12, fontweight='bold')
ax6.grid(True, alpha=0.3)
plt.colorbar(scatter6, ax=ax6, label='Scientific Return')

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

# Statistical analysis
print("\n" + "=" * 60)
print("Statistical Analysis of Pareto-Optimal Solutions")
print("=" * 60)
print(f"\nCost (Million USD):")
print(f"  Min: ${pareto_objectives_display[:, 0].min():.2f}M")
print(f"  Max: ${pareto_objectives_display[:, 0].max():.2f}M")
print(f"  Mean: ${pareto_objectives_display[:, 0].mean():.2f}M")
print(f"  Std: ${pareto_objectives_display[:, 0].std():.2f}M")

print(f"\nRisk Score:")
print(f"  Min: {pareto_objectives_display[:, 1].min():.2f}")
print(f"  Max: {pareto_objectives_display[:, 1].max():.2f}")
print(f"  Mean: {pareto_objectives_display[:, 1].mean():.2f}")
print(f"  Std: {pareto_objectives_display[:, 1].std():.2f}")

print(f"\nScientific Return:")
print(f"  Min: {pareto_objectives_display[:, 2].min():.2f}")
print(f"  Max: {pareto_objectives_display[:, 2].max():.2f}")
print(f"  Mean: {pareto_objectives_display[:, 2].mean():.2f}")
print(f"  Std: {pareto_objectives_display[:, 2].std():.2f}")

# Find interesting solutions
print("\n" + "=" * 60)
print("Representative Mission Configurations")
print("=" * 60)

# Lowest cost
idx_low_cost = np.argmin(pareto_objectives_display[:, 0])
print("\n1. LOWEST COST MISSION:")
print(f"   Cost: ${pareto_objectives_display[idx_low_cost, 0]:.2f}M")
print(f"   Risk: {pareto_objectives_display[idx_low_cost, 1]:.2f}")
print(f"   Science: {pareto_objectives_display[idx_low_cost, 2]:.2f}")
print(f"   Configuration:")
print(f"     - Cameras: {int(pareto_solutions[idx_low_cost, 0])}")
print(f"     - Spectrometers: {int(pareto_solutions[idx_low_cost, 1])}")
print(f"     - Drills: {int(pareto_solutions[idx_low_cost, 2])}")
print(f"     - Duration: {int(pareto_solutions[idx_low_cost, 3])} days")
print(f"     - Solar panels: {pareto_solutions[idx_low_cost, 4]:.1f} m²")
print(f"     - Bandwidth: {pareto_solutions[idx_low_cost, 5]:.1f} Mbps")

# Highest science
idx_high_science = np.argmax(pareto_objectives_display[:, 2])
print("\n2. HIGHEST SCIENCE MISSION:")
print(f"   Cost: ${pareto_objectives_display[idx_high_science, 0]:.2f}M")
print(f"   Risk: {pareto_objectives_display[idx_high_science, 1]:.2f}")
print(f"   Science: {pareto_objectives_display[idx_high_science, 2]:.2f}")
print(f"   Configuration:")
print(f"     - Cameras: {int(pareto_solutions[idx_high_science, 0])}")
print(f"     - Spectrometers: {int(pareto_solutions[idx_high_science, 1])}")
print(f"     - Drills: {int(pareto_solutions[idx_high_science, 2])}")
print(f"     - Duration: {int(pareto_solutions[idx_high_science, 3])} days")
print(f"     - Solar panels: {pareto_solutions[idx_high_science, 4]:.1f} m²")
print(f"     - Bandwidth: {pareto_solutions[idx_high_science, 5]:.1f} Mbps")

# Lowest risk
idx_low_risk = np.argmin(pareto_objectives_display[:, 1])
print("\n3. LOWEST RISK MISSION:")
print(f"   Cost: ${pareto_objectives_display[idx_low_risk, 0]:.2f}M")
print(f"   Risk: {pareto_objectives_display[idx_low_risk, 1]:.2f}")
print(f"   Science: {pareto_objectives_display[idx_low_risk, 2]:.2f}")
print(f"   Configuration:")
print(f"     - Cameras: {int(pareto_solutions[idx_low_risk, 0])}")
print(f"     - Spectrometers: {int(pareto_solutions[idx_low_risk, 1])}")
print(f"     - Drills: {int(pareto_solutions[idx_low_risk, 2])}")
print(f"     - Duration: {int(pareto_solutions[idx_low_risk, 3])} days")
print(f"     - Solar panels: {pareto_solutions[idx_low_risk, 4]:.1f} m²")
print(f"     - Bandwidth: {pareto_solutions[idx_low_risk, 5]:.1f} Mbps")

# Balanced solution (closest to ideal point)
normalized_obj = (pareto_objectives_display - pareto_objectives_display.min(axis=0)) / \
                 (pareto_objectives_display.max(axis=0) - pareto_objectives_display.min(axis=0))
normalized_obj[:, 2] = 1 - normalized_obj[:, 2]  # For science (higher is better)
distances = np.sqrt(np.sum(normalized_obj**2, axis=1))
idx_balanced = np.argmin(distances)

print("\n4. BALANCED MISSION (Best compromise):")
print(f"   Cost: ${pareto_objectives_display[idx_balanced, 0]:.2f}M")
print(f"   Risk: {pareto_objectives_display[idx_balanced, 1]:.2f}")
print(f"   Science: {pareto_objectives_display[idx_balanced, 2]:.2f}")
print(f"   Configuration:")
print(f"     - Cameras: {int(pareto_solutions[idx_balanced, 0])}")
print(f"     - Spectrometers: {int(pareto_solutions[idx_balanced, 1])}")
print(f"     - Drills: {int(pareto_solutions[idx_balanced, 2])}")
print(f"     - Duration: {int(pareto_solutions[idx_balanced, 3])} days")
print(f"     - Solar panels: {pareto_solutions[idx_balanced, 4]:.1f} m²")
print(f"     - Bandwidth: {pareto_solutions[idx_balanced, 5]:.1f} Mbps")

print("\n" + "=" * 60)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import beta, binom
from scipy.special import betaln
import seaborn as sns

# Set style for better visualizations
sns.set_style("whitegrid")
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 11

class BayesianLifeProbability:
    """
    Bayesian estimation of life existence probability on exoplanets
    """
    
    def __init__(self, prior_alpha, prior_beta):
        """
        Initialize with prior Beta distribution parameters
        
        Parameters:
        -----------
        prior_alpha : float
            Alpha parameter of prior Beta distribution (successes + 1)
        prior_beta : float
            Beta parameter of prior Beta distribution (failures + 1)
        """
        self.prior_alpha = prior_alpha
        self.prior_beta = prior_beta
        
    def prior_distribution(self, theta):
        """Calculate prior probability density"""
        return beta.pdf(theta, self.prior_alpha, self.prior_beta)
    
    def likelihood(self, theta, n_obs, n_detections):
        """
        Calculate likelihood of observing n_detections in n_obs observations
        
        Parameters:
        -----------
        theta : array-like
            Probability values
        n_obs : int
            Number of observations
        n_detections : int
            Number of positive biosignature detections
        """
        return binom.pmf(n_detections, n_obs, theta)
    
    def posterior_distribution(self, theta, n_obs, n_detections):
        """
        Calculate posterior probability density after observing data
        Using conjugate prior property: Beta + Binomial = Beta
        """
        post_alpha = self.prior_alpha + n_detections
        post_beta = self.prior_beta + (n_obs - n_detections)
        return beta.pdf(theta, post_alpha, post_beta)
    
    def update_posterior_params(self, n_obs, n_detections):
        """Get updated posterior parameters"""
        post_alpha = self.prior_alpha + n_detections
        post_beta = self.prior_beta + (n_obs - n_detections)
        return post_alpha, post_beta
    
    def posterior_statistics(self, n_obs, n_detections):
        """Calculate posterior mean, mode, variance, and credible interval"""
        post_alpha, post_beta = self.update_posterior_params(n_obs, n_detections)
        
        # Mean
        mean = post_alpha / (post_alpha + post_beta)
        
        # Mode (for alpha, beta > 1)
        if post_alpha > 1 and post_beta > 1:
            mode = (post_alpha - 1) / (post_alpha + post_beta - 2)
        else:
            mode = None
        
        # Variance
        variance = (post_alpha * post_beta) / \
                   ((post_alpha + post_beta)**2 * (post_alpha + post_beta + 1))
        std = np.sqrt(variance)
        
        # 95% Credible Interval
        ci_lower = beta.ppf(0.025, post_alpha, post_beta)
        ci_upper = beta.ppf(0.975, post_alpha, post_beta)
        
        return {
            'mean': mean,
            'mode': mode,
            'std': std,
            'variance': variance,
            'ci_95': (ci_lower, ci_upper)
        }
    
    def information_gain(self, n_obs, n_detections):
        """
        Calculate Kullback-Leibler divergence from prior to posterior
        Measures information gain from observations
        """
        post_alpha, post_beta = self.update_posterior_params(n_obs, n_detections)
        
        # KL divergence for Beta distributions
        kl = betaln(self.prior_alpha, self.prior_beta) - \
             betaln(post_alpha, post_beta) + \
             (post_alpha - self.prior_alpha) * \
             (np.euler_gamma + np.log(post_alpha + post_beta - 2) - 
              np.log(post_alpha - 1) if post_alpha > 1 else 0) + \
             (post_beta - self.prior_beta) * \
             (np.euler_gamma + np.log(post_alpha + post_beta - 2) - 
              np.log(post_beta - 1) if post_beta > 1 else 0)
        
        return abs(kl)

# Example: Exoplanet biosignature detection
# Prior: Moderately skeptical (alpha=2, beta=10)
# This represents prior belief that ~16.7% of habitable zone planets have life

bayesian_model = BayesianLifeProbability(prior_alpha=2, prior_beta=10)

# Observational data: 15 observations, 5 positive biosignature detections
n_observations = 15
n_detections = 5

# Generate probability values for plotting
theta_values = np.linspace(0, 1, 1000)

# Calculate distributions
prior = bayesian_model.prior_distribution(theta_values)
likelihood = bayesian_model.likelihood(theta_values, n_observations, n_detections)
posterior = bayesian_model.posterior_distribution(theta_values, n_observations, n_detections)

# Normalize likelihood for visualization
likelihood_normalized = likelihood / np.trapz(likelihood, theta_values)

# Get posterior statistics
stats = bayesian_model.posterior_statistics(n_observations, n_detections)

# Calculate information gain
info_gain = bayesian_model.information_gain(n_observations, n_detections)

# Print results
print("="*70)
print("BAYESIAN ESTIMATION OF LIFE EXISTENCE PROBABILITY")
print("="*70)
print(f"\nPrior Distribution: Beta(α={bayesian_model.prior_alpha}, β={bayesian_model.prior_beta})")
print(f"Prior Mean: {bayesian_model.prior_alpha/(bayesian_model.prior_alpha + bayesian_model.prior_beta):.4f}")
print(f"\nObservational Data:")
print(f"  - Number of observations: {n_observations}")
print(f"  - Positive detections: {n_detections}")
print(f"  - Detection rate: {n_detections/n_observations:.2%}")
print(f"\nPosterior Statistics:")
print(f"  - Mean: {stats['mean']:.4f}")
if stats['mode']:
    print(f"  - Mode: {stats['mode']:.4f}")
print(f"  - Standard Deviation: {stats['std']:.4f}")
print(f"  - 95% Credible Interval: [{stats['ci_95'][0]:.4f}, {stats['ci_95'][1]:.4f}]")
print(f"\nInformation Gain (KL Divergence): {info_gain:.4f}")
print(f"Uncertainty Reduction: {(1 - stats['std']/np.sqrt(bayesian_model.prior_alpha * bayesian_model.prior_beta / 
      ((bayesian_model.prior_alpha + bayesian_model.prior_beta)**2 * 
       (bayesian_model.prior_alpha + bayesian_model.prior_beta + 1))))*100:.2f}%")
print("="*70)

# Visualization 1: Prior, Likelihood, and Posterior
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# Plot 1: Prior Distribution
axes[0, 0].fill_between(theta_values, prior, alpha=0.6, color='blue', label='Prior')
axes[0, 0].axvline(bayesian_model.prior_alpha/(bayesian_model.prior_alpha + bayesian_model.prior_beta), 
                    color='blue', linestyle='--', linewidth=2, label='Prior Mean')
axes[0, 0].set_xlabel(r'Probability of Life ($\theta$)', fontsize=12)
axes[0, 0].set_ylabel('Density', fontsize=12)
axes[0, 0].set_title('Prior Distribution: Beta(2, 10)', fontsize=14, fontweight='bold')
axes[0, 0].legend(fontsize=11)
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Likelihood Function
axes[0, 1].fill_between(theta_values, likelihood_normalized, alpha=0.6, color='green', label='Likelihood')
axes[0, 1].axvline(n_detections/n_observations, color='green', linestyle='--', 
                    linewidth=2, label='MLE')
axes[0, 1].set_xlabel(r'Probability of Life ($\theta$)', fontsize=12)
axes[0, 1].set_ylabel('Normalized Density', fontsize=12)
axes[0, 1].set_title(f'Likelihood: {n_detections}/{n_observations} Detections', fontsize=14, fontweight='bold')
axes[0, 1].legend(fontsize=11)
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Posterior Distribution
axes[1, 0].fill_between(theta_values, posterior, alpha=0.6, color='red', label='Posterior')
axes[1, 0].axvline(stats['mean'], color='red', linestyle='--', linewidth=2, label='Posterior Mean')
axes[1, 0].axvspan(stats['ci_95'][0], stats['ci_95'][1], alpha=0.2, color='red', 
                    label='95% Credible Interval')
axes[1, 0].set_xlabel(r'Probability of Life ($\theta$)', fontsize=12)
axes[1, 0].set_ylabel('Density', fontsize=12)
axes[1, 0].set_title('Posterior Distribution', fontsize=14, fontweight='bold')
axes[1, 0].legend(fontsize=11)
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: All distributions together
axes[1, 1].plot(theta_values, prior, 'b-', linewidth=2.5, label='Prior', alpha=0.8)
axes[1, 1].plot(theta_values, likelihood_normalized, 'g-', linewidth=2.5, label='Likelihood', alpha=0.8)
axes[1, 1].plot(theta_values, posterior, 'r-', linewidth=2.5, label='Posterior', alpha=0.8)
axes[1, 1].set_xlabel(r'Probability of Life ($\theta$)', fontsize=12)
axes[1, 1].set_ylabel('Density', fontsize=12)
axes[1, 1].set_title('Bayesian Update: Prior → Posterior', fontsize=14, fontweight='bold')
axes[1, 1].legend(fontsize=11)
axes[1, 1].grid(True, alpha=0.3)

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

# Visualization 2: Sequential Bayesian Update
fig, ax = plt.subplots(figsize=(14, 8))

# Simulate sequential observations
detections_sequence = [1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(detections_sequence) + 1))

# Start with prior
current_alpha = bayesian_model.prior_alpha
current_beta = bayesian_model.prior_beta

ax.plot(theta_values, beta.pdf(theta_values, current_alpha, current_beta), 
        linewidth=2.5, label=f'Initial Prior (n=0)', color=colors[0])

# Sequential update
cumulative_detections = 0
for i, detection in enumerate(detections_sequence):
    cumulative_detections += detection
    current_alpha += detection
    current_beta += (1 - detection)
    
    if (i + 1) % 3 == 0 or i == len(detections_sequence) - 1:
        ax.plot(theta_values, beta.pdf(theta_values, current_alpha, current_beta),
                linewidth=2, label=f'After n={i+1} obs ({cumulative_detections} det)', 
                color=colors[i+1], alpha=0.8)

ax.set_xlabel(r'Probability of Life ($\theta$)', fontsize=13)
ax.set_ylabel('Probability Density', fontsize=13)
ax.set_title('Sequential Bayesian Update: How Posterior Evolves with Data', 
             fontsize=15, fontweight='bold')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)

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

# Visualization 3: 3D Surface Plot - Posterior as Function of Observations
fig = plt.figure(figsize=(16, 12))

# Create 3D subplots
ax1 = fig.add_subplot(221, projection='3d')
ax2 = fig.add_subplot(222, projection='3d')
ax3 = fig.add_subplot(223, projection='3d')
ax4 = fig.add_subplot(224)

# Create meshgrid for 3D surface
n_obs_range = np.arange(1, 31)
theta_range = np.linspace(0, 1, 100)
N_OBS, THETA = np.meshgrid(n_obs_range, theta_range)

# Surface 1: Posterior mean as function of n_obs and n_detections
detection_ratios = [0.2, 0.33, 0.5, 0.67]
for ratio in detection_ratios:
    posterior_means = np.zeros_like(N_OBS, dtype=float)
    for i, n in enumerate(n_obs_range):
        n_det = int(n * ratio)
        post_alpha = bayesian_model.prior_alpha + n_det
        post_beta = bayesian_model.prior_beta + (n - n_det)
        mean = post_alpha / (post_alpha + post_beta)
        posterior_means[:, i] = mean
    
    ax1.plot_surface(N_OBS, THETA * 0 + ratio, posterior_means, alpha=0.4, 
                     label=f'Detection rate={ratio:.0%}')

ax1.set_xlabel('Number of Observations', fontsize=10)
ax1.set_ylabel('Detection Rate', fontsize=10)
ax1.set_zlabel('Posterior Mean', fontsize=10)
ax1.set_title('Posterior Mean vs. Observations', fontsize=12, fontweight='bold')
ax1.view_init(elev=25, azim=45)

# Surface 2: Posterior variance (uncertainty)
Z_variance = np.zeros_like(N_OBS, dtype=float)
for i, n in enumerate(n_obs_range):
    for j, theta in enumerate(theta_range):
        n_det = int(n * theta)
        post_alpha = bayesian_model.prior_alpha + n_det
        post_beta = bayesian_model.prior_beta + (n - n_det)
        variance = (post_alpha * post_beta) / \
                   ((post_alpha + post_beta)**2 * (post_alpha + post_beta + 1))
        Z_variance[j, i] = variance

surf2 = ax2.plot_surface(N_OBS, THETA, Z_variance, cmap='coolwarm', alpha=0.8)
ax2.set_xlabel('Number of Observations', fontsize=10)
ax2.set_ylabel(r'True $\theta$', fontsize=10)
ax2.set_zlabel('Posterior Variance', fontsize=10)
ax2.set_title('Uncertainty Reduction with Data', fontsize=12, fontweight='bold')
ax2.view_init(elev=25, azim=135)

# Surface 3: 95% Credible Interval Width
Z_ci_width = np.zeros_like(N_OBS, dtype=float)
for i, n in enumerate(n_obs_range):
    for j, theta in enumerate(theta_range):
        n_det = int(n * theta)
        post_alpha = bayesian_model.prior_alpha + n_det
        post_beta = bayesian_model.prior_beta + (n - n_det)
        ci_lower = beta.ppf(0.025, post_alpha, post_beta)
        ci_upper = beta.ppf(0.975, post_alpha, post_beta)
        Z_ci_width[j, i] = ci_upper - ci_lower

surf3 = ax3.plot_surface(N_OBS, THETA, Z_ci_width, cmap='viridis', alpha=0.8)
ax3.set_xlabel('Number of Observations', fontsize=10)
ax3.set_ylabel(r'True $\theta$', fontsize=10)
ax3.set_zlabel('95% CI Width', fontsize=10)
ax3.set_title('Credible Interval Width', fontsize=12, fontweight='bold')
ax3.view_init(elev=25, azim=45)

# Plot 4: Heatmap of posterior mean
n_obs_heatmap = np.arange(1, 51)
theta_heatmap = np.linspace(0, 1, 50)
posterior_mean_grid = np.zeros((len(theta_heatmap), len(n_obs_heatmap)))

for i, n in enumerate(n_obs_heatmap):
    for j, theta in enumerate(theta_heatmap):
        n_det = int(n * theta)
        post_alpha = bayesian_model.prior_alpha + n_det
        post_beta = bayesian_model.prior_beta + (n - n_det)
        posterior_mean_grid[j, i] = post_alpha / (post_alpha + post_beta)

im = ax4.imshow(posterior_mean_grid, aspect='auto', origin='lower', 
                extent=[1, 50, 0, 1], cmap='RdYlGn', vmin=0, vmax=1)
ax4.set_xlabel('Number of Observations', fontsize=11)
ax4.set_ylabel('Detection Rate', fontsize=11)
ax4.set_title('Posterior Mean Heatmap', fontsize=12, fontweight='bold')
cbar = plt.colorbar(im, ax=ax4)
cbar.set_label('Posterior Mean', fontsize=10)

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

# Visualization 4: Sensitivity Analysis
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# Different prior beliefs
priors_to_test = [
    (1, 1, "Uniform (uninformative)"),
    (2, 10, "Skeptical (our example)"),
    (5, 5, "Neutral"),
    (10, 2, "Optimistic")
]

for prior_params in priors_to_test:
    alpha, beta_param, label = prior_params
    model = BayesianLifeProbability(alpha, beta_param)
    posterior = model.posterior_distribution(theta_values, n_observations, n_detections)
    axes[0, 0].plot(theta_values, posterior, linewidth=2.5, label=label, alpha=0.8)

axes[0, 0].set_xlabel(r'Probability of Life ($\theta$)', fontsize=12)
axes[0, 0].set_ylabel('Posterior Density', fontsize=12)
axes[0, 0].set_title('Sensitivity to Prior Beliefs', fontsize=14, fontweight='bold')
axes[0, 0].legend(fontsize=11)
axes[0, 0].grid(True, alpha=0.3)

# Different amounts of data (same detection rate)
data_amounts = [(5, 2), (10, 3), (15, 5), (30, 10), (60, 20)]
for n, k in data_amounts:
    posterior = bayesian_model.posterior_distribution(theta_values, n, k)
    axes[0, 1].plot(theta_values, posterior, linewidth=2.5, 
                    label=f'n={n}, k={k}', alpha=0.8)

axes[0, 1].set_xlabel(r'Probability of Life ($\theta$)', fontsize=12)
axes[0, 1].set_ylabel('Posterior Density', fontsize=12)
axes[0, 1].set_title('Convergence with More Data (~33% detection rate)', 
                     fontsize=14, fontweight='bold')
axes[0, 1].legend(fontsize=11)
axes[0, 1].grid(True, alpha=0.3)

# Posterior mean vs. sample size
sample_sizes = np.arange(1, 101)
posterior_means = []
ci_widths = []

for n in sample_sizes:
    k = int(n * 0.33)  # 33% detection rate
    stats = bayesian_model.posterior_statistics(n, k)
    posterior_means.append(stats['mean'])
    ci_widths.append(stats['ci_95'][1] - stats['ci_95'][0])

axes[1, 0].plot(sample_sizes, posterior_means, 'b-', linewidth=2.5, label='Posterior Mean')
axes[1, 0].axhline(0.33, color='red', linestyle='--', linewidth=2, label='True Rate (MLE)')
axes[1, 0].fill_between(sample_sizes, 
                        np.array(posterior_means) - np.array(ci_widths)/2,
                        np.array(posterior_means) + np.array(ci_widths)/2,
                        alpha=0.3, color='blue', label='95% CI')
axes[1, 0].set_xlabel('Sample Size', fontsize=12)
axes[1, 0].set_ylabel('Posterior Mean', fontsize=12)
axes[1, 0].set_title('Convergence of Posterior Mean', fontsize=14, fontweight='bold')
axes[1, 0].legend(fontsize=11)
axes[1, 0].grid(True, alpha=0.3)

# Uncertainty reduction
axes[1, 1].plot(sample_sizes, ci_widths, 'purple', linewidth=2.5)
axes[1, 1].set_xlabel('Sample Size', fontsize=12)
axes[1, 1].set_ylabel('95% Credible Interval Width', fontsize=12)
axes[1, 1].set_title('Uncertainty Reduction with Sample Size', fontsize=14, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)

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

print("\nVisualization complete! All graphs have been generated.")