import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from scipy.spatial.distance import euclidean
import random
from typing import List, Tuple
import copy

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

class RoverPathOptimizer:
    """
    Optimizes rover exploration routes using genetic algorithms
    to maximize scientific value while minimizing fuel and time.
    """
    
    def __init__(self, start_pos: Tuple[float, float], 
                 poi_positions: List[Tuple[float, float]],
                 poi_values: List[float],
                 max_distance: float,
                 max_time: float,
                 sampling_time: float = 2.0):
        """
        Initialize the optimizer.
        
        Parameters:
        - start_pos: Starting position (x, y)
        - poi_positions: List of points of interest coordinates
        - poi_values: Scientific value of each POI (0-100)
        - max_distance: Maximum travel distance (fuel constraint)
        - max_time: Maximum mission time
        - sampling_time: Time required to sample each POI
        """
        self.start = start_pos
        self.pois = poi_positions
        self.values = np.array(poi_values)
        self.max_dist = max_distance
        self.max_time = max_time
        self.sampling_time = sampling_time
        self.n_pois = len(poi_positions)
        self.rover_speed = 1.0  # units per time unit
        
    def calculate_distance(self, route: List[int]) -> float:
        """Calculate total distance for a given route."""
        if len(route) == 0:
            return 0.0
        
        total_dist = euclidean(self.start, self.pois[route[0]])
        for i in range(len(route) - 1):
            total_dist += euclidean(self.pois[route[i]], self.pois[route[i+1]])
        return total_dist
    
    def calculate_time(self, route: List[int]) -> float:
        """Calculate total mission time."""
        travel_time = self.calculate_distance(route) / self.rover_speed
        sampling_time = len(route) * self.sampling_time
        return travel_time + sampling_time
    
    def calculate_value(self, route: List[int]) -> float:
        """Calculate total scientific value."""
        if len(route) == 0:
            return 0.0
        return np.sum(self.values[route])
    
    def is_feasible(self, route: List[int]) -> bool:
        """Check if route satisfies constraints."""
        distance = self.calculate_distance(route)
        time = self.calculate_time(route)
        return distance <= self.max_dist and time <= self.max_time
    
    def fitness(self, route: List[int]) -> float:
        """
        Calculate fitness score (higher is better).
        Combines scientific value with penalties for constraint violations.
        """
        if len(route) == 0:
            return 0.0
        
        value = self.calculate_value(route)
        distance = self.calculate_distance(route)
        time = self.calculate_time(route)
        
        # Penalty for constraint violations
        distance_penalty = max(0, distance - self.max_dist) * 100
        time_penalty = max(0, time - self.max_time) * 100
        
        # Fitness = value - normalized_distance - penalties
        fitness = value - (distance / self.max_dist) * 20 - distance_penalty - time_penalty
        
        return fitness
    
    def generate_random_route(self) -> List[int]:
        """Generate a random feasible route."""
        all_indices = list(range(self.n_pois))
        random.shuffle(all_indices)
        
        # Start with empty route and add POIs until constraints violated
        route = []
        for idx in all_indices:
            test_route = route + [idx]
            if self.is_feasible(test_route):
                route = test_route
            else:
                break
        
        return route
    
    def mutate(self, route: List[int]) -> List[int]:
        """Apply mutation to a route."""
        if len(route) < 2:
            return route
        
        new_route = route.copy()
        mutation_type = random.choice(['swap', 'reverse', 'insert', 'remove'])
        
        if mutation_type == 'swap' and len(new_route) >= 2:
            i, j = random.sample(range(len(new_route)), 2)
            new_route[i], new_route[j] = new_route[j], new_route[i]
        
        elif mutation_type == 'reverse' and len(new_route) >= 2:
            i, j = sorted(random.sample(range(len(new_route)), 2))
            new_route[i:j+1] = reversed(new_route[i:j+1])
        
        elif mutation_type == 'insert':
            # Try to insert a new POI
            available = [i for i in range(self.n_pois) if i not in new_route]
            if available:
                new_poi = random.choice(available)
                pos = random.randint(0, len(new_route))
                new_route.insert(pos, new_poi)
        
        elif mutation_type == 'remove' and len(new_route) > 0:
            pos = random.randint(0, len(new_route) - 1)
            new_route.pop(pos)
        
        return new_route if self.is_feasible(new_route) else route
    
    def crossover(self, parent1: List[int], parent2: List[int]) -> List[int]:
        """Perform ordered crossover."""
        if len(parent1) == 0 or len(parent2) == 0:
            return parent1 if len(parent1) > len(parent2) else parent2
        
        # Take a subset from parent1
        size = min(len(parent1), len(parent2))
        start = random.randint(0, size - 1)
        end = random.randint(start + 1, size)
        
        child = parent1[start:end]
        
        # Add elements from parent2 that are not in child
        for gene in parent2:
            if gene not in child:
                child.append(gene)
        
        return child if self.is_feasible(child) else parent1
    
    def genetic_algorithm(self, population_size: int = 100, 
                         generations: int = 200,
                         mutation_rate: float = 0.2) -> Tuple[List[int], List[float]]:
        """
        Run genetic algorithm to find optimal route.
        
        Returns:
        - best_route: Optimal route found
        - history: Fitness history over generations
        """
        # Initialize population
        population = [self.generate_random_route() for _ in range(population_size)]
        best_fitness_history = []
        
        for gen in range(generations):
            # Evaluate fitness
            fitness_scores = [self.fitness(route) for route in population]
            
            # Track best solution
            best_idx = np.argmax(fitness_scores)
            best_fitness_history.append(fitness_scores[best_idx])
            
            # Selection (tournament selection)
            new_population = []
            for _ in range(population_size):
                # Tournament
                tournament_size = 5
                tournament = random.sample(list(zip(population, fitness_scores)), 
                                         min(tournament_size, len(population)))
                winner = max(tournament, key=lambda x: x[1])[0]
                new_population.append(winner.copy())
            
            # Crossover and mutation
            next_generation = [population[best_idx].copy()]  # Elitism
            
            for i in range(1, population_size):
                parent1, parent2 = random.sample(new_population, 2)
                
                # Crossover
                if random.random() < 0.7:
                    child = self.crossover(parent1, parent2)
                else:
                    child = parent1.copy()
                
                # Mutation
                if random.random() < mutation_rate:
                    child = self.mutate(child)
                
                next_generation.append(child)
            
            population = next_generation
        
        # Return best solution
        final_fitness = [self.fitness(route) for route in population]
        best_idx = np.argmax(final_fitness)
        
        return population[best_idx], best_fitness_history


# ============================================================================
# Example Problem Setup
# ============================================================================

print("=" * 70)
print("ROVER EXPLORATION ROUTE OPTIMIZATION")
print("=" * 70)
print()

# Define the problem
start_position = (0, 0)

# 15 Points of Interest with varying scientific values
poi_positions = [
    (5, 8), (12, 15), (18, 6), (25, 20), (8, 3),
    (15, 10), (22, 14), (10, 18), (28, 8), (6, 12),
    (20, 4), (14, 22), (30, 16), (4, 15), (26, 11)
]

# Scientific values (0-100): higher means more interesting
poi_values = [85, 92, 70, 95, 60, 88, 75, 80, 65, 78, 
              72, 90, 98, 68, 82]

# Constraints
max_distance = 120.0  # Maximum fuel allows 120 distance units
max_time = 150.0      # Maximum mission time

print("Problem Configuration:")
print(f"  Start Position: {start_position}")
print(f"  Number of POIs: {len(poi_positions)}")
print(f"  Maximum Distance: {max_distance} units")
print(f"  Maximum Mission Time: {max_time} time units")
print(f"  Sampling Time per POI: 2.0 time units")
print()

# Create optimizer
optimizer = RoverPathOptimizer(
    start_pos=start_position,
    poi_positions=poi_positions,
    poi_values=poi_values,
    max_distance=max_distance,
    max_time=max_time,
    sampling_time=2.0
)

# Run optimization
print("Running Genetic Algorithm...")
print("  Population Size: 100")
print("  Generations: 200")
print("  Mutation Rate: 0.2")
print()

best_route, fitness_history = optimizer.genetic_algorithm(
    population_size=100,
    generations=200,
    mutation_rate=0.2
)

# Calculate metrics for best route
total_distance = optimizer.calculate_distance(best_route)
total_time = optimizer.calculate_time(best_route)
total_value = optimizer.calculate_value(best_route)

print("=" * 70)
print("OPTIMIZATION RESULTS")
print("=" * 70)
print()
print(f"Best Route Found: {best_route}")
print(f"Number of POIs Visited: {len(best_route)}")
print(f"Total Scientific Value: {total_value:.2f} / {sum(poi_values):.2f}")
print(f"Total Distance: {total_distance:.2f} / {max_distance:.2f} units")
print(f"Total Mission Time: {total_time:.2f} / {max_time:.2f} units")
print(f"Distance Utilization: {(total_distance/max_distance)*100:.1f}%")
print(f"Time Utilization: {(total_time/max_time)*100:.1f}%")
print(f"Value Efficiency: {(total_value/sum(poi_values))*100:.1f}%")
print()

# ============================================================================
# Visualization
# ============================================================================

fig, axes = plt.subplots(2, 2, figsize=(16, 14))

# Plot 1: Route Visualization
ax1 = axes[0, 0]
ax1.set_title('Optimal Rover Exploration Route', fontsize=14, fontweight='bold')
ax1.set_xlabel('X Coordinate (km)', fontsize=11)
ax1.set_ylabel('Y Coordinate (km)', fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')

# Plot all POIs
for i, (pos, value) in enumerate(zip(poi_positions, poi_values)):
    if i in best_route:
        color = 'green'
        size = 200
        alpha = 0.7
    else:
        color = 'lightgray'
        size = 100
        alpha = 0.4
    
    ax1.scatter(pos[0], pos[1], c=color, s=size, alpha=alpha, edgecolors='black', linewidth=1.5)
    ax1.text(pos[0], pos[1]-1.2, f'POI{i}\n({value})', 
             ha='center', va='top', fontsize=8, fontweight='bold')

# Plot start position
ax1.scatter(start_position[0], start_position[1], c='red', s=300, marker='*', 
           edgecolors='black', linewidth=2, label='Start', zorder=5)

# Plot route
if len(best_route) > 0:
    route_coords = [start_position] + [poi_positions[i] for i in best_route]
    route_x = [coord[0] for coord in route_coords]
    route_y = [coord[1] for coord in route_coords]
    ax1.plot(route_x, route_y, 'b-', linewidth=2, alpha=0.6, label='Route')
    
    # Add arrows
    for i in range(len(route_coords) - 1):
        dx = route_coords[i+1][0] - route_coords[i][0]
        dy = route_coords[i+1][1] - route_coords[i][1]
        ax1.arrow(route_coords[i][0], route_coords[i][1], dx*0.8, dy*0.8,
                 head_width=0.8, head_length=0.5, fc='blue', ec='blue', alpha=0.5)

ax1.legend(loc='upper right', fontsize=10)

# Plot 2: Fitness Evolution
ax2 = axes[0, 1]
ax2.set_title('Genetic Algorithm Convergence', fontsize=14, fontweight='bold')
ax2.set_xlabel('Generation', fontsize=11)
ax2.set_ylabel('Fitness Score', fontsize=11)
ax2.plot(fitness_history, linewidth=2, color='darkblue')
ax2.grid(True, alpha=0.3)
ax2.fill_between(range(len(fitness_history)), fitness_history, alpha=0.3, color='lightblue')

# Plot 3: POI Values Comparison
ax3 = axes[1, 0]
ax3.set_title('Scientific Value: Visited vs Unvisited POIs', fontsize=14, fontweight='bold')
ax3.set_xlabel('POI Index', fontsize=11)
ax3.set_ylabel('Scientific Value', fontsize=11)

visited_values = [poi_values[i] if i in best_route else 0 for i in range(len(poi_values))]
unvisited_values = [poi_values[i] if i not in best_route else 0 for i in range(len(poi_values))]

x = np.arange(len(poi_values))
width = 0.8

ax3.bar(x, visited_values, width, label='Visited', color='green', alpha=0.7)
ax3.bar(x, unvisited_values, width, bottom=visited_values, label='Unvisited', 
        color='lightgray', alpha=0.7)
ax3.set_xticks(x)
ax3.set_xticklabels([f'POI{i}' for i in range(len(poi_values))], rotation=45, ha='right')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')

# Plot 4: Resource Utilization
ax4 = axes[1, 1]
ax4.set_title('Resource Utilization Analysis', fontsize=14, fontweight='bold')

categories = ['Distance\n(Fuel)', 'Time', 'Scientific\nValue']
utilization = [
    (total_distance / max_distance) * 100,
    (total_time / max_time) * 100,
    (total_value / sum(poi_values)) * 100
]
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']

bars = ax4.barh(categories, utilization, color=colors, alpha=0.7, edgecolor='black', linewidth=1.5)
ax4.set_xlabel('Utilization (%)', fontsize=11)
ax4.set_xlim(0, 110)
ax4.axvline(x=100, color='red', linestyle='--', linewidth=2, label='Maximum Capacity')
ax4.grid(True, alpha=0.3, axis='x')

# Add percentage labels
for i, (bar, val) in enumerate(zip(bars, utilization)):
    ax4.text(val + 2, bar.get_y() + bar.get_height()/2, f'{val:.1f}%', 
            va='center', fontsize=11, fontweight='bold')

ax4.legend(fontsize=10)

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

print()
print("=" * 70)
print("Visualization saved as 'rover_optimization_results.png'")
print("=" * 70)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import pandas as pd

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

# ============================================================================
# STEP 1: Generate Exoplanet Dataset
# ============================================================================

n_planets = 20  # Number of candidate exoplanets

# Planet names
planet_names = [f"Kepler-{100+i}b" for i in range(n_planets)]

# Generate planet characteristics
# Distance from Earth (light-years): affects observation time
distances = np.random.uniform(10, 100, n_planets)

# Star type factor (M-dwarf=1.5, G-type=1.0, K-type=1.2)
star_types = np.random.choice([1.0, 1.2, 1.5], n_planets, p=[0.3, 0.4, 0.3])

# Planet size (Earth radii): affects signal strength
planet_sizes = np.random.uniform(0.8, 2.5, n_planets)

# Habitable zone score (0-1): 1 = center of HZ, 0 = edge
hz_score = np.random.uniform(0.3, 1.0, n_planets)

# Atmospheric thickness indicator (0-1)
atmosphere_score = np.random.uniform(0.4, 1.0, n_planets)

# Calculate observation time required (hours)
# Farther planets and smaller planets need more time
observation_time = (distances / 10) * (2.0 / planet_sizes) * star_types
observation_time = observation_time.round(1)

# Calculate life probability based on multiple factors
# This is a simplified model combining various habitability indicators
life_probability = (
    hz_score * 0.4 +  # Habitable zone position
    atmosphere_score * 0.3 +  # Atmospheric presence
    (1 / (planet_sizes ** 0.5)) * 0.15 +  # Size factor (Earth-like preferred)
    (1 / (distances / 50)) * 0.15  # Proximity bonus
)
# Normalize to 0-1 range
life_probability = (life_probability - life_probability.min()) / (life_probability.max() - life_probability.min())
life_probability = (life_probability * 0.6 + 0.1).round(3)  # Scale to 0.1-0.7 range

# ============================================================================
# STEP 2: Optimization Algorithm (Dynamic Programming for 0/1 Knapsack)
# ============================================================================

def knapsack_optimization(life_probs, obs_times, max_time):
    """
    Solve the 0/1 knapsack problem using dynamic programming.
    
    Parameters:
    -----------
    life_probs : array-like
        Life probability for each planet
    obs_times : array-like
        Observation time required for each planet (hours)
    max_time : float
        Maximum available telescope time (hours)
    
    Returns:
    --------
    selected_indices : list
        Indices of selected planets
    total_prob : float
        Total expected life detection probability
    total_time : float
        Total observation time used
    """
    n = len(life_probs)
    # Convert observation times to integers (multiply by 10 to preserve 1 decimal)
    times_int = (obs_times * 10).astype(int)
    max_time_int = int(max_time * 10)
    
    # Scale probabilities to avoid floating point issues
    probs_scaled = (life_probs * 1000).astype(int)
    
    # DP table: dp[i][w] = maximum probability using first i items with capacity w
    dp = np.zeros((n + 1, max_time_int + 1), dtype=int)
    
    # Fill the DP table
    for i in range(1, n + 1):
        for w in range(max_time_int + 1):
            # Don't include item i-1
            dp[i][w] = dp[i-1][w]
            
            # Include item i-1 if it fits
            if times_int[i-1] <= w:
                include_value = dp[i-1][w - times_int[i-1]] + probs_scaled[i-1]
                dp[i][w] = max(dp[i][w], include_value)
    
    # Backtrack to find selected items
    selected = []
    w = max_time_int
    for i in range(n, 0, -1):
        if dp[i][w] != dp[i-1][w]:
            selected.append(i - 1)
            w -= times_int[i - 1]
    
    selected.reverse()
    
    total_prob = sum(life_probs[selected])
    total_time = sum(obs_times[selected])
    
    return selected, total_prob, total_time

# ============================================================================
# STEP 3: Greedy Algorithm for Comparison
# ============================================================================

def greedy_optimization(life_probs, obs_times, max_time):
    """
    Greedy algorithm: select planets by efficiency (prob/time ratio).
    """
    n = len(life_probs)
    efficiency = life_probs / obs_times
    sorted_indices = np.argsort(efficiency)[::-1]
    
    selected = []
    total_time = 0
    
    for idx in sorted_indices:
        if total_time + obs_times[idx] <= max_time:
            selected.append(idx)
            total_time += obs_times[idx]
    
    total_prob = sum(life_probs[selected])
    
    return selected, total_prob, total_time

# ============================================================================
# STEP 4: Run Optimization
# ============================================================================

# Available telescope time (hours)
MAX_TELESCOPE_TIME = 100.0

print("=" * 80)
print("EXOPLANET OBSERVATION OPTIMIZATION PROBLEM")
print("=" * 80)
print(f"\nTotal Available Telescope Time: {MAX_TELESCOPE_TIME} hours")
print(f"Number of Candidate Exoplanets: {n_planets}")
print(f"Total Time Needed to Observe All: {observation_time.sum():.1f} hours")
print("\n")

# Solve using Dynamic Programming (Optimal)
selected_optimal, prob_optimal, time_optimal = knapsack_optimization(
    life_probability, observation_time, MAX_TELESCOPE_TIME
)

# Solve using Greedy Algorithm (for comparison)
selected_greedy, prob_greedy, time_greedy = greedy_optimization(
    life_probability, observation_time, MAX_TELESCOPE_TIME
)

print("OPTIMIZATION RESULTS")
print("-" * 80)
print(f"\n[OPTIMAL SOLUTION - Dynamic Programming]")
print(f"  Number of planets selected: {len(selected_optimal)}")
print(f"  Total expected life probability: {prob_optimal:.4f}")
print(f"  Total observation time used: {time_optimal:.1f} / {MAX_TELESCOPE_TIME} hours")
print(f"  Telescope utilization: {(time_optimal/MAX_TELESCOPE_TIME)*100:.1f}%")

print(f"\n[GREEDY SOLUTION - For Comparison]")
print(f"  Number of planets selected: {len(selected_greedy)}")
print(f"  Total expected life probability: {prob_greedy:.4f}")
print(f"  Total observation time used: {time_greedy:.1f} / {MAX_TELESCOPE_TIME} hours")
print(f"  Telescope utilization: {(time_greedy/MAX_TELESCOPE_TIME)*100:.1f}%")

print(f"\n[PERFORMANCE IMPROVEMENT]")
print(f"  Probability gain: {((prob_optimal - prob_greedy) / prob_greedy * 100):.2f}%")
print("\n")

# ============================================================================
# STEP 5: Detailed Results Table
# ============================================================================

# Create comprehensive dataframe
df = pd.DataFrame({
    'Planet': planet_names,
    'Distance (ly)': distances.round(1),
    'Size (R_Earth)': planet_sizes.round(2),
    'HZ Score': hz_score.round(2),
    'Atmosphere': atmosphere_score.round(2),
    'Life Prob': life_probability,
    'Obs Time (h)': observation_time,
    'Efficiency': (life_probability / observation_time).round(4)
})

# Add selection flags
df['Optimal'] = ['✓' if i in selected_optimal else '' for i in range(n_planets)]
df['Greedy'] = ['✓' if i in selected_greedy else '' for i in range(n_planets)]

print("CANDIDATE EXOPLANETS - DETAILED DATA")
print("-" * 80)
print(df.to_string(index=False))
print("\n")

# ============================================================================
# STEP 6: Visualization
# ============================================================================

fig = plt.figure(figsize=(16, 12))

# Plot 1: Life Probability vs Observation Time (Scatter)
ax1 = plt.subplot(2, 3, 1)
colors = ['red' if i in selected_optimal else 'lightgray' for i in range(n_planets)]
sizes = [200 if i in selected_optimal else 50 for i in range(n_planets)]
scatter = ax1.scatter(observation_time, life_probability, c=colors, s=sizes, alpha=0.6, edgecolors='black')

# Add labels for selected planets
for i in selected_optimal:
    ax1.annotate(planet_names[i], (observation_time[i], life_probability[i]), 
                fontsize=8, xytext=(5, 5), textcoords='offset points')

ax1.set_xlabel('Observation Time Required (hours)', fontsize=11, fontweight='bold')
ax1.set_ylabel('Life Detection Probability', fontsize=11, fontweight='bold')
ax1.set_title('Planet Selection Map\n(Red = Selected by Optimal Algorithm)', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Plot 2: Efficiency Ranking
ax2 = plt.subplot(2, 3, 2)
efficiency = life_probability / observation_time
sorted_idx = np.argsort(efficiency)[::-1]
colors_eff = ['green' if i in selected_optimal else 'lightblue' for i in sorted_idx]
bars = ax2.barh(range(n_planets), efficiency[sorted_idx], color=colors_eff, edgecolor='black')
ax2.set_yticks(range(n_planets))
ax2.set_yticklabels([planet_names[i] for i in sorted_idx], fontsize=8)
ax2.set_xlabel('Efficiency (Probability / Hour)', fontsize=11, fontweight='bold')
ax2.set_title('Planet Efficiency Ranking\n(Green = Selected)', fontsize=12, fontweight='bold')
ax2.invert_yaxis()
ax2.grid(True, axis='x', alpha=0.3)

# Plot 3: Time Allocation Comparison
ax3 = plt.subplot(2, 3, 3)
methods = ['Optimal\n(DP)', 'Greedy', 'Unused']
times = [time_optimal, time_greedy, MAX_TELESCOPE_TIME - time_optimal]
colors_bar = ['#2ecc71', '#3498db', '#ecf0f1']
bars = ax3.bar(methods, [time_optimal, time_greedy, 0], color=colors_bar[:2], edgecolor='black', linewidth=2)
ax3.bar(['Unused'], [MAX_TELESCOPE_TIME - time_optimal], bottom=[time_optimal], 
        color=colors_bar[2], edgecolor='black', linewidth=2)
ax3.bar(['Unused'], [MAX_TELESCOPE_TIME - time_greedy], bottom=[time_greedy], 
        color=colors_bar[2], edgecolor='black', linewidth=2)
ax3.axhline(MAX_TELESCOPE_TIME, color='red', linestyle='--', linewidth=2, label='Max Time')
ax3.set_ylabel('Telescope Time (hours)', fontsize=11, fontweight='bold')
ax3.set_title('Telescope Time Utilization', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, axis='y', alpha=0.3)

# Plot 4: Expected Life Probability Comparison
ax4 = plt.subplot(2, 3, 4)
methods = ['Optimal (DP)', 'Greedy']
probs = [prob_optimal, prob_greedy]
colors_prob = ['#e74c3c', '#f39c12']
bars = ax4.bar(methods, probs, color=colors_prob, edgecolor='black', linewidth=2)
ax4.set_ylabel('Total Expected Life Probability', fontsize=11, fontweight='bold')
ax4.set_title('Total Expected Life Detection Probability', fontsize=12, fontweight='bold')
for i, (bar, prob) in enumerate(zip(bars, probs)):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
            f'{prob:.4f}', ha='center', va='bottom', fontsize=12, fontweight='bold')
ax4.grid(True, axis='y', alpha=0.3)

# Plot 5: Selected Planets Breakdown
ax5 = plt.subplot(2, 3, 5)
selected_names = [planet_names[i] for i in selected_optimal]
selected_probs = [life_probability[i] for i in selected_optimal]
selected_times = [observation_time[i] for i in selected_optimal]

# Create stacked time bar
cumulative_time = 0
colors_stack = plt.cm.viridis(np.linspace(0.2, 0.9, len(selected_optimal)))
for i, (name, time_val, prob) in enumerate(zip(selected_names, selected_times, selected_probs)):
    ax5.barh(0, time_val, left=cumulative_time, color=colors_stack[i], 
            edgecolor='black', linewidth=1.5, label=f'{name}\n(P={prob:.3f})')
    # Add text in the middle of each segment
    ax5.text(cumulative_time + time_val/2, 0, f'{time_val:.1f}h', 
            ha='center', va='center', fontsize=9, fontweight='bold')
    cumulative_time += time_val

ax5.set_xlim(0, MAX_TELESCOPE_TIME)
ax5.set_ylim(-0.5, 0.5)
ax5.set_xlabel('Cumulative Observation Time (hours)', fontsize=11, fontweight='bold')
ax5.set_yticks([])
ax5.set_title('Optimal Observation Schedule Timeline', fontsize=12, fontweight='bold')
ax5.legend(loc='upper left', bbox_to_anchor=(1.02, 1), fontsize=8, ncol=1)
ax5.grid(True, axis='x', alpha=0.3)

# Plot 6: Planet Characteristics Heatmap
ax6 = plt.subplot(2, 3, 6)
selected_chars = np.array([
    [hz_score[i] for i in selected_optimal],
    [atmosphere_score[i] for i in selected_optimal],
    [1 - (distances[i] / 100) for i in selected_optimal],  # Proximity (inverted distance)
    [planet_sizes[i] / 2.5 for i in selected_optimal],  # Size (normalized)
]).T

im = ax6.imshow(selected_chars, cmap='YlOrRd', aspect='auto', vmin=0, vmax=1)
ax6.set_xticks(range(4))
ax6.set_xticklabels(['HZ\nScore', 'Atmosphere', 'Proximity', 'Size\n(norm)'], fontsize=9)
ax6.set_yticks(range(len(selected_optimal)))
ax6.set_yticklabels([planet_names[i] for i in selected_optimal], fontsize=9)
ax6.set_title('Selected Planets Characteristics Heatmap', fontsize=12, fontweight='bold')

# Add colorbar
cbar = plt.colorbar(im, ax=ax6)
cbar.set_label('Normalized Score', fontsize=10)

# Add values in cells
for i in range(len(selected_optimal)):
    for j in range(4):
        text = ax6.text(j, i, f'{selected_chars[i, j]:.2f}',
                       ha="center", va="center", color="black", fontsize=8, fontweight='bold')

plt.tight_layout()
plt.savefig('exoplanet_optimization_results.png', dpi=150, bbox_inches='tight')
plt.show()

print("=" * 80)
print("VISUALIZATION COMPLETE")
print("=" * 80)
print("Graph saved as: exoplanet_optimization_results.png")
print("\n")

# ============================================================================
# STEP 7: Summary Statistics
# ============================================================================

print("SUMMARY STATISTICS")
print("-" * 80)
print(f"\nSelected Planets (Optimal Solution):")
for idx in selected_optimal:
    print(f"  • {planet_names[idx]}: "
          f"Life Prob = {life_probability[idx]:.3f}, "
          f"Time = {observation_time[idx]:.1f}h, "
          f"Efficiency = {(life_probability[idx]/observation_time[idx]):.4f}")

print(f"\nPlanets Not Selected (Top 5 by Efficiency):")
not_selected = [i for i in range(n_planets) if i not in selected_optimal]
not_selected_eff = [(i, life_probability[i]/observation_time[i]) for i in not_selected]
not_selected_eff.sort(key=lambda x: x[1], reverse=True)
for idx, eff in not_selected_eff[:5]:
    print(f"  • {planet_names[idx]}: "
          f"Life Prob = {life_probability[idx]:.3f}, "
          f"Time = {observation_time[idx]:.1f}h, "
          f"Efficiency = {eff:.4f}")

print("\n" + "=" * 80)
print("ANALYSIS COMPLETE")
print("=" * 80)

dp[i][w] = max(dp[i-1][w], dp[i-1][w - times_int[i-1]] + probs_scaled[i-1])

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.patches import Rectangle
import seaborn as sns

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

# Define the number of observation targets
n_targets = 20

# Generate synthetic data for astronomical targets
target_names = [f"Target-{i+1:02d}" for i in range(n_targets)]

# Observation costs (in hours, ranging from 0.5 to 8 hours)
observation_costs = np.random.uniform(0.5, 8.0, n_targets)

# Scientific importance scores (0-100 scale)
scientific_importance = np.random.uniform(30, 100, n_targets)

# Visibility scores (0-1 scale, representing altitude and atmospheric quality)
# Higher values mean better observing conditions
visibility_scores = np.random.uniform(0.3, 1.0, n_targets)

# Cost penalty parameter
alpha = 1.0

# Calculate the value function for each target
# Value = (Scientific_Importance × Visibility) / Cost^alpha
def calculate_value(science, visibility, cost, alpha=1.0):
    """
    Calculate the optimization value for an observation target.
    
    Parameters:
    -----------
    science : float or array
        Scientific importance score (0-100)
    visibility : float or array
        Visibility score (0-1)
    cost : float or array
        Observation cost (hours)
    alpha : float
        Cost penalty exponent
    
    Returns:
    --------
    value : float or array
        Calculated value metric
    """
    return (science * visibility) / (cost ** alpha)

values = calculate_value(scientific_importance, visibility_scores, 
                         observation_costs, alpha)

# Create a DataFrame for easy analysis
df = pd.DataFrame({
    'Target': target_names,
    'Cost (hours)': observation_costs,
    'Scientific_Importance': scientific_importance,
    'Visibility': visibility_scores,
    'Value': values
})

# Sort by value (descending)
df_sorted = df.sort_values('Value', ascending=False).reset_index(drop=True)
df_sorted['Rank'] = range(1, len(df_sorted) + 1)

# Display top 10 targets
print("=" * 80)
print("ASTRONOMICAL OBSERVATION TARGET SELECTION OPTIMIZATION")
print("=" * 80)
print("\nTop 10 Recommended Targets (sorted by optimization value):\n")
print(df_sorted[['Rank', 'Target', 'Cost (hours)', 'Scientific_Importance', 
                  'Visibility', 'Value']].head(10).to_string(index=False))
print("\n" + "=" * 80)

# Statistical summary
print("\nStatistical Summary:")
print("-" * 80)
print(f"Total number of targets: {n_targets}")
print(f"Cost penalty parameter (α): {alpha}")
print(f"\nValue metric statistics:")
print(f"  Mean:   {values.mean():.2f}")
print(f"  Median: {np.median(values):.2f}")
print(f"  Std:    {values.std():.2f}")
print(f"  Max:    {values.max():.2f} (Target: {df_sorted.iloc[0]['Target']})")
print(f"  Min:    {values.min():.2f} (Target: {df_sorted.iloc[-1]['Target']})")
print("=" * 80)

# ============================================================================
# VISUALIZATION
# ============================================================================

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

# Create figure with multiple subplots
fig = plt.figure(figsize=(16, 12))

# -------------------------------------------------------------------------
# Plot 1: 3D Scatter Plot - Cost vs Scientific Importance vs Visibility
# -------------------------------------------------------------------------
ax1 = fig.add_subplot(2, 3, 1, projection='3d')

scatter = ax1.scatter(observation_costs, scientific_importance, visibility_scores,
                      c=values, cmap='viridis', s=100, alpha=0.7, edgecolors='k')

ax1.set_xlabel('Observation Cost (hours)', fontsize=10, labelpad=10)
ax1.set_ylabel('Scientific Importance', fontsize=10, labelpad=10)
ax1.set_zlabel('Visibility Score', fontsize=10, labelpad=10)
ax1.set_title('3D Parameter Space\n(colored by Value)', fontsize=12, fontweight='bold')

cbar1 = plt.colorbar(scatter, ax=ax1, pad=0.1, shrink=0.6)
cbar1.set_label('Value', fontsize=9)

# -------------------------------------------------------------------------
# Plot 2: Bar Chart - Top 10 Targets by Value
# -------------------------------------------------------------------------
ax2 = fig.add_subplot(2, 3, 2)

top10 = df_sorted.head(10)
colors = plt.cm.viridis(np.linspace(0.3, 0.9, 10))

bars = ax2.barh(range(10), top10['Value'], color=colors, edgecolor='black', linewidth=1.2)
ax2.set_yticks(range(10))
ax2.set_yticklabels(top10['Target'], fontsize=9)
ax2.set_xlabel('Optimization Value', fontsize=11, fontweight='bold')
ax2.set_title('Top 10 Targets by Value', fontsize=12, fontweight='bold')
ax2.invert_yaxis()
ax2.grid(axis='x', alpha=0.3)

# Add value labels on bars
for i, (bar, val) in enumerate(zip(bars, top10['Value'])):
    ax2.text(val + 0.5, i, f'{val:.1f}', va='center', fontsize=8)

# -------------------------------------------------------------------------
# Plot 3: Scatter Plot - Cost vs Value
# -------------------------------------------------------------------------
ax3 = fig.add_subplot(2, 3, 3)

scatter3 = ax3.scatter(observation_costs, values, c=visibility_scores, 
                       cmap='coolwarm', s=150, alpha=0.7, edgecolors='k', linewidth=1)

ax3.set_xlabel('Observation Cost (hours)', fontsize=11, fontweight='bold')
ax3.set_ylabel('Optimization Value', fontsize=11, fontweight='bold')
ax3.set_title('Cost vs Value\n(colored by Visibility)', fontsize=12, fontweight='bold')
ax3.grid(alpha=0.3)

cbar3 = plt.colorbar(scatter3, ax=ax3)
cbar3.set_label('Visibility Score', fontsize=9)

# Highlight top 3 targets
top3_indices = df_sorted.head(3).index
for idx in top3_indices:
    ax3.annotate(df.iloc[idx]['Target'], 
                xy=(df.iloc[idx]['Cost (hours)'], df.iloc[idx]['Value']),
                xytext=(10, 10), textcoords='offset points',
                bbox=dict(boxstyle='round,pad=0.5', fc='yellow', alpha=0.7),
                arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0', lw=1.5),
                fontsize=8, fontweight='bold')

# -------------------------------------------------------------------------
# Plot 4: Heatmap - Correlation Matrix
# -------------------------------------------------------------------------
ax4 = fig.add_subplot(2, 3, 4)

correlation_matrix = df[['Cost (hours)', 'Scientific_Importance', 
                          'Visibility', 'Value']].corr()

sns.heatmap(correlation_matrix, annot=True, fmt='.2f', cmap='coolwarm', 
            center=0, square=True, linewidths=2, cbar_kws={"shrink": 0.8},
            ax=ax4, vmin=-1, vmax=1)

ax4.set_title('Correlation Matrix', fontsize=12, fontweight='bold')
ax4.set_xticklabels(['Cost', 'Science', 'Visibility', 'Value'], rotation=45, ha='right')
ax4.set_yticklabels(['Cost', 'Science', 'Visibility', 'Value'], rotation=0)

# -------------------------------------------------------------------------
# Plot 5: Distribution Plots
# -------------------------------------------------------------------------
ax5 = fig.add_subplot(2, 3, 5)

ax5.hist(values, bins=15, color='steelblue', alpha=0.7, edgecolor='black', linewidth=1.2)
ax5.axvline(values.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean = {values.mean():.2f}')
ax5.axvline(np.median(values), color='green', linestyle='--', linewidth=2, label=f'Median = {np.median(values):.2f}')

ax5.set_xlabel('Optimization Value', fontsize=11, fontweight='bold')
ax5.set_ylabel('Frequency', fontsize=11, fontweight='bold')
ax5.set_title('Value Distribution', fontsize=12, fontweight='bold')
ax5.legend(fontsize=9)
ax5.grid(axis='y', alpha=0.3)

# -------------------------------------------------------------------------
# Plot 6: Pareto Front Analysis (Cost vs Combined Score)
# -------------------------------------------------------------------------
ax6 = fig.add_subplot(2, 3, 6)

combined_score = scientific_importance * visibility_scores

# Plot all points
ax6.scatter(observation_costs, combined_score, c='lightgray', s=100, 
            alpha=0.5, edgecolors='k', linewidth=1, label='All Targets')

# Identify and plot Pareto-optimal points
# A point is Pareto-optimal if no other point has both lower cost AND higher score
pareto_points = []
for i in range(n_targets):
    is_pareto = True
    for j in range(n_targets):
        if i != j:
            # Check if point j dominates point i
            if (observation_costs[j] <= observation_costs[i] and 
                combined_score[j] >= combined_score[i] and
                (observation_costs[j] < observation_costs[i] or 
                 combined_score[j] > combined_score[i])):
                is_pareto = False
                break
    if is_pareto:
        pareto_points.append(i)

pareto_costs = observation_costs[pareto_points]
pareto_scores = combined_score[pareto_points]

# Sort for plotting the frontier
sort_idx = np.argsort(pareto_costs)
pareto_costs_sorted = pareto_costs[sort_idx]
pareto_scores_sorted = pareto_scores[sort_idx]

ax6.scatter(pareto_costs, pareto_scores, c='red', s=200, 
            marker='*', edgecolors='black', linewidth=1.5, 
            label='Pareto Optimal', zorder=5)

ax6.plot(pareto_costs_sorted, pareto_scores_sorted, 'r--', 
         linewidth=2, alpha=0.5, label='Pareto Front')

ax6.set_xlabel('Observation Cost (hours)', fontsize=11, fontweight='bold')
ax6.set_ylabel('Combined Score (Science × Visibility)', fontsize=11, fontweight='bold')
ax6.set_title('Pareto Front Analysis', fontsize=12, fontweight='bold')
ax6.legend(fontsize=9)
ax6.grid(alpha=0.3)

# Annotate top value target
top_target_idx = df_sorted.iloc[0].name
ax6.annotate(f"Best Value:\n{df.iloc[top_target_idx]['Target']}", 
            xy=(df.iloc[top_target_idx]['Cost (hours)'], 
                combined_score[top_target_idx]),
            xytext=(20, -20), textcoords='offset points',
            bbox=dict(boxstyle='round,pad=0.5', fc='lightgreen', alpha=0.8),
            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.3', lw=2),
            fontsize=8, fontweight='bold')

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

# ============================================================================
# Sensitivity Analysis: Effect of alpha parameter
# ============================================================================

print("\n" + "=" * 80)
print("SENSITIVITY ANALYSIS: Effect of Cost Penalty Parameter (α)")
print("=" * 80)

alphas = [0.5, 1.0, 1.5, 2.0]
fig2, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

for idx, alpha_val in enumerate(alphas):
    values_alpha = calculate_value(scientific_importance, visibility_scores, 
                                    observation_costs, alpha_val)
    df_alpha = df.copy()
    df_alpha['Value'] = values_alpha
    df_alpha_sorted = df_alpha.sort_values('Value', ascending=False).reset_index(drop=True)
    
    top10_alpha = df_alpha_sorted.head(10)
    colors_alpha = plt.cm.plasma(np.linspace(0.2, 0.9, 10))
    
    axes[idx].barh(range(10), top10_alpha['Value'], color=colors_alpha, 
                   edgecolor='black', linewidth=1.2)
    axes[idx].set_yticks(range(10))
    axes[idx].set_yticklabels(top10_alpha['Target'], fontsize=9)
    axes[idx].set_xlabel('Value', fontsize=10, fontweight='bold')
    axes[idx].set_title(f'α = {alpha_val} (Top 10 Targets)', 
                       fontsize=11, fontweight='bold')
    axes[idx].invert_yaxis()
    axes[idx].grid(axis='x', alpha=0.3)
    
    print(f"\nα = {alpha_val}:")
    print(f"  Top target: {df_alpha_sorted.iloc[0]['Target']}")
    print(f"  Top value:  {df_alpha_sorted.iloc[0]['Value']:.2f}")

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

print("\n" + "=" * 80)
print("Analysis complete! Figures saved.")
print("=" * 80)

observation_costs = np.random.uniform(0.5, 8.0, n_targets)
scientific_importance = np.random.uniform(30, 100, n_targets)
visibility_scores = np.random.uniform(0.3, 1.0, n_targets)

def calculate_value(science, visibility, cost, alpha=1.0):
    return (science * visibility) / (cost ** alpha)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
import pandas as pd
from typing import List, Tuple, Dict
import seaborn as sns

# Set style for better visualizations
sns.set_style("whitegrid")
plt.rcParams['figure.figsize'] = (14, 10)

class ScientificDataOptimizer:
    """
    Optimizes data compression and transmission strategy for scientific missions
    under bandwidth and battery constraints.
    """
    
    def __init__(self, 
                 bandwidth_budget: float,
                 energy_budget: float,
                 compression_ratio: float = 0.3,
                 compression_quality: float = 0.85,
                 compression_energy_factor: float = 0.1):
        """
        Initialize the optimizer.
        
        Parameters:
        -----------
        bandwidth_budget : float
            Total available bandwidth (MB)
        energy_budget : float
            Total available energy (Joules)
        compression_ratio : float
            Size reduction after compression (0.3 = 70% size reduction)
        compression_quality : float
            Quality retention after compression (0.85 = 85% quality retained)
        compression_energy_factor : float
            Energy cost for compression relative to transmission
        """
        self.bandwidth_budget = bandwidth_budget
        self.energy_budget = energy_budget
        self.compression_ratio = compression_ratio
        self.compression_quality = compression_quality
        self.compression_energy_factor = compression_energy_factor
        
    def generate_sample_data(self, n_samples: int = 50) -> pd.DataFrame:
        """
        Generate realistic scientific data samples with varying characteristics.
        
        Returns:
        --------
        pd.DataFrame with columns: sample_id, data_type, size_mb, scientific_value,
                                    transmission_energy, priority_score
        """
        np.random.seed(42)
        
        # Define different types of scientific data
        data_types = ['Image', 'Spectrometry', 'Temperature', 'Seismic', 'Atmospheric']
        type_characteristics = {
            'Image': {'size_range': (5, 20), 'value_range': (60, 95), 'energy_factor': 1.2},
            'Spectrometry': {'size_range': (2, 8), 'value_range': (70, 100), 'energy_factor': 0.8},
            'Temperature': {'size_range': (0.1, 0.5), 'value_range': (40, 70), 'energy_factor': 0.3},
            'Seismic': {'size_range': (3, 12), 'value_range': (50, 85), 'energy_factor': 1.0},
            'Atmospheric': {'size_range': (1, 5), 'value_range': (55, 80), 'energy_factor': 0.7}
        }
        
        samples = []
        for i in range(n_samples):
            data_type = np.random.choice(data_types)
            char = type_characteristics[data_type]
            
            size = np.random.uniform(*char['size_range'])
            value = np.random.uniform(*char['value_range'])
            # Energy proportional to size with type-specific factor
            energy = size * char['energy_factor'] * np.random.uniform(0.8, 1.2)
            
            # Priority score combines value and urgency
            priority = value * np.random.uniform(0.7, 1.3)
            
            samples.append({
                'sample_id': i,
                'data_type': data_type,
                'size_mb': size,
                'scientific_value': value,
                'transmission_energy': energy,
                'priority_score': priority
            })
        
        return pd.DataFrame(samples)
    
    def optimize_strategy(self, data: pd.DataFrame) -> Dict:
        """
        Optimize data transmission strategy using linear programming.
        
        This solves a multi-objective optimization problem:
        - Maximize: scientific value of transmitted data
        - Subject to: bandwidth and energy constraints
        - Decide for each sample: transmit as-is, compress+transmit, or discard
        
        Returns:
        --------
        Dictionary containing optimization results
        """
        n = len(data)
        
        # Decision variables:
        # x[0:n] = transmit uncompressed (binary)
        # x[n:2n] = transmit compressed (binary)
        # We'll use LP relaxation and round for approximate solution
        
        # Scientific value coefficients (negative because linprog minimizes)
        c = np.concatenate([
            -data['scientific_value'].values,  # uncompressed value
            -data['scientific_value'].values * self.compression_quality  # compressed value
        ])
        
        # Inequality constraints: Ax <= b
        # Constraint 1: Bandwidth
        bandwidth_coef = np.concatenate([
            data['size_mb'].values,  # uncompressed size
            data['size_mb'].values * self.compression_ratio  # compressed size
        ])
        
        # Constraint 2: Energy
        energy_coef = np.concatenate([
            data['transmission_energy'].values,  # transmission energy
            data['transmission_energy'].values * self.compression_ratio + 
            data['size_mb'].values * self.compression_energy_factor  # compressed transmission + compression cost
        ])
        
        # Constraint 3: Each sample can only be sent once (uncompressed OR compressed)
        # For each sample i: x[i] + x[n+i] <= 1
        sample_constraints = np.zeros((n, 2*n))
        for i in range(n):
            sample_constraints[i, i] = 1
            sample_constraints[i, n+i] = 1
        
        # Combine all constraints
        A_ub = np.vstack([
            bandwidth_coef.reshape(1, -1),
            energy_coef.reshape(1, -1),
            sample_constraints
        ])
        
        b_ub = np.concatenate([
            [self.bandwidth_budget],
            [self.energy_budget],
            np.ones(n)
        ])
        
        # Bounds: 0 <= x[i] <= 1 (LP relaxation)
        bounds = [(0, 1) for _ in range(2*n)]
        
        # Solve linear program
        result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')
        
        if not result.success:
            print(f"Warning: Optimization did not converge. Status: {result.message}")
        
        # Extract and round solution
        x_uncompressed = np.round(result.x[:n])
        x_compressed = np.round(result.x[n:])
        
        # Calculate metrics
        transmitted_data = data.copy()
        transmitted_data['transmit_uncompressed'] = x_uncompressed
        transmitted_data['transmit_compressed'] = x_compressed
        transmitted_data['action'] = 'Discard'
        transmitted_data.loc[x_uncompressed == 1, 'action'] = 'Transmit'
        transmitted_data.loc[x_compressed == 1, 'action'] = 'Compress+Transmit'
        
        # Calculate actual values
        transmitted_data['actual_value'] = 0.0
        transmitted_data.loc[x_uncompressed == 1, 'actual_value'] = \
            transmitted_data.loc[x_uncompressed == 1, 'scientific_value']
        transmitted_data.loc[x_compressed == 1, 'actual_value'] = \
            transmitted_data.loc[x_compressed == 1, 'scientific_value'] * self.compression_quality
        
        transmitted_data['actual_size'] = 0.0
        transmitted_data.loc[x_uncompressed == 1, 'actual_size'] = \
            transmitted_data.loc[x_uncompressed == 1, 'size_mb']
        transmitted_data.loc[x_compressed == 1, 'actual_size'] = \
            transmitted_data.loc[x_compressed == 1, 'size_mb'] * self.compression_ratio
        
        transmitted_data['actual_energy'] = 0.0
        transmitted_data.loc[x_uncompressed == 1, 'actual_energy'] = \
            transmitted_data.loc[x_uncompressed == 1, 'transmission_energy']
        transmitted_data.loc[x_compressed == 1, 'actual_energy'] = \
            transmitted_data.loc[x_compressed == 1, 'transmission_energy'] * self.compression_ratio + \
            transmitted_data.loc[x_compressed == 1, 'size_mb'] * self.compression_energy_factor
        
        total_value = transmitted_data['actual_value'].sum()
        total_bandwidth = transmitted_data['actual_size'].sum()
        total_energy = transmitted_data['actual_energy'].sum()
        
        return {
            'data': transmitted_data,
            'total_value': total_value,
            'total_bandwidth_used': total_bandwidth,
            'total_energy_used': total_energy,
            'bandwidth_utilization': total_bandwidth / self.bandwidth_budget * 100,
            'energy_utilization': total_energy / self.energy_budget * 100,
            'samples_transmitted': int(x_uncompressed.sum() + x_compressed.sum()),
            'samples_compressed': int(x_compressed.sum()),
            'samples_uncompressed': int(x_uncompressed.sum()),
            'samples_discarded': n - int(x_uncompressed.sum() + x_compressed.sum())
        }
    
    def compare_strategies(self, data: pd.DataFrame) -> Dict:
        """
        Compare optimized strategy with naive strategies:
        1. Greedy by value (highest value first)
        2. Greedy by size (smallest first)
        3. Compress everything
        4. Optimized strategy
        """
        results = {}
        
        # Strategy 1: Greedy by value
        sorted_by_value = data.sort_values('scientific_value', ascending=False)
        results['greedy_value'] = self._greedy_strategy(sorted_by_value)
        
        # Strategy 2: Greedy by size
        sorted_by_size = data.sort_values('size_mb', ascending=True)
        results['greedy_size'] = self._greedy_strategy(sorted_by_size)
        
        # Strategy 3: Compress everything
        results['compress_all'] = self._compress_all_strategy(data)
        
        # Strategy 4: Optimized
        results['optimized'] = self.optimize_strategy(data)
        
        return results
    
    def _greedy_strategy(self, sorted_data: pd.DataFrame) -> Dict:
        """Greedy strategy: select samples in order until budget exhausted"""
        bandwidth_used = 0
        energy_used = 0
        value_gained = 0
        transmitted = []
        
        for _, row in sorted_data.iterrows():
            if (bandwidth_used + row['size_mb'] <= self.bandwidth_budget and
                energy_used + row['transmission_energy'] <= self.energy_budget):
                bandwidth_used += row['size_mb']
                energy_used += row['transmission_energy']
                value_gained += row['scientific_value']
                transmitted.append(row['sample_id'])
        
        return {
            'total_value': value_gained,
            'total_bandwidth_used': bandwidth_used,
            'total_energy_used': energy_used,
            'samples_transmitted': len(transmitted),
            'transmitted_ids': transmitted
        }
    
    def _compress_all_strategy(self, data: pd.DataFrame) -> Dict:
        """Strategy: compress everything and transmit as much as possible"""
        bandwidth_used = 0
        energy_used = 0
        value_gained = 0
        transmitted = []
        
        sorted_data = data.sort_values('scientific_value', ascending=False)
        
        for _, row in sorted_data.iterrows():
            compressed_size = row['size_mb'] * self.compression_ratio
            compressed_energy = (row['transmission_energy'] * self.compression_ratio + 
                               row['size_mb'] * self.compression_energy_factor)
            
            if (bandwidth_used + compressed_size <= self.bandwidth_budget and
                energy_used + compressed_energy <= self.energy_budget):
                bandwidth_used += compressed_size
                energy_used += compressed_energy
                value_gained += row['scientific_value'] * self.compression_quality
                transmitted.append(row['sample_id'])
        
        return {
            'total_value': value_gained,
            'total_bandwidth_used': bandwidth_used,
            'total_energy_used': energy_used,
            'samples_transmitted': len(transmitted),
            'transmitted_ids': transmitted
        }


def visualize_results(optimizer: ScientificDataOptimizer, 
                     data: pd.DataFrame,
                     comparison: Dict):
    """
    Create comprehensive visualizations of the optimization results.
    """
    fig = plt.figure(figsize=(16, 12))
    
    # 1. Strategy Comparison Bar Chart
    ax1 = plt.subplot(2, 3, 1)
    strategies = ['Greedy\n(Value)', 'Greedy\n(Size)', 'Compress\nAll', 'Optimized']
    values = [comparison['greedy_value']['total_value'],
              comparison['greedy_size']['total_value'],
              comparison['compress_all']['total_value'],
              comparison['optimized']['total_value']]
    
    bars = ax1.bar(strategies, values, color=['#ff9999', '#ffcc99', '#99ccff', '#99ff99'])
    ax1.set_ylabel('Total Scientific Value', fontsize=11, fontweight='bold')
    ax1.set_title('Strategy Comparison: Total Value Achieved', fontsize=12, fontweight='bold')
    ax1.grid(axis='y', alpha=0.3)
    
    # Add value labels on bars
    for bar, val in zip(bars, values):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.1f}', ha='center', va='bottom', fontweight='bold')
    
    # 2. Resource Utilization
    ax2 = plt.subplot(2, 3, 2)
    opt_result = comparison['optimized']
    resources = ['Bandwidth', 'Energy']
    utilization = [opt_result['bandwidth_utilization'], opt_result['energy_utilization']]
    
    bars = ax2.barh(resources, utilization, color=['#3498db', '#e74c3c'])
    ax2.set_xlabel('Utilization (%)', fontsize=11, fontweight='bold')
    ax2.set_title('Optimized Strategy: Resource Utilization', fontsize=12, fontweight='bold')
    ax2.set_xlim(0, 110)
    ax2.grid(axis='x', alpha=0.3)
    
    # Add percentage labels
    for bar, val in zip(bars, utilization):
        width = bar.get_width()
        ax2.text(width + 2, bar.get_y() + bar.get_height()/2.,
                f'{val:.1f}%', ha='left', va='center', fontweight='bold')
    
    # 3. Data Type Distribution in Optimized Strategy
    ax3 = plt.subplot(2, 3, 3)
    opt_data = opt_result['data']
    transmitted_data = opt_data[opt_data['action'] != 'Discard']
    type_counts = transmitted_data['data_type'].value_counts()
    
    colors_pie = ['#ff9999', '#66b3ff', '#99ff99', '#ffcc99', '#ff99cc']
    wedges, texts, autotexts = ax3.pie(type_counts.values, labels=type_counts.index,
                                        autopct='%1.1f%%', colors=colors_pie,
                                        startangle=90)
    ax3.set_title('Transmitted Data by Type', fontsize=12, fontweight='bold')
    
    for autotext in autotexts:
        autotext.set_color('white')
        autotext.set_fontweight('bold')
    
    # 4. Action Distribution
    ax4 = plt.subplot(2, 3, 4)
    action_counts = opt_data['action'].value_counts()
    colors_action = {'Transmit': '#2ecc71', 'Compress+Transmit': '#f39c12', 'Discard': '#e74c3c'}
    
    bars = ax4.bar(range(len(action_counts)), action_counts.values,
                   color=[colors_action[action] for action in action_counts.index])
    ax4.set_xticks(range(len(action_counts)))
    ax4.set_xticklabels(action_counts.index, rotation=0)
    ax4.set_ylabel('Number of Samples', fontsize=11, fontweight='bold')
    ax4.set_title('Decision Distribution: Sample Actions', fontsize=12, fontweight='bold')
    ax4.grid(axis='y', alpha=0.3)
    
    # Add count labels
    for i, (bar, val) in enumerate(zip(bars, action_counts.values)):
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{int(val)}', ha='center', va='bottom', fontweight='bold')
    
    # 5. Value vs Size Scatter Plot
    ax5 = plt.subplot(2, 3, 5)
    colors_scatter = {'Transmit': '#2ecc71', 'Compress+Transmit': '#f39c12', 'Discard': '#e74c3c'}
    
    for action in ['Discard', 'Compress+Transmit', 'Transmit']:
        subset = opt_data[opt_data['action'] == action]
        ax5.scatter(subset['size_mb'], subset['scientific_value'],
                   label=action, alpha=0.6, s=100, color=colors_scatter[action],
                   edgecolors='black', linewidth=0.5)
    
    ax5.set_xlabel('Data Size (MB)', fontsize=11, fontweight='bold')
    ax5.set_ylabel('Scientific Value', fontsize=11, fontweight='bold')
    ax5.set_title('Sample Selection: Value vs Size Trade-off', fontsize=12, fontweight='bold')
    ax5.legend(loc='best', framealpha=0.9)
    ax5.grid(alpha=0.3)
    
    # 6. Efficiency Comparison
    ax6 = plt.subplot(2, 3, 6)
    
    strategies_full = ['Greedy\n(Value)', 'Greedy\n(Size)', 'Compress\nAll', 'Optimized']
    samples_transmitted = [
        comparison['greedy_value']['samples_transmitted'],
        comparison['greedy_size']['samples_transmitted'],
        comparison['compress_all']['samples_transmitted'],
        comparison['optimized']['samples_transmitted']
    ]
    
    x_pos = np.arange(len(strategies_full))
    bars1 = ax6.bar(x_pos - 0.2, values, 0.4, label='Total Value', color='#3498db', alpha=0.8)
    ax6_twin = ax6.twinx()
    bars2 = ax6_twin.bar(x_pos + 0.2, samples_transmitted, 0.4, label='Samples Sent',
                         color='#e67e22', alpha=0.8)
    
    ax6.set_xlabel('Strategy', fontsize=11, fontweight='bold')
    ax6.set_ylabel('Total Scientific Value', fontsize=10, fontweight='bold', color='#3498db')
    ax6_twin.set_ylabel('Samples Transmitted', fontsize=10, fontweight='bold', color='#e67e22')
    ax6.set_title('Efficiency: Value vs. Sample Count', fontsize=12, fontweight='bold')
    ax6.set_xticks(x_pos)
    ax6.set_xticklabels(strategies_full)
    ax6.tick_params(axis='y', labelcolor='#3498db')
    ax6_twin.tick_params(axis='y', labelcolor='#e67e22')
    ax6.grid(axis='y', alpha=0.3)
    
    # Add legend
    lines1, labels1 = ax6.get_legend_handles_labels()
    lines2, labels2 = ax6_twin.get_legend_handles_labels()
    ax6.legend(lines1 + lines2, labels1 + labels2, loc='upper left', framealpha=0.9)
    
    plt.tight_layout()
    plt.savefig('optimization_results.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # Print detailed statistics
    print("\n" + "="*80)
    print("DETAILED OPTIMIZATION RESULTS")
    print("="*80)
    print(f"\nBandwidth Budget: {optimizer.bandwidth_budget} MB")
    print(f"Energy Budget: {optimizer.energy_budget} J")
    print(f"Compression Ratio: {optimizer.compression_ratio} (reduces size to {optimizer.compression_ratio*100:.0f}%)")
    print(f"Compression Quality: {optimizer.compression_quality} (retains {optimizer.compression_quality*100:.0f}% value)")
    
    print("\n" + "-"*80)
    print("STRATEGY COMPARISON")
    print("-"*80)
    print(f"{'Strategy':<20} {'Total Value':>12} {'Samples':>10} {'Bandwidth':>12} {'Energy':>12}")
    print("-"*80)
    
    for name, result in comparison.items():
        strategy_name = name.replace('_', ' ').title()
        print(f"{strategy_name:<20} {result['total_value']:>12.2f} {result['samples_transmitted']:>10} "
              f"{result['total_bandwidth_used']:>11.2f}M {result['total_energy_used']:>11.2f}J")
    
    print("\n" + "-"*80)
    print("OPTIMIZED STRATEGY BREAKDOWN")
    print("-"*80)
    print(f"Samples Transmitted (Uncompressed): {opt_result['samples_uncompressed']}")
    print(f"Samples Transmitted (Compressed):   {opt_result['samples_compressed']}")
    print(f"Samples Discarded:                  {opt_result['samples_discarded']}")
    print(f"Total Samples:                      {len(data)}")
    print(f"\nBandwidth Utilization: {opt_result['bandwidth_utilization']:.2f}%")
    print(f"Energy Utilization:    {opt_result['energy_utilization']:.2f}%")
    
    # Calculate improvement
    best_baseline = max(comparison['greedy_value']['total_value'],
                       comparison['greedy_size']['total_value'],
                       comparison['compress_all']['total_value'])
    improvement = (opt_result['total_value'] - best_baseline) / best_baseline * 100
    
    print(f"\nImprovement over best baseline: {improvement:.2f}%")
    print("="*80)


# Main execution
if __name__ == "__main__":
    print("="*80)
    print("SCIENTIFIC DATA COMPRESSION AND TRANSMISSION OPTIMIZER")
    print("="*80)
    print("\nSimulating a Mars rover mission with limited bandwidth and battery...")
    
    # Initialize optimizer with realistic constraints
    optimizer = ScientificDataOptimizer(
        bandwidth_budget=100.0,      # 100 MB available bandwidth
        energy_budget=250.0,          # 250 Joules available energy
        compression_ratio=0.3,        # Compression reduces size to 30%
        compression_quality=0.85,     # Compression retains 85% of scientific value
        compression_energy_factor=0.1 # Compression costs 0.1 J per MB
    )
    
    # Generate sample scientific data
    print("\nGenerating 50 scientific data samples...")
    data = optimizer.generate_sample_data(n_samples=50)
    
    print("\nSample Data Overview:")
    print(data.head(10).to_string(index=False))
    
    # Run optimization and comparison
    print("\n\nRunning optimization and comparing strategies...")
    comparison_results = optimizer.compare_strategies(data)
    
    # Visualize results
    print("\nGenerating visualizations...")
    visualize_results(optimizer, data, comparison_results)
    
    print("\n✓ Analysis complete! Check the generated plots above.")
    print("\n" + "="*80)
    print("EXECUTION RESULTS WILL APPEAR BELOW")
    print("="*80)
    print("\n\n\n")