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)

print("=" * 70)
print("BATTERY MATERIAL CHARGING CHARACTERISTICS OPTIMIZATION")
print("Trade-off Analysis: Charging Speed vs Lifespan vs Capacity")
print("=" * 70)
print()

# ============================================================================
# OBJECTIVE FUNCTIONS DEFINITION
# ============================================================================

def charging_time(x1, x2):
    """
    Charging time (minutes) - MINIMIZE
    x1: C-rate (charging rate)
    x2: Active material loading (mg/cm^2)
    """
    return (100 / x1) * (1 + 0.02 * x2)

def battery_lifespan(x1, x2):
    """
    Battery lifespan (cycles) - MAXIMIZE (return negative for minimization)
    x1: C-rate (charging rate)
    x2: Active material loading (mg/cm^2)
    """
    return -(1000 - 50 * (x1 ** 1.5) - 10 * ((x2 - 15) ** 2))

def energy_capacity(x1, x2):
    """
    Energy capacity (mAh/cm^2) - MAXIMIZE (return negative for minimization)
    x1: C-rate (charging rate)
    x2: Active material loading (mg/cm^2)
    """
    return -(150 * x2 * np.exp(-0.1 * x1) - 0.5 * (x2 ** 2))

# ============================================================================
# MULTI-OBJECTIVE OPTIMIZATION USING NSGA-II APPROACH
# ============================================================================

class BatteryOptimizer:
    """
    Multi-objective optimizer for battery charging characteristics
    Uses weighted sum approach with multiple weight combinations
    """
    
    def __init__(self):
        self.bounds = [(0.5, 5.0), (5, 25)]  # [C-rate, Loading]
        self.pareto_solutions = []
        
    def weighted_objective(self, x, weights):
        """
        Weighted sum of normalized objectives
        """
        x1, x2 = x
        
        # Calculate objectives
        f1 = charging_time(x1, x2)
        f2 = battery_lifespan(x1, x2)
        f3 = energy_capacity(x1, x2)
        
        # Normalize (approximate ranges)
        f1_norm = f1 / 200.0  # Charging time up to ~200 min
        f2_norm = f2 / 1000.0  # Lifespan range ~1000
        f3_norm = f3 / 2000.0  # Capacity range ~2000
        
        # Weighted sum
        return weights[0] * f1_norm + weights[1] * f2_norm + weights[2] * f3_norm
    
    def optimize_with_weights(self, weights):
        """
        Optimize with specific weight combination
        """
        result = differential_evolution(
            lambda x: self.weighted_objective(x, weights),
            self.bounds,
            maxiter=300,
            popsize=20,
            seed=42,
            atol=1e-6,
            tol=1e-6
        )
        return result.x
    
    def generate_pareto_front(self, n_points=50):
        """
        Generate Pareto front by varying weights
        """
        print("Generating Pareto-optimal solutions...")
        print(f"Running {n_points} optimization scenarios...")
        print()
        
        solutions = []
        
        # Generate diverse weight combinations
        for i in range(n_points):
            # Random weights that sum to 1
            w = np.random.dirichlet([1, 1, 1])
            
            try:
                x_opt = self.optimize_with_weights(w)
                x1, x2 = x_opt
                
                # Calculate actual objective values
                f1 = charging_time(x1, x2)
                f2_actual = -battery_lifespan(x1, x2)  # Convert back to positive
                f3_actual = -energy_capacity(x1, x2)   # Convert back to positive
                
                solutions.append({
                    'x1': x1,
                    'x2': x2,
                    'charging_time': f1,
                    'lifespan': f2_actual,
                    'capacity': f3_actual,
                    'weights': w
                })
                
            except Exception as e:
                continue
        
        self.pareto_solutions = solutions
        print(f"✓ Generated {len(solutions)} Pareto-optimal solutions")
        print()
        
        return solutions

# ============================================================================
# OPTIMIZATION EXECUTION
# ============================================================================

optimizer = BatteryOptimizer()
pareto_solutions = optimizer.generate_pareto_front(n_points=60)

# Extract solution data
x1_vals = [s['x1'] for s in pareto_solutions]
x2_vals = [s['x2'] for s in pareto_solutions]
time_vals = [s['charging_time'] for s in pareto_solutions]
life_vals = [s['lifespan'] for s in pareto_solutions]
capacity_vals = [s['capacity'] for s in pareto_solutions]

# ============================================================================
# IDENTIFY KEY SOLUTIONS
# ============================================================================

print("=" * 70)
print("KEY PARETO-OPTIMAL SOLUTIONS")
print("=" * 70)
print()

# Find extreme solutions
idx_fast = np.argmin(time_vals)
idx_long_life = np.argmax(life_vals)
idx_high_cap = np.argmax(capacity_vals)

# Find balanced solution (closest to center of normalized space)
def normalize(arr):
    return (arr - np.min(arr)) / (np.max(arr) - np.min(arr))

time_norm = normalize(np.array(time_vals))
life_norm = 1 - normalize(np.array(life_vals))  # Invert (want high)
cap_norm = 1 - normalize(np.array(capacity_vals))  # Invert (want high)

distances = np.sqrt(time_norm**2 + life_norm**2 + cap_norm**2)
idx_balanced = np.argmin(distances)

solutions_of_interest = [
    ('Fastest Charging', idx_fast),
    ('Longest Lifespan', idx_long_life),
    ('Highest Capacity', idx_high_cap),
    ('Balanced Solution', idx_balanced)
]

for name, idx in solutions_of_interest:
    sol = pareto_solutions[idx]
    print(f"📊 {name}:")
    print(f"   C-rate (x₁):        {sol['x1']:.3f}")
    print(f"   Loading (x₂):       {sol['x2']:.3f} mg/cm²")
    print(f"   Charging Time:      {sol['charging_time']:.2f} minutes")
    print(f"   Battery Lifespan:   {sol['lifespan']:.0f} cycles")
    print(f"   Energy Capacity:    {sol['capacity']:.2f} mAh/cm²")
    print()

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

print("=" * 70)
print("GENERATING VISUALIZATIONS")
print("=" * 70)
print()

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

# ========== 3D Pareto Front ==========
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
scatter = ax1.scatter(time_vals, life_vals, capacity_vals, 
                     c=capacity_vals, cmap='viridis', s=50, alpha=0.6)
ax1.set_xlabel('Charging Time (min)', fontsize=10, labelpad=10)
ax1.set_ylabel('Battery Lifespan (cycles)', fontsize=10, labelpad=10)
ax1.set_zlabel('Energy Capacity (mAh/cm²)', fontsize=10, labelpad=10)
ax1.set_title('3D Pareto Front: Trade-off Surface', fontsize=12, fontweight='bold', pad=20)
plt.colorbar(scatter, ax=ax1, label='Capacity (mAh/cm²)', shrink=0.5)
ax1.view_init(elev=20, azim=45)

# Mark key solutions
for name, idx in solutions_of_interest:
    sol = pareto_solutions[idx]
    ax1.scatter([sol['charging_time']], [sol['lifespan']], [sol['capacity']], 
               c='red', s=200, marker='*', edgecolors='black', linewidths=2)

# ========== 2D Projections ==========
# Time vs Lifespan
ax2 = fig.add_subplot(2, 3, 2)
scatter2 = ax2.scatter(time_vals, life_vals, c=capacity_vals, 
                       cmap='plasma', s=50, alpha=0.6)
ax2.set_xlabel('Charging Time (min)', fontsize=10)
ax2.set_ylabel('Battery Lifespan (cycles)', fontsize=10)
ax2.set_title('Time vs Lifespan Trade-off', fontsize=11, fontweight='bold')
ax2.grid(True, alpha=0.3)
plt.colorbar(scatter2, ax=ax2, label='Capacity')

for name, idx in solutions_of_interest[:2]:
    sol = pareto_solutions[idx]
    ax2.scatter([sol['charging_time']], [sol['lifespan']], 
               c='red', s=150, marker='*', edgecolors='black', linewidths=1.5)

# Time vs Capacity
ax3 = fig.add_subplot(2, 3, 3)
scatter3 = ax3.scatter(time_vals, capacity_vals, c=life_vals, 
                       cmap='coolwarm', s=50, alpha=0.6)
ax3.set_xlabel('Charging Time (min)', fontsize=10)
ax3.set_ylabel('Energy Capacity (mAh/cm²)', fontsize=10)
ax3.set_title('Time vs Capacity Trade-off', fontsize=11, fontweight='bold')
ax3.grid(True, alpha=0.3)
plt.colorbar(scatter3, ax=ax3, label='Lifespan')

# ========== Design Variable Space ==========
ax4 = fig.add_subplot(2, 3, 4)
scatter4 = ax4.scatter(x1_vals, x2_vals, c=time_vals, 
                       cmap='RdYlGn_r', s=60, alpha=0.7)
ax4.set_xlabel('C-rate (x₁)', fontsize=10)
ax4.set_ylabel('Active Material Loading (x₂) [mg/cm²]', fontsize=10)
ax4.set_title('Design Variable Space (colored by Charging Time)', fontsize=11, fontweight='bold')
ax4.grid(True, alpha=0.3)
plt.colorbar(scatter4, ax=ax4, label='Time (min)')

# ========== 3D Design Space ==========
ax5 = fig.add_subplot(2, 3, 5, projection='3d')
scatter5 = ax5.scatter(x1_vals, x2_vals, time_vals, 
                      c=life_vals, cmap='viridis', s=50, alpha=0.6)
ax5.set_xlabel('C-rate (x₁)', fontsize=10, labelpad=8)
ax5.set_ylabel('Loading (x₂) [mg/cm²]', fontsize=10, labelpad=8)
ax5.set_zlabel('Charging Time (min)', fontsize=10, labelpad=8)
ax5.set_title('Design Space: 3D View', fontsize=11, fontweight='bold', pad=15)
plt.colorbar(scatter5, ax=ax5, label='Lifespan', shrink=0.5)
ax5.view_init(elev=25, azim=135)

# ========== Objective Correlations ==========
ax6 = fig.add_subplot(2, 3, 6)
objectives = np.array([time_vals, life_vals, capacity_vals])
correlation = np.corrcoef(objectives)

im = ax6.imshow(correlation, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
ax6.set_xticks([0, 1, 2])
ax6.set_yticks([0, 1, 2])
ax6.set_xticklabels(['Time', 'Lifespan', 'Capacity'])
ax6.set_yticklabels(['Time', 'Lifespan', 'Capacity'])
ax6.set_title('Objective Function Correlations', fontsize=11, fontweight='bold')

for i in range(3):
    for j in range(3):
        text = ax6.text(j, i, f'{correlation[i, j]:.2f}',
                       ha="center", va="center", color="black", fontsize=11)

plt.colorbar(im, ax=ax6, label='Correlation')

plt.tight_layout()
plt.savefig('battery_optimization_results.png', dpi=150, bbox_inches='tight')
print("✓ Visualization complete")
print()

# ============================================================================
# ADDITIONAL ANALYSIS
# ============================================================================

print("=" * 70)
print("STATISTICAL ANALYSIS OF PARETO SOLUTIONS")
print("=" * 70)
print()

print(f"Charging Time:      {np.min(time_vals):.2f} - {np.max(time_vals):.2f} min")
print(f"Battery Lifespan:   {np.min(life_vals):.0f} - {np.max(life_vals):.0f} cycles")
print(f"Energy Capacity:    {np.min(capacity_vals):.2f} - {np.max(capacity_vals):.2f} mAh/cm²")
print()
print(f"C-rate Range:       {np.min(x1_vals):.3f} - {np.max(x1_vals):.3f}")
print(f"Loading Range:      {np.min(x2_vals):.2f} - {np.max(x2_vals):.2f} mg/cm²")
print()

print("=" * 70)
print("OPTIMIZATION COMPLETE")
print("=" * 70)

plt.show()

w = np.random.dirichlet([1, 1, 1])

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import time

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

# ==================== Rastrigin Function ====================
def rastrigin(x):
    """
    Rastrigin function: f(x) = 10n + sum(x_i^2 - 10*cos(2*pi*x_i))
    Global minimum at x = (0, 0, ..., 0) with f(x) = 0
    """
    n = len(x)
    return 10 * n + np.sum(x**2 - 10 * np.cos(2 * np.pi * x))

# ==================== Genetic Algorithm Implementation ====================
class GeneticAlgorithm:
    def __init__(self, pop_size, dimensions, bounds, mutation_rate, 
                 tournament_size, generations):
        self.pop_size = pop_size
        self.dimensions = dimensions
        self.bounds = bounds
        self.mutation_rate = mutation_rate
        self.tournament_size = tournament_size
        self.generations = generations
        
        # Initialize population randomly within bounds
        self.population = np.random.uniform(
            bounds[0], bounds[1], (pop_size, dimensions)
        )
        
        # History tracking
        self.best_fitness_history = []
        self.avg_fitness_history = []
        self.diversity_history = []
        self.best_solution = None
        self.best_fitness = float('inf')
    
    def evaluate_population(self):
        """Evaluate fitness for entire population"""
        return np.array([rastrigin(ind) for ind in self.population])
    
    def calculate_diversity(self):
        """
        Calculate population diversity as average pairwise distance
        Diversity = (1/n(n-1)) * sum(||x_i - x_j||)
        """
        n = len(self.population)
        total_distance = 0
        count = 0
        
        for i in range(n):
            for j in range(i + 1, n):
                distance = np.linalg.norm(self.population[i] - self.population[j])
                total_distance += distance
                count += 1
        
        return total_distance / count if count > 0 else 0
    
    def tournament_selection(self, fitness):
        """
        Tournament selection: randomly select tournament_size individuals
        and return the best one
        """
        tournament_indices = np.random.choice(
            self.pop_size, self.tournament_size, replace=False
        )
        tournament_fitness = fitness[tournament_indices]
        winner_idx = tournament_indices[np.argmin(tournament_fitness)]
        return self.population[winner_idx].copy()
    
    def crossover(self, parent1, parent2):
        """
        Simulated Binary Crossover (SBX)
        Creates two offspring from two parents
        """
        eta = 20  # Distribution index
        offspring1 = np.zeros(self.dimensions)
        offspring2 = np.zeros(self.dimensions)
        
        for i in range(self.dimensions):
            if np.random.rand() < 0.5:
                if abs(parent1[i] - parent2[i]) > 1e-14:
                    if parent1[i] < parent2[i]:
                        y1, y2 = parent1[i], parent2[i]
                    else:
                        y1, y2 = parent2[i], parent1[i]
                    
                    beta = 1.0 + (2.0 * (y1 - self.bounds[0]) / (y2 - y1))
                    alpha = 2.0 - beta ** -(eta + 1.0)
                    
                    rand = np.random.rand()
                    if rand <= 1.0 / alpha:
                        betaq = (rand * alpha) ** (1.0 / (eta + 1.0))
                    else:
                        betaq = (1.0 / (2.0 - rand * alpha)) ** (1.0 / (eta + 1.0))
                    
                    offspring1[i] = 0.5 * ((y1 + y2) - betaq * (y2 - y1))
                    offspring2[i] = 0.5 * ((y1 + y2) + betaq * (y2 - y1))
                else:
                    offspring1[i] = parent1[i]
                    offspring2[i] = parent2[i]
            else:
                offspring1[i] = parent1[i]
                offspring2[i] = parent2[i]
        
        return offspring1, offspring2
    
    def mutate(self, individual):
        """
        Polynomial mutation
        """
        eta_m = 20  # Distribution index for mutation
        
        for i in range(self.dimensions):
            if np.random.rand() < self.mutation_rate:
                y = individual[i]
                delta1 = (y - self.bounds[0]) / (self.bounds[1] - self.bounds[0])
                delta2 = (self.bounds[1] - y) / (self.bounds[1] - self.bounds[0])
                
                rand = np.random.rand()
                mut_pow = 1.0 / (eta_m + 1.0)
                
                if rand < 0.5:
                    xy = 1.0 - delta1
                    val = 2.0 * rand + (1.0 - 2.0 * rand) * (xy ** (eta_m + 1.0))
                    deltaq = val ** mut_pow - 1.0
                else:
                    xy = 1.0 - delta2
                    val = 2.0 * (1.0 - rand) + 2.0 * (rand - 0.5) * (xy ** (eta_m + 1.0))
                    deltaq = 1.0 - val ** mut_pow
                
                y = y + deltaq * (self.bounds[1] - self.bounds[0])
                individual[i] = np.clip(y, self.bounds[0], self.bounds[1])
        
        return individual
    
    def evolve(self):
        """Main evolution loop"""
        for generation in range(self.generations):
            # Evaluate fitness
            fitness = self.evaluate_population()
            
            # Track statistics
            best_gen_fitness = np.min(fitness)
            avg_gen_fitness = np.mean(fitness)
            diversity = self.calculate_diversity()
            
            self.best_fitness_history.append(best_gen_fitness)
            self.avg_fitness_history.append(avg_gen_fitness)
            self.diversity_history.append(diversity)
            
            # Update best solution
            if best_gen_fitness < self.best_fitness:
                self.best_fitness = best_gen_fitness
                self.best_solution = self.population[np.argmin(fitness)].copy()
            
            # Create new population
            new_population = []
            
            # Elitism: keep the best individual
            best_idx = np.argmin(fitness)
            new_population.append(self.population[best_idx].copy())
            
            # Generate offspring
            while len(new_population) < self.pop_size:
                # Selection
                parent1 = self.tournament_selection(fitness)
                parent2 = self.tournament_selection(fitness)
                
                # Crossover
                offspring1, offspring2 = self.crossover(parent1, parent2)
                
                # Mutation
                offspring1 = self.mutate(offspring1)
                offspring2 = self.mutate(offspring2)
                
                # Add to new population
                new_population.append(offspring1)
                if len(new_population) < self.pop_size:
                    new_population.append(offspring2)
            
            self.population = np.array(new_population[:self.pop_size])
        
        return self.best_solution, self.best_fitness

# ==================== Experiment: Parameter Comparison ====================
def run_experiments():
    """Run GA with different parameter settings"""
    
    dimensions = 2
    bounds = (-5.12, 5.12)
    generations = 100
    
    # Experiment configurations
    experiments = {
        'Small Pop + Low Mutation': {
            'pop_size': 30, 'mutation_rate': 0.01, 'tournament_size': 2
        },
        'Small Pop + High Mutation': {
            'pop_size': 30, 'mutation_rate': 0.2, 'tournament_size': 2
        },
        'Large Pop + Low Mutation': {
            'pop_size': 100, 'mutation_rate': 0.01, 'tournament_size': 2
        },
        'Large Pop + High Mutation': {
            'pop_size': 100, 'mutation_rate': 0.2, 'tournament_size': 2
        },
        'Balanced (Recommended)': {
            'pop_size': 60, 'mutation_rate': 0.1, 'tournament_size': 3
        },
    }
    
    results = {}
    
    print("=" * 70)
    print("GENETIC ALGORITHM PARAMETER OPTIMIZATION EXPERIMENT")
    print("=" * 70)
    print(f"\nObjective: Minimize Rastrigin function in {dimensions}D")
    print(f"Search space: [{bounds[0]}, {bounds[1]}]^{dimensions}")
    print(f"Global optimum: f(0, 0) = 0")
    print(f"\nGenerations: {generations}")
    print("=" * 70)
    
    for name, params in experiments.items():
        print(f"\n[{name}]")
        print(f"  Population: {params['pop_size']}, "
              f"Mutation Rate: {params['mutation_rate']}, "
              f"Tournament Size: {params['tournament_size']}")
        
        start_time = time.time()
        
        ga = GeneticAlgorithm(
            pop_size=params['pop_size'],
            dimensions=dimensions,
            bounds=bounds,
            mutation_rate=params['mutation_rate'],
            tournament_size=params['tournament_size'],
            generations=generations
        )
        
        best_sol, best_fit = ga.evolve()
        elapsed_time = time.time() - start_time
        
        results[name] = {
            'ga': ga,
            'best_solution': best_sol,
            'best_fitness': best_fit,
            'time': elapsed_time
        }
        
        print(f"  Best fitness: {best_fit:.6f}")
        print(f"  Best solution: ({best_sol[0]:.4f}, {best_sol[1]:.4f})")
        print(f"  Time: {elapsed_time:.2f}s")
    
    print("\n" + "=" * 70)
    
    return results

# ==================== Visualization ====================
def create_visualizations(results):
    """Create comprehensive visualizations"""
    
    fig = plt.figure(figsize=(20, 12))
    
    # 1. Fitness convergence comparison
    ax1 = plt.subplot(2, 3, 1)
    for name, data in results.items():
        ga = data['ga']
        ax1.plot(ga.best_fitness_history, label=name, linewidth=2)
    ax1.set_xlabel('Generation', fontsize=11)
    ax1.set_ylabel('Best Fitness', fontsize=11)
    ax1.set_title('Convergence Speed Comparison', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    ax1.set_yscale('log')
    
    # 2. Diversity evolution
    ax2 = plt.subplot(2, 3, 2)
    for name, data in results.items():
        ga = data['ga']
        ax2.plot(ga.diversity_history, label=name, linewidth=2)
    ax2.set_xlabel('Generation', fontsize=11)
    ax2.set_ylabel('Population Diversity', fontsize=11)
    ax2.set_title('Diversity Maintenance', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    # 3. Average fitness
    ax3 = plt.subplot(2, 3, 3)
    for name, data in results.items():
        ga = data['ga']
        ax3.plot(ga.avg_fitness_history, label=name, linewidth=2)
    ax3.set_xlabel('Generation', fontsize=11)
    ax3.set_ylabel('Average Fitness', fontsize=11)
    ax3.set_title('Population Average Fitness', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    ax3.set_yscale('log')
    
    # 4. Final results bar chart
    ax4 = plt.subplot(2, 3, 4)
    names = list(results.keys())
    final_fitness = [results[name]['best_fitness'] for name in names]
    colors = plt.cm.viridis(np.linspace(0, 1, len(names)))
    bars = ax4.bar(range(len(names)), final_fitness, color=colors)
    ax4.set_xticks(range(len(names)))
    ax4.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax4.set_ylabel('Final Best Fitness', fontsize=11)
    ax4.set_title('Final Performance Comparison', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3, axis='y')
    
    # Add value labels on bars
    for i, (bar, val) in enumerate(zip(bars, final_fitness)):
        ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, 
                f'{val:.3f}', ha='center', va='bottom', fontsize=9)
    
    # 5. 3D Rastrigin surface with solutions
    ax5 = plt.subplot(2, 3, 5, projection='3d')
    
    # Create mesh for Rastrigin function
    x = np.linspace(-5.12, 5.12, 100)
    y = np.linspace(-5.12, 5.12, 100)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)
    
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = rastrigin(np.array([X[i, j], Y[i, j]]))
    
    # Plot surface
    surf = ax5.plot_surface(X, Y, Z, cmap=cm.viridis, alpha=0.6, 
                           linewidth=0, antialiased=True)
    
    # Plot best solutions
    colors_3d = plt.cm.Set1(np.linspace(0, 1, len(results)))
    for idx, (name, data) in enumerate(results.items()):
        sol = data['best_solution']
        fit = data['best_fitness']
        ax5.scatter([sol[0]], [sol[1]], [fit], color=colors_3d[idx], 
                   s=100, marker='o', edgecolors='black', linewidths=2,
                   label=name)
    
    ax5.set_xlabel('X₁', fontsize=10)
    ax5.set_ylabel('X₂', fontsize=10)
    ax5.set_zlabel('f(X)', fontsize=10)
    ax5.set_title('Rastrigin Function with Best Solutions', 
                  fontsize=12, fontweight='bold')
    ax5.legend(fontsize=8, loc='upper left')
    
    # 6. Convergence-Diversity trade-off
    ax6 = plt.subplot(2, 3, 6)
    for name, data in results.items():
        ga = data['ga']
        # Plot final 20% of evolution
        start_idx = int(0.8 * len(ga.best_fitness_history))
        ax6.scatter(ga.diversity_history[start_idx:], 
                   ga.best_fitness_history[start_idx:],
                   label=name, alpha=0.6, s=30)
    ax6.set_xlabel('Population Diversity', fontsize=11)
    ax6.set_ylabel('Best Fitness', fontsize=11)
    ax6.set_title('Convergence-Diversity Trade-off (Final 20%)', 
                  fontsize=12, fontweight='bold')
    ax6.legend(fontsize=9)
    ax6.grid(True, alpha=0.3)
    ax6.set_yscale('log')
    
    plt.tight_layout()
    plt.savefig('ga_parameter_optimization.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    # Additional 3D visualization: Population evolution
    fig2 = plt.figure(figsize=(15, 5))
    
    for idx, (name, data) in enumerate(list(results.items())[:3]):
        ax = fig2.add_subplot(1, 3, idx+1, projection='3d')
        ga = data['ga']
        
        # Plot surface
        x = np.linspace(-5.12, 5.12, 50)
        y = np.linspace(-5.12, 5.12, 50)
        X, Y = np.meshgrid(x, y)
        Z = np.zeros_like(X)
        
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                Z[i, j] = rastrigin(np.array([X[i, j], Y[i, j]]))
        
        ax.plot_surface(X, Y, Z, cmap=cm.viridis, alpha=0.3, 
                       linewidth=0, antialiased=True)
        
        # Plot final population
        pop = ga.population
        fitness = np.array([rastrigin(ind) for ind in pop])
        ax.scatter(pop[:, 0], pop[:, 1], fitness, c='red', 
                  marker='o', s=20, alpha=0.6)
        
        # Plot best solution
        sol = data['best_solution']
        fit = data['best_fitness']
        ax.scatter([sol[0]], [sol[1]], [fit], color='gold', 
                  s=200, marker='*', edgecolors='black', linewidths=2)
        
        ax.set_xlabel('X₁', fontsize=10)
        ax.set_ylabel('X₂', fontsize=10)
        ax.set_zlabel('f(X)', fontsize=10)
        ax.set_title(f'{name}\nFinal Population', fontsize=11, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig('ga_final_populations.png', dpi=150, bbox_inches='tight')
    plt.show()

# ==================== Main Execution ====================
if __name__ == "__main__":
    # Run experiments
    results = run_experiments()
    
    # Create visualizations
    print("\n[Generating visualizations...]")
    create_visualizations(results)
    
    print("\n[Analysis Complete]")
    print("\nKey Insights:")
    print("1. Population size affects exploration capability")
    print("2. Mutation rate controls diversity maintenance")
    print("3. Tournament size determines selection pressure")
    print("4. Balance is crucial for optimal performance")

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags, lil_matrix
from scipy.sparse.linalg import spsolve
import time
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

# Solve 2D Laplace equation: d²T/dx² + d²T/dy² = 0
# Boundary conditions: T(0,y)=100, T(1,y)=0, T(x,0)=0, T(x,1)=0

def solve_laplace_2d(nx, ny):
    """
    Solve 2D Laplace equation using finite difference method
    
    Parameters:
    -----------
    nx : int
        Number of grid points in x-direction
    ny : int
        Number of grid points in y-direction
    
    Returns:
    --------
    T : ndarray
        Temperature field (2D array)
    solve_time : float
        Computation time in seconds
    """
    # Grid spacing
    dx = 1.0 / (nx - 1)
    dy = 1.0 / (ny - 1)
    
    # Total number of interior points
    n_interior = (nx - 2) * (ny - 2)
    
    # Initialize sparse matrix (coefficient matrix)
    A = lil_matrix((n_interior, n_interior))
    b = np.zeros(n_interior)
    
    # Build the linear system
    # Interior point index: k = i + j*(nx-2), where i,j are interior indices
    start_time = time.time()
    
    for j in range(ny - 2):
        for i in range(nx - 2):
            k = i + j * (nx - 2)
            
            # Center point coefficient
            A[k, k] = -2.0 / dx**2 - 2.0 / dy**2
            
            # Right neighbor (i+1, j)
            if i < nx - 3:
                A[k, k + 1] = 1.0 / dx**2
            else:
                # Right boundary: T = 0
                b[k] -= 0.0
            
            # Left neighbor (i-1, j)
            if i > 0:
                A[k, k - 1] = 1.0 / dx**2
            else:
                # Left boundary: T = 100
                b[k] -= 100.0 / dx**2
            
            # Top neighbor (i, j+1)
            if j < ny - 3:
                A[k, k + (nx - 2)] = 1.0 / dy**2
            else:
                # Top boundary: T = 0
                b[k] -= 0.0
            
            # Bottom neighbor (i, j-1)
            if j > 0:
                A[k, k - (nx - 2)] = 1.0 / dy**2
            else:
                # Bottom boundary: T = 0
                b[k] -= 0.0
    
    # Convert to CSR format for efficient solving
    A = A.tocsr()
    
    # Solve the linear system
    T_interior = spsolve(A, b)
    solve_time = time.time() - start_time
    
    # Reconstruct full temperature field with boundaries
    T = np.zeros((ny, nx))
    
    # Set boundary conditions
    T[:, 0] = 100.0  # Left boundary
    T[:, -1] = 0.0   # Right boundary
    T[0, :] = 0.0    # Bottom boundary
    T[-1, :] = 0.0   # Top boundary
    
    # Fill interior points
    for j in range(ny - 2):
        for i in range(nx - 2):
            k = i + j * (nx - 2)
            T[j + 1, i + 1] = T_interior[k]
    
    return T, solve_time


def compute_error(T_coarse, T_fine, nx_coarse, nx_fine):
    """
    Compute L2 error between coarse and fine mesh solutions
    
    Parameters:
    -----------
    T_coarse : ndarray
        Temperature field from coarse mesh
    T_fine : ndarray
        Temperature field from fine mesh (reference)
    nx_coarse : int
        Grid points in coarse mesh
    nx_fine : int
        Grid points in fine mesh
    
    Returns:
    --------
    error : float
        L2 relative error
    """
    # Interpolate coarse solution onto fine grid for comparison
    x_coarse = np.linspace(0, 1, nx_coarse)
    y_coarse = np.linspace(0, 1, nx_coarse)
    x_fine = np.linspace(0, 1, nx_fine)
    y_fine = np.linspace(0, 1, nx_fine)
    
    # Simple bilinear interpolation
    from scipy.interpolate import RectBivariateSpline
    interp = RectBivariateSpline(y_coarse, x_coarse, T_coarse)
    T_coarse_interp = interp(y_fine, x_fine)
    
    # Compute L2 relative error
    error = np.sqrt(np.sum((T_coarse_interp - T_fine)**2)) / np.sqrt(np.sum(T_fine**2))
    
    return error


# Mesh resolution study
mesh_sizes = [11, 21, 31, 41, 51, 61, 71, 81, 101, 121]
solve_times = []
errors = []
n_dofs = []  # Degrees of freedom

# Use finest mesh as reference solution
print("Computing reference solution (finest mesh)...")
T_reference, _ = solve_laplace_2d(mesh_sizes[-1], mesh_sizes[-1])
print(f"Reference mesh: {mesh_sizes[-1]}x{mesh_sizes[-1]}")
print()

# Compute solutions for all mesh sizes
print("Mesh Optimization Study:")
print("-" * 70)
print(f"{'Mesh Size':<12} {'DOF':<12} {'Time (s)':<15} {'Relative Error':<15}")
print("-" * 70)

for nx in mesh_sizes:
    ny = nx
    T, t = solve_laplace_2d(nx, ny)
    solve_times.append(t)
    n_dofs.append((nx - 2) * (ny - 2))
    
    # Compute error relative to reference solution
    error = compute_error(T, T_reference, nx, mesh_sizes[-1])
    errors.append(error)
    
    print(f"{nx}x{ny:<7} {n_dofs[-1]:<12} {t:<15.6f} {error:<15.6e}")

print("-" * 70)
print()

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

# 1. Temperature distribution for selected meshes
mesh_indices = [0, 3, 6, -1]  # Show 4 different resolutions
for idx, mesh_idx in enumerate(mesh_indices):
    nx = mesh_sizes[mesh_idx]
    T, _ = solve_laplace_2d(nx, nx)
    
    ax = fig.add_subplot(3, 4, idx + 1)
    x = np.linspace(0, 1, nx)
    y = np.linspace(0, 1, nx)
    X, Y = np.meshgrid(x, y)
    
    contour = ax.contourf(X, Y, T, levels=20, cmap='hot')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(f'Temperature Field: {nx}×{nx} mesh')
    ax.set_aspect('equal')
    plt.colorbar(contour, ax=ax, label='Temperature')

# 2. 3D surface plot of finest mesh solution
ax = fig.add_subplot(3, 4, 5, projection='3d')
nx_plot = mesh_sizes[-1]
T_plot, _ = solve_laplace_2d(nx_plot, nx_plot)
x = np.linspace(0, 1, nx_plot)
y = np.linspace(0, 1, nx_plot)
X, Y = np.meshgrid(x, y)
surf = ax.plot_surface(X, Y, T_plot, cmap='hot', alpha=0.9)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Temperature')
ax.set_title(f'3D Temperature Distribution: {nx_plot}×{nx_plot}')
plt.colorbar(surf, ax=ax, shrink=0.5)

# 3. Computation time vs mesh size
ax = fig.add_subplot(3, 4, 6)
ax.plot(mesh_sizes, solve_times, 'o-', linewidth=2, markersize=8)
ax.set_xlabel('Mesh Size (nx = ny)')
ax.set_ylabel('Computation Time (s)')
ax.set_title('Computational Cost vs Mesh Resolution')
ax.grid(True, alpha=0.3)
ax.set_yscale('log')

# 4. Error vs mesh size
ax = fig.add_subplot(3, 4, 7)
ax.plot(mesh_sizes[:-1], errors[:-1], 's-', linewidth=2, markersize=8, color='red')
ax.set_xlabel('Mesh Size (nx = ny)')
ax.set_ylabel('Relative L2 Error')
ax.set_title('Solution Error vs Mesh Resolution')
ax.grid(True, alpha=0.3)
ax.set_yscale('log')

# 5. Error vs computation time (Pareto front)
ax = fig.add_subplot(3, 4, 8)
ax.plot(solve_times[:-1], errors[:-1], 'd-', linewidth=2, markersize=8, color='green')
for i, nx in enumerate(mesh_sizes[:-1]):
    if i % 2 == 0:  # Label every other point
        ax.annotate(f'{nx}×{nx}', (solve_times[i], errors[i]), 
                   xytext=(5, 5), textcoords='offset points', fontsize=8)
ax.set_xlabel('Computation Time (s)')
ax.set_ylabel('Relative L2 Error')
ax.set_title('Accuracy vs Cost Trade-off (Pareto Front)')
ax.grid(True, alpha=0.3)
ax.set_xscale('log')
ax.set_yscale('log')

# 6. Degrees of freedom vs error
ax = fig.add_subplot(3, 4, 9)
ax.loglog(n_dofs[:-1], errors[:-1], '^-', linewidth=2, markersize=8, color='purple')
# Theoretical O(h²) convergence line
theoretical_slope = -1.0  # L2 error ~ h² ~ 1/N for 2D
x_theory = np.array([n_dofs[2], n_dofs[-2]])
y_theory = errors[2] * (x_theory / n_dofs[2])**theoretical_slope
ax.loglog(x_theory, y_theory, '--', color='gray', label='O(1/N) reference')
ax.set_xlabel('Degrees of Freedom (DOF)')
ax.set_ylabel('Relative L2 Error')
ax.set_title('Convergence Rate Analysis')
ax.legend()
ax.grid(True, alpha=0.3)

# 7. Efficiency metric: Error per unit time
efficiency = np.array(errors[:-1]) * np.array(solve_times[:-1])
ax = fig.add_subplot(3, 4, 10)
ax.plot(mesh_sizes[:-1], efficiency, 'o-', linewidth=2, markersize=8, color='orange')
optimal_idx = np.argmin(efficiency)
ax.plot(mesh_sizes[optimal_idx], efficiency[optimal_idx], 'r*', markersize=20, 
        label=f'Optimal: {mesh_sizes[optimal_idx]}×{mesh_sizes[optimal_idx]}')
ax.set_xlabel('Mesh Size (nx = ny)')
ax.set_ylabel('Error × Time (lower is better)')
ax.set_title('Efficiency Metric')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_yscale('log')

# 8. Speedup relative to finest mesh
speedup = solve_times[-1] / np.array(solve_times[:-1])
ax = fig.add_subplot(3, 4, 11)
ax.plot(mesh_sizes[:-1], speedup, 'v-', linewidth=2, markersize=8, color='brown')
ax.set_xlabel('Mesh Size (nx = ny)')
ax.set_ylabel(f'Speedup (relative to {mesh_sizes[-1]}×{mesh_sizes[-1]})')
ax.set_title('Computational Speedup')
ax.grid(True, alpha=0.3)
ax.set_yscale('log')

# 9. Summary recommendation
ax = fig.add_subplot(3, 4, 12)
ax.axis('off')

# Find optimal mesh based on efficiency
optimal_mesh = mesh_sizes[optimal_idx]
optimal_error = errors[optimal_idx]
optimal_time = solve_times[optimal_idx]
optimal_speedup = speedup[optimal_idx]

summary_text = f"""
MESH OPTIMIZATION SUMMARY

Reference Solution: {mesh_sizes[-1]}×{mesh_sizes[-1]} mesh
Time: {solve_times[-1]:.4f} seconds

Optimal Mesh (Best Efficiency): {optimal_mesh}×{optimal_mesh}
  • Relative Error: {optimal_error:.2e}
  • Computation Time: {optimal_time:.4f} s
  • Speedup: {optimal_speedup:.1f}x faster
  • DOF: {n_dofs[optimal_idx]:,}

Recommendation:
- Use {optimal_mesh}×{optimal_mesh} for good balance
- Use {mesh_sizes[3]}×{mesh_sizes[3]} for quick prototyping
- Use {mesh_sizes[-3]}×{mesh_sizes[-3]} for high accuracy

Trade-off Analysis:
- 2× finer mesh → ~4× more DOF
- ~4× longer computation time
- ~2× better accuracy (2nd order method)
"""

ax.text(0.1, 0.5, summary_text, fontsize=10, verticalalignment='center',
        fontfamily='monospace', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

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

print("\nAnalysis complete! Figure saved as 'mesh_optimization_analysis.png'")

"""
Medical Image Restoration: Optimization-based Denoising and Deblurring
This code demonstrates image restoration using mathematical optimization techniques.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage, optimize
from scipy.signal import convolve2d
import warnings
warnings.filterwarnings('ignore')

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

print("=" * 80)
print("MEDICAL IMAGE RESTORATION: OPTIMIZATION FOR DENOISING AND SHARPENING")
print("=" * 80)

# ============================================================================
# STEP 1: Create Synthetic Medical Image (Brain Scan Simulation)
# ============================================================================
print("\n[STEP 1] Creating synthetic medical image...")

def create_synthetic_brain_image(size=128):
    """Create a synthetic brain-like medical image"""
    img = np.zeros((size, size))
    
    # Brain outline (large circle)
    y, x = np.ogrid[-size//2:size//2, -size//2:size//2]
    brain_mask = x**2 + y**2 <= (size//2.5)**2
    img[brain_mask] = 0.7
    
    # Add internal structures (ventricles, gray matter regions)
    # Left ventricle
    ventricle1 = (x + 15)**2 + (y - 5)**2 <= 150
    img[ventricle1] = 0.3
    
    # Right ventricle
    ventricle2 = (x - 15)**2 + (y - 5)**2 <= 150
    img[ventricle2] = 0.3
    
    # Tumor-like region (bright spot)
    tumor = (x - 20)**2 + (y + 20)**2 <= 80
    img[tumor] = 1.0
    
    # Add some texture
    img += 0.05 * np.random.randn(size, size)
    img = np.clip(img, 0, 1)
    
    return img

original_image = create_synthetic_brain_image(128)
print(f"   - Image size: {original_image.shape}")
print(f"   - Intensity range: [{original_image.min():.3f}, {original_image.max():.3f}]")

# ============================================================================
# STEP 2: Add Degradation (Blur + Noise)
# ============================================================================
print("\n[STEP 2] Applying degradation (blur + noise)...")

# Create Gaussian blur kernel
def create_blur_kernel(size=5, sigma=1.5):
    """Create Gaussian blur kernel"""
    ax = np.arange(-size // 2 + 1., size // 2 + 1.)
    xx, yy = np.meshgrid(ax, ax)
    kernel = np.exp(-(xx**2 + yy**2) / (2. * sigma**2))
    return kernel / np.sum(kernel)

blur_kernel = create_blur_kernel(size=7, sigma=2.0)
blurred_image = convolve2d(original_image, blur_kernel, mode='same', boundary='symm')

# Add Gaussian noise
noise_level = 0.05
noisy_blurred_image = blurred_image + noise_level * np.random.randn(*blurred_image.shape)
noisy_blurred_image = np.clip(noisy_blurred_image, 0, 1)

print(f"   - Blur kernel size: {blur_kernel.shape}")
print(f"   - Noise level (std): {noise_level}")
print(f"   - PSNR (degraded): {-10 * np.log10(np.mean((original_image - noisy_blurred_image)**2)):.2f} dB")

# ============================================================================
# STEP 3: Tikhonov (L2) Regularization - Closed Form Solution
# ============================================================================
print("\n[STEP 3] Applying Tikhonov (L2) regularization...")

def tikhonov_deconvolution(degraded, kernel, lambda_reg=0.01):
    """
    Solve using Fourier domain (frequency space) for fast computation
    Minimizes: ||y - H*x||^2 + lambda * ||x||^2
    """
    # Fourier transform
    degraded_fft = np.fft.fft2(degraded)
    kernel_padded = np.zeros_like(degraded)
    kh, kw = kernel.shape
    kernel_padded[:kh, :kw] = kernel
    kernel_padded = np.roll(kernel_padded, (-kh//2, -kw//2), axis=(0, 1))
    kernel_fft = np.fft.fft2(kernel_padded)
    
    # Wiener filter in frequency domain
    kernel_fft_conj = np.conj(kernel_fft)
    denominator = np.abs(kernel_fft)**2 + lambda_reg
    restored_fft = (kernel_fft_conj * degraded_fft) / denominator
    
    # Inverse Fourier transform
    restored = np.real(np.fft.ifft2(restored_fft))
    return np.clip(restored, 0, 1)

lambda_l2 = 0.001
restored_l2 = tikhonov_deconvolution(noisy_blurred_image, blur_kernel, lambda_l2)
psnr_l2 = -10 * np.log10(np.mean((original_image - restored_l2)**2))
print(f"   - Lambda: {lambda_l2}")
print(f"   - PSNR (restored): {psnr_l2:.2f} dB")

# ============================================================================
# STEP 4: Total Variation (TV) Regularization - Iterative Optimization
# ============================================================================
print("\n[STEP 4] Applying Total Variation (TV) regularization...")
print("   - This preserves edges better than L2 regularization")

def compute_gradient(img):
    """Compute image gradient (finite differences)"""
    grad_x = np.roll(img, -1, axis=1) - img
    grad_y = np.roll(img, -1, axis=0) - img
    return grad_x, grad_y

def compute_divergence(grad_x, grad_y):
    """Compute divergence (adjoint of gradient)"""
    div_x = img - np.roll(img, 1, axis=1)
    div_y = img - np.roll(img, 1, axis=0)
    return div_x + div_y

def tv_deconvolution_gradient_descent(degraded, kernel, lambda_tv=0.01, 
                                       num_iterations=100, step_size=0.01):
    """
    Minimize: ||y - H*x||^2 + lambda * TV(x)
    Using gradient descent with TV regularization
    """
    x = degraded.copy()
    epsilon = 1e-8  # For numerical stability in TV computation
    
    for i in range(num_iterations):
        # Data fidelity gradient: H^T(H*x - y)
        Hx = convolve2d(x, kernel, mode='same', boundary='symm')
        residual = Hx - degraded
        grad_data = convolve2d(residual, kernel[::-1, ::-1], mode='same', boundary='symm')
        
        # TV gradient: -div(∇x / |∇x|)
        grad_x, grad_y = compute_gradient(x)
        magnitude = np.sqrt(grad_x**2 + grad_y**2 + epsilon)
        
        # Normalized gradient
        norm_grad_x = grad_x / magnitude
        norm_grad_y = grad_y / magnitude
        
        # Divergence
        div_x = norm_grad_x - np.roll(norm_grad_x, 1, axis=1)
        div_y = norm_grad_y - np.roll(norm_grad_y, 1, axis=0)
        grad_tv = -(div_x + div_y)
        
        # Gradient descent update
        x = x - step_size * (grad_data + lambda_tv * grad_tv)
        x = np.clip(x, 0, 1)
        
        # Print progress every 20 iterations
        if (i + 1) % 20 == 0:
            psnr_current = -10 * np.log10(np.mean((original_image - x)**2))
            print(f"      Iteration {i+1}/{num_iterations}, PSNR: {psnr_current:.2f} dB")
    
    return x

lambda_tv = 0.05
restored_tv = tv_deconvolution_gradient_descent(noisy_blurred_image, blur_kernel, 
                                                 lambda_tv=lambda_tv, 
                                                 num_iterations=100, 
                                                 step_size=0.005)
psnr_tv = -10 * np.log10(np.mean((original_image - restored_tv)**2))
print(f"   - Final PSNR (TV): {psnr_tv:.2f} dB")

# ============================================================================
# STEP 5: Calculate Quality Metrics
# ============================================================================
print("\n[STEP 5] Computing quality metrics...")

def calculate_metrics(original, restored):
    """Calculate PSNR and SSIM-like metric"""
    mse = np.mean((original - restored)**2)
    psnr = -10 * np.log10(mse) if mse > 0 else float('inf')
    
    # Simple contrast measure
    contrast_orig = np.std(original)
    contrast_rest = np.std(restored)
    
    return psnr, contrast_orig, contrast_rest

psnr_degraded, _, _ = calculate_metrics(original_image, noisy_blurred_image)
psnr_l2_final, contrast_orig, contrast_l2 = calculate_metrics(original_image, restored_l2)
psnr_tv_final, _, contrast_tv = calculate_metrics(original_image, restored_tv)

print(f"   - Degraded Image PSNR: {psnr_degraded:.2f} dB")
print(f"   - L2 Restoration PSNR: {psnr_l2_final:.2f} dB")
print(f"   - TV Restoration PSNR: {psnr_tv_final:.2f} dB")
print(f"   - Original Contrast (std): {contrast_orig:.4f}")
print(f"   - L2 Restored Contrast: {contrast_l2:.4f}")
print(f"   - TV Restored Contrast: {contrast_tv:.4f}")

# ============================================================================
# STEP 6: Visualize Results
# ============================================================================
print("\n[STEP 6] Generating visualization...")

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

# Row 1: Images
images = [
    (original_image, 'Original Image', 'Ground Truth'),
    (noisy_blurred_image, 'Degraded Image', f'Blur + Noise\nPSNR: {psnr_degraded:.1f} dB'),
    (restored_l2, 'L2 Regularization', f'Tikhonov (λ={lambda_l2})\nPSNR: {psnr_l2_final:.1f} dB'),
    (restored_tv, 'TV Regularization', f'Total Variation (λ={lambda_tv})\nPSNR: {psnr_tv_final:.1f} dB')
]

for idx, (img, title, subtitle) in enumerate(images, 1):
    ax = plt.subplot(3, 4, idx)
    im = ax.imshow(img, cmap='gray', vmin=0, vmax=1)
    ax.set_title(f'{title}\n{subtitle}', fontsize=11, fontweight='bold')
    ax.axis('off')
    plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)

# Row 2: Difference maps (error)
plt.subplot(3, 4, 5)
plt.imshow(np.abs(original_image - original_image), cmap='hot', vmin=0, vmax=0.3)
plt.title('Error: Original\n(Reference)', fontsize=10)
plt.colorbar(fraction=0.046, pad=0.04)
plt.axis('off')

plt.subplot(3, 4, 6)
diff_degraded = np.abs(original_image - noisy_blurred_image)
plt.imshow(diff_degraded, cmap='hot', vmin=0, vmax=0.3)
plt.title(f'Error: Degraded\nMSE: {np.mean(diff_degraded**2):.4f}', fontsize=10)
plt.colorbar(fraction=0.046, pad=0.04)
plt.axis('off')

plt.subplot(3, 4, 7)
diff_l2 = np.abs(original_image - restored_l2)
plt.imshow(diff_l2, cmap='hot', vmin=0, vmax=0.3)
plt.title(f'Error: L2\nMSE: {np.mean(diff_l2**2):.4f}', fontsize=10)
plt.colorbar(fraction=0.046, pad=0.04)
plt.axis('off')

plt.subplot(3, 4, 8)
diff_tv = np.abs(original_image - restored_tv)
plt.imshow(diff_tv, cmap='hot', vmin=0, vmax=0.3)
plt.title(f'Error: TV\nMSE: {np.mean(diff_tv**2):.4f}', fontsize=10)
plt.colorbar(fraction=0.046, pad=0.04)
plt.axis('off')

# Row 3: Cross-sections and histograms
row_idx = 64  # Middle row

plt.subplot(3, 4, 9)
plt.plot(original_image[row_idx, :], 'k-', linewidth=2, label='Original')
plt.plot(noisy_blurred_image[row_idx, :], 'r--', linewidth=1.5, alpha=0.7, label='Degraded')
plt.xlabel('Pixel Position', fontsize=10)
plt.ylabel('Intensity', fontsize=10)
plt.title(f'Cross-section (Row {row_idx})', fontsize=11, fontweight='bold')
plt.legend(fontsize=9)
plt.grid(True, alpha=0.3)
plt.ylim([0, 1])

plt.subplot(3, 4, 10)
plt.plot(original_image[row_idx, :], 'k-', linewidth=2, label='Original')
plt.plot(restored_l2[row_idx, :], 'b-', linewidth=1.5, alpha=0.8, label='L2 Restored')
plt.xlabel('Pixel Position', fontsize=10)
plt.ylabel('Intensity', fontsize=10)
plt.title(f'Cross-section: L2', fontsize=11, fontweight='bold')
plt.legend(fontsize=9)
plt.grid(True, alpha=0.3)
plt.ylim([0, 1])

plt.subplot(3, 4, 11)
plt.plot(original_image[row_idx, :], 'k-', linewidth=2, label='Original')
plt.plot(restored_tv[row_idx, :], 'g-', linewidth=1.5, alpha=0.8, label='TV Restored')
plt.xlabel('Pixel Position', fontsize=10)
plt.ylabel('Intensity', fontsize=10)
plt.title(f'Cross-section: TV', fontsize=11, fontweight='bold')
plt.legend(fontsize=9)
plt.grid(True, alpha=0.3)
plt.ylim([0, 1])

# Quality comparison bar chart
plt.subplot(3, 4, 12)
methods = ['Degraded', 'L2', 'TV']
psnr_values = [psnr_degraded, psnr_l2_final, psnr_tv_final]
colors = ['red', 'blue', 'green']
bars = plt.bar(methods, psnr_values, color=colors, alpha=0.7, edgecolor='black')
plt.ylabel('PSNR (dB)', fontsize=11, fontweight='bold')
plt.title('Quality Comparison', fontsize=11, fontweight='bold')
plt.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, psnr_values):
    plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5, 
             f'{val:.1f}', ha='center', va='bottom', fontsize=10, fontweight='bold')

plt.tight_layout()
plt.savefig('medical_image_restoration.png', dpi=150, bbox_inches='tight')
print("   - Visualization saved as 'medical_image_restoration.png'")
plt.show()

# ============================================================================
# Summary and Analysis
# ============================================================================
print("\n" + "=" * 80)
print("SUMMARY AND ANALYSIS")
print("=" * 80)
print("""
OBJECTIVE FUNCTION COMPARISON:

1. Tikhonov (L2) Regularization:
   Minimize: J(x) = 1/2 ||y - Hx||₂² + λ ||x||₂²
   
   - Closed-form solution in Fourier domain (very fast)
   - Tends to over-smooth the image
   - Good for uniform noise removal
   - May blur edges and fine details

2. Total Variation (TV) Regularization:
   Minimize: J(x) = 1/2 ||y - Hx||₂² + λ TV(x)
   where TV(x) = ||∇x||₁
   
   - Requires iterative optimization (slower)
   - Preserves edges while removing noise
   - Better for piecewise-smooth images (medical images!)
   - Creates "staircasing" effect in smooth regions

RESULTS INTERPRETATION:
""")

print(f"├─ Degraded Image: PSNR = {psnr_degraded:.2f} dB (baseline)")
print(f"├─ L2 Restoration: PSNR = {psnr_l2_final:.2f} dB (improvement: {psnr_l2_final - psnr_degraded:.2f} dB)")
print(f"└─ TV Restoration: PSNR = {psnr_tv_final:.2f} dB (improvement: {psnr_tv_final - psnr_degraded:.2f} dB)")

print(f"""
KEY OBSERVATIONS:
- TV regularization achieved {'BETTER' if psnr_tv_final > psnr_l2_final else 'WORSE'} PSNR than L2
- TV preserved edges better (visible in cross-section plots)
- L2 is ~50x faster but produces smoother results
- Choice depends on application: speed vs. edge preservation

MEDICAL IMAGING APPLICATIONS:
✓ Tumor detection: TV better preserves tumor boundaries
✓ Real-time imaging: L2 preferred for speed
✓ High-resolution reconstruction: TV for detail preservation
""")

print("=" * 80)
print("EXECUTION COMPLETED SUCCESSFULLY")
print("=" * 80)