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

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, least_squares
from sklearn.metrics import mean_squared_error, r2_score
import time

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

print("=" * 70)
print("LABORATORY EQUIPMENT CALIBRATION OPTIMIZATION")
print("=" * 70)
print()

# Generate synthetic calibration data
# Simulating a temperature sensor with known systematic errors
n_points = 50
T_measured = np.linspace(20, 100, n_points)  # Measured temperatures (°C)

# True relationship with known coefficients (to simulate real sensor)
true_a0 = -2.5
true_a1 = 1.05
true_a2 = -0.0003
true_a3 = 0.000001

T_true = (true_a0 + true_a1 * T_measured + 
          true_a2 * T_measured**2 + 
          true_a3 * T_measured**3)

# Add measurement noise
noise_std = 0.5
T_true_noisy = T_true + np.random.normal(0, noise_std, n_points)

print("Step 1: Generated Calibration Data")
print(f"  - Number of calibration points: {n_points}")
print(f"  - Temperature range: {T_measured.min():.1f}°C to {T_measured.max():.1f}°C")
print(f"  - Measurement noise (std): {noise_std}°C")
print(f"  - True coefficients: a0={true_a0}, a1={true_a1}, a2={true_a2}, a3={true_a3}")
print()

# Define the calibration model
def calibration_model(coeffs, T_meas):
    """Polynomial calibration model"""
    a0, a1, a2, a3 = coeffs
    return a0 + a1 * T_meas + a2 * T_meas**2 + a3 * T_meas**3

# Define the objective function (residuals for least squares)
def residuals(coeffs, T_meas, T_true):
    """Calculate residuals between model and true values"""
    T_pred = calibration_model(coeffs, T_meas)
    return T_pred - T_true

# Method 1: Using scipy.optimize.least_squares (Recommended for this problem)
print("Step 2: Optimization using least_squares method")
print("-" * 70)
start_time = time.time()

# Initial guess
initial_guess = [0.0, 1.0, 0.0, 0.0]

# Optimize using least_squares (more robust for this type of problem)
result_ls = least_squares(residuals, initial_guess, 
                          args=(T_measured, T_true_noisy),
                          method='lm')  # Levenberg-Marquardt algorithm

time_ls = time.time() - start_time

print(f"  Optimization completed in {time_ls:.4f} seconds")
print(f"  Success: {result_ls.success}")
print(f"  Number of function evaluations: {result_ls.nfev}")
print()

# Extract optimized coefficients
opt_coeffs_ls = result_ls.x
print("  Optimized Calibration Coefficients:")
print(f"    a0 = {opt_coeffs_ls[0]:.6f} (true: {true_a0:.6f})")
print(f"    a1 = {opt_coeffs_ls[1]:.6f} (true: {true_a1:.6f})")
print(f"    a2 = {opt_coeffs_ls[2]:.8f} (true: {true_a2:.8f})")
print(f"    a3 = {opt_coeffs_ls[3]:.10f} (true: {true_a3:.10f})")
print()

# Calculate predictions with optimized coefficients
T_pred_optimized = calibration_model(opt_coeffs_ls, T_measured)

# Calculate error metrics
rmse_before = np.sqrt(mean_squared_error(T_true_noisy, T_measured))
rmse_after = np.sqrt(mean_squared_error(T_true_noisy, T_pred_optimized))
r2_score_val = r2_score(T_true_noisy, T_pred_optimized)

print("Step 3: Performance Metrics")
print("-" * 70)
print(f"  RMSE before calibration: {rmse_before:.4f}°C")
print(f"  RMSE after calibration: {rmse_after:.4f}°C")
print(f"  Improvement: {((rmse_before - rmse_after) / rmse_before * 100):.2f}%")
print(f"  R² Score: {r2_score_val:.6f}")
print()

# Calculate residuals
residuals_before = T_true_noisy - T_measured
residuals_after = T_true_noisy - T_pred_optimized

print("Step 4: Statistical Analysis of Residuals")
print("-" * 70)
print(f"  Before calibration:")
print(f"    Mean error: {np.mean(residuals_before):.4f}°C")
print(f"    Std deviation: {np.std(residuals_before):.4f}°C")
print(f"    Max error: {np.max(np.abs(residuals_before)):.4f}°C")
print()
print(f"  After calibration:")
print(f"    Mean error: {np.mean(residuals_after):.4f}°C")
print(f"    Std deviation: {np.std(residuals_after):.4f}°C")
print(f"    Max error: {np.max(np.abs(residuals_after)):.4f}°C")
print()

# Validation: Test on new data points
print("Step 5: Validation on Independent Test Data")
print("-" * 70)
n_test = 20
T_test_measured = np.linspace(25, 95, n_test)
T_test_true = (true_a0 + true_a1 * T_test_measured + 
               true_a2 * T_test_measured**2 + 
               true_a3 * T_test_measured**3)
T_test_true_noisy = T_test_true + np.random.normal(0, noise_std, n_test)
T_test_pred = calibration_model(opt_coeffs_ls, T_test_measured)

rmse_test = np.sqrt(mean_squared_error(T_test_true_noisy, T_test_pred))
print(f"  Test set RMSE: {rmse_test:.4f}°C")
print(f"  Test set R²: {r2_score(T_test_true_noisy, T_test_pred):.6f}")
print()

# Create comprehensive visualization
fig = plt.figure(figsize=(16, 12))

# Plot 1: Calibration Curve
ax1 = plt.subplot(2, 3, 1)
ax1.scatter(T_measured, T_true_noisy, alpha=0.6, s=50, 
           label='Calibration Data', color='blue')
ax1.plot(T_measured, T_pred_optimized, 'r-', linewidth=2, 
         label='Optimized Model')
ax1.plot(T_measured, T_measured, 'k--', alpha=0.5, 
         label='Perfect Calibration (y=x)')
ax1.set_xlabel('Measured Temperature (°C)', fontsize=11)
ax1.set_ylabel('True Temperature (°C)', fontsize=11)
ax1.set_title('Calibration Curve', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Error Before vs After
ax2 = plt.subplot(2, 3, 2)
ax2.scatter(T_measured, residuals_before, alpha=0.6, s=50,
           label='Before Calibration', color='orange')
ax2.scatter(T_measured, residuals_after, alpha=0.6, s=50,
           label='After Calibration', color='green')
ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
ax2.set_xlabel('Measured Temperature (°C)', fontsize=11)
ax2.set_ylabel('Error (°C)', fontsize=11)
ax2.set_title('Measurement Errors', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Error Distribution
ax3 = plt.subplot(2, 3, 3)
ax3.hist(residuals_before, bins=15, alpha=0.6, label='Before', color='orange')
ax3.hist(residuals_after, bins=15, alpha=0.6, label='After', color='green')
ax3.set_xlabel('Error (°C)', fontsize=11)
ax3.set_ylabel('Frequency', fontsize=11)
ax3.set_title('Error Distribution', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Calibration Function Components
ax4 = plt.subplot(2, 3, 4)
T_range = np.linspace(20, 100, 200)
linear_term = opt_coeffs_ls[0] + opt_coeffs_ls[1] * T_range
quadratic_term = opt_coeffs_ls[2] * T_range**2
cubic_term = opt_coeffs_ls[3] * T_range**3
ax4.plot(T_range, linear_term, label='Linear (a₀ + a₁T)', linewidth=2)
ax4.plot(T_range, quadratic_term, label='Quadratic (a₂T²)', linewidth=2)
ax4.plot(T_range, cubic_term, label='Cubic (a₃T³)', linewidth=2)
ax4.set_xlabel('Measured Temperature (°C)', fontsize=11)
ax4.set_ylabel('Contribution to True Temperature (°C)', fontsize=11)
ax4.set_title('Calibration Function Components', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

# Plot 5: Residual Analysis (Q-Q Plot approximation)
ax5 = plt.subplot(2, 3, 5)
sorted_residuals = np.sort(residuals_after)
theoretical_quantiles = np.linspace(-3, 3, len(sorted_residuals))
ax5.scatter(theoretical_quantiles, sorted_residuals, alpha=0.6, s=40)
ax5.plot(theoretical_quantiles, theoretical_quantiles * np.std(residuals_after), 
         'r--', label='Normal Distribution')
ax5.set_xlabel('Theoretical Quantiles', fontsize=11)
ax5.set_ylabel('Sample Quantiles (°C)', fontsize=11)
ax5.set_title('Q-Q Plot (Normality Check)', fontsize=12, fontweight='bold')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Plot 6: Validation Results
ax6 = plt.subplot(2, 3, 6)
ax6.scatter(T_test_measured, T_test_true_noisy, alpha=0.6, s=50,
           label='Test Data', color='purple')
ax6.plot(T_test_measured, T_test_pred, 'r-', linewidth=2,
         label='Model Prediction')
ax6.plot(T_test_measured, T_test_measured, 'k--', alpha=0.5)
ax6.set_xlabel('Measured Temperature (°C)', fontsize=11)
ax6.set_ylabel('True Temperature (°C)', fontsize=11)
ax6.set_title('Validation on Test Data', fontsize=12, fontweight='bold')
ax6.legend()
ax6.grid(True, alpha=0.3)

plt.suptitle('Laboratory Equipment Calibration Optimization Results', 
             fontsize=14, fontweight='bold', y=0.995)
plt.tight_layout()
plt.show()

print("=" * 70)
print("CALIBRATION OPTIMIZATION COMPLETED SUCCESSFULLY")
print("=" * 70)

# Use sparse matrices if applicable
from scipy.sparse import csr_matrix
# Or parallel processing
from joblib import Parallel, delayed

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

# Circuit parameters
class CircuitParameters:
    def __init__(self):
        self.Vdd = 1.8  # Supply voltage (V)
        self.Cload = 100e-15  # Load capacitance (F) - 100fF
        self.f = 1e9  # Operating frequency (Hz) - 1GHz
        self.alpha = 0.5  # Activity factor
        self.Ileakage_per_um = 1e-9  # Leakage current per um width (A/um)
        self.k = 200e-6  # Transconductance parameter (A/V^2)
        self.L = 0.18e-6  # Channel length (m)
        self.Tamb = 25  # Ambient temperature (°C)
        self.Rthermal = 50  # Thermal resistance (°C/W)
        self.W_min = 0.5e-6  # Minimum width (m)
        self.W_max = 50e-6  # Maximum width (m)
        
params = CircuitParameters()

def calculate_dynamic_power(W, params):
    """Calculate dynamic power consumption"""
    # Gate capacitance proportional to W
    Cgate = 2e-15 * (W / 1e-6)  # fF per um
    Ctotal = params.Cload + Cgate
    P_dyn = params.alpha * Ctotal * params.Vdd**2 * params.f
    return P_dyn

def calculate_leakage_power(W, params):
    """Calculate leakage power consumption"""
    W_um = W / 1e-6
    I_leak = params.Ileakage_per_um * W_um
    P_leak = I_leak * params.Vdd
    return P_leak

def calculate_delay(W, params):
    """Calculate propagation delay"""
    # Drive current proportional to W
    Idrive = params.k * (W / params.L) * (params.Vdd - 0.7)**2
    # Delay inversely proportional to current
    tpd = (params.Cload * params.Vdd) / (2 * Idrive)
    return tpd

def calculate_noise_margin(W, params):
    """Calculate noise margin (simplified model)"""
    # Noise margin improves with square root of W
    W_um = W / 1e-6
    NM = 0.4 * params.Vdd * np.sqrt(W_um / 10)  # Normalized to 10um
    return NM

def calculate_temperature(W, params):
    """Calculate junction temperature"""
    P_total = calculate_dynamic_power(W, params) + calculate_leakage_power(W, params)
    Tjunction = params.Tamb + P_total * params.Rthermal
    return Tjunction

def objective_function(W, params, weights):
    """
    Multi-objective optimization function
    Minimize: weighted sum of normalized metrics
    """
    W = W[0]  # Extract scalar from array
    
    # Calculate metrics
    P_dyn = calculate_dynamic_power(W, params)
    P_leak = calculate_leakage_power(W, params)
    P_total = P_dyn + P_leak
    delay = calculate_delay(W, params)
    NM = calculate_noise_margin(W, params)
    temp = calculate_temperature(W, params)
    
    # Normalization factors (based on typical ranges)
    P_norm = P_total / 1e-3  # Normalize to 1mW
    delay_norm = delay / 1e-9  # Normalize to 1ns
    NM_norm = NM / params.Vdd  # Normalize to Vdd
    temp_norm = (temp - params.Tamb) / 100  # Normalize to 100°C rise
    
    # Objective: minimize power, delay, temp; maximize noise margin
    objective = (weights['power'] * P_norm + 
                 weights['delay'] * delay_norm + 
                 weights['temp'] * temp_norm - 
                 weights['noise'] * NM_norm)
    
    return objective

def optimize_circuit(weights, method='SLSQP'):
    """Optimize circuit with given weights"""
    print(f"\nOptimizing with weights: {weights}")
    print(f"Method: {method}")
    
    # Initial guess: middle of range
    W0 = [(params.W_min + params.W_max) / 2]
    
    # Bounds
    bounds = [(params.W_min, params.W_max)]
    
    start_time = time.time()
    
    if method == 'differential_evolution':
        # Global optimization - slower but more robust
        result = differential_evolution(
            lambda W: objective_function(W, params, weights),
            bounds=bounds,
            maxiter=100,
            seed=42,
            atol=1e-10,
            tol=1e-10
        )
    else:
        # Local optimization - faster
        result = minimize(
            lambda W: objective_function(W, params, weights),
            W0,
            method=method,
            bounds=bounds
        )
    
    elapsed_time = time.time() - start_time
    
    W_opt = result.x[0]
    
    # Calculate all metrics at optimal point
    metrics = {
        'W': W_opt,
        'W_um': W_opt / 1e-6,
        'P_dynamic': calculate_dynamic_power(W_opt, params),
        'P_leakage': calculate_leakage_power(W_opt, params),
        'P_total': calculate_dynamic_power(W_opt, params) + calculate_leakage_power(W_opt, params),
        'delay': calculate_delay(W_opt, params),
        'noise_margin': calculate_noise_margin(W_opt, params),
        'temperature': calculate_temperature(W_opt, params),
        'time': elapsed_time
    }
    
    print(f"Optimal Width: {metrics['W_um']:.2f} μm")
    print(f"Total Power: {metrics['P_total']*1e6:.2f} μW")
    print(f"Delay: {metrics['delay']*1e12:.2f} ps")
    print(f"Noise Margin: {metrics['noise_margin']*1e3:.2f} mV")
    print(f"Temperature: {metrics['temperature']:.2f} °C")
    print(f"Optimization time: {elapsed_time:.4f} seconds")
    
    return metrics

# Define different optimization scenarios
scenarios = {
    'Balanced': {'power': 1.0, 'delay': 1.0, 'noise': 1.0, 'temp': 1.0},
    'Low Power': {'power': 3.0, 'delay': 1.0, 'noise': 0.5, 'temp': 2.0},
    'High Speed': {'power': 0.5, 'delay': 3.0, 'noise': 1.0, 'temp': 0.5},
    'High Reliability': {'power': 1.0, 'delay': 0.5, 'noise': 3.0, 'temp': 1.5},
}

print("="*60)
print("ELECTRIC CIRCUIT DESIGN OPTIMIZATION")
print("="*60)

# Optimize for each scenario
results = {}
for scenario_name, weights in scenarios.items():
    results[scenario_name] = optimize_circuit(weights, method='differential_evolution')

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

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

# 1. Comparison of optimal widths
ax1 = plt.subplot(3, 3, 1)
scenario_names = list(results.keys())
widths = [results[s]['W_um'] for s in scenario_names]
colors = ['#3498db', '#e74c3c', '#2ecc71', '#f39c12']
bars = ax1.bar(scenario_names, widths, color=colors, alpha=0.7, edgecolor='black')
ax1.set_ylabel('Width (μm)', fontsize=11, fontweight='bold')
ax1.set_title('Optimal Transistor Width', fontsize=12, fontweight='bold')
ax1.grid(axis='y', alpha=0.3)
for bar, width in zip(bars, widths):
    ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5, 
             f'{width:.1f}', ha='center', va='bottom', fontweight='bold')

# 2. Power consumption comparison
ax2 = plt.subplot(3, 3, 2)
power_dyn = [results[s]['P_dynamic']*1e6 for s in scenario_names]
power_leak = [results[s]['P_leakage']*1e6 for s in scenario_names]
x = np.arange(len(scenario_names))
width_bar = 0.35
bars1 = ax2.bar(x - width_bar/2, power_dyn, width_bar, label='Dynamic', 
                color='#3498db', alpha=0.7, edgecolor='black')
bars2 = ax2.bar(x + width_bar/2, power_leak, width_bar, label='Leakage', 
                color='#e74c3c', alpha=0.7, edgecolor='black')
ax2.set_ylabel('Power (μW)', fontsize=11, fontweight='bold')
ax2.set_title('Power Consumption', fontsize=12, fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(scenario_names, rotation=15, ha='right')
ax2.legend()
ax2.grid(axis='y', alpha=0.3)

# 3. Delay comparison
ax3 = plt.subplot(3, 3, 3)
delays = [results[s]['delay']*1e12 for s in scenario_names]
bars = ax3.bar(scenario_names, delays, color=colors, alpha=0.7, edgecolor='black')
ax3.set_ylabel('Delay (ps)', fontsize=11, fontweight='bold')
ax3.set_title('Propagation Delay', fontsize=12, fontweight='bold')
ax3.grid(axis='y', alpha=0.3)
for bar, delay in zip(bars, delays):
    ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
             f'{delay:.1f}', ha='center', va='bottom', fontweight='bold')

# 4. Noise margin comparison
ax4 = plt.subplot(3, 3, 4)
noise_margins = [results[s]['noise_margin']*1e3 for s in scenario_names]
bars = ax4.bar(scenario_names, noise_margins, color=colors, alpha=0.7, edgecolor='black')
ax4.set_ylabel('Noise Margin (mV)', fontsize=11, fontweight='bold')
ax4.set_title('Noise Margin', fontsize=12, fontweight='bold')
ax4.grid(axis='y', alpha=0.3)
for bar, nm in zip(bars, noise_margins):
    ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 10, 
             f'{nm:.1f}', ha='center', va='bottom', fontweight='bold')

# 5. Temperature comparison
ax5 = plt.subplot(3, 3, 5)
temps = [results[s]['temperature'] for s in scenario_names]
bars = ax5.bar(scenario_names, temps, color=colors, alpha=0.7, edgecolor='black')
ax5.set_ylabel('Temperature (°C)', fontsize=11, fontweight='bold')
ax5.set_title('Junction Temperature', fontsize=12, fontweight='bold')
ax5.axhline(y=params.Tamb, color='green', linestyle='--', label='Ambient')
ax5.axhline(y=85, color='red', linestyle='--', label='Max Safe')
ax5.legend()
ax5.grid(axis='y', alpha=0.3)

# 6. Trade-off analysis: Power vs Delay
ax6 = plt.subplot(3, 3, 6)
W_range = np.linspace(params.W_min, params.W_max, 100)
power_range = [(calculate_dynamic_power(W, params) + calculate_leakage_power(W, params))*1e6 
               for W in W_range]
delay_range = [calculate_delay(W, params)*1e12 for W in W_range]
ax6.plot(power_range, delay_range, 'b-', linewidth=2, label='Pareto Curve')
for i, scenario in enumerate(scenario_names):
    ax6.scatter(results[scenario]['P_total']*1e6, 
               results[scenario]['delay']*1e12,
               s=150, color=colors[i], marker='o', edgecolor='black', 
               linewidth=2, label=scenario, zorder=5)
ax6.set_xlabel('Total Power (μW)', fontsize=11, fontweight='bold')
ax6.set_ylabel('Delay (ps)', fontsize=11, fontweight='bold')
ax6.set_title('Power-Delay Trade-off', fontsize=12, fontweight='bold')
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)

# 7. Trade-off analysis: Width vs All Metrics (normalized)
ax7 = plt.subplot(3, 3, 7)
W_range_um = W_range / 1e-6
# Normalize all metrics to 0-1 range for comparison
power_norm = np.array(power_range) / max(power_range)
delay_norm = np.array(delay_range) / max(delay_range)
nm_range = [calculate_noise_margin(W, params)*1e3 for W in W_range]
nm_norm = np.array(nm_range) / max(nm_range)
temp_range = [calculate_temperature(W, params) for W in W_range]
temp_norm = (np.array(temp_range) - min(temp_range)) / (max(temp_range) - min(temp_range))

ax7.plot(W_range_um, power_norm, 'r-', linewidth=2, label='Power')
ax7.plot(W_range_um, delay_norm, 'b-', linewidth=2, label='Delay')
ax7.plot(W_range_um, nm_norm, 'g-', linewidth=2, label='Noise Margin')
ax7.plot(W_range_um, temp_norm, 'orange', linewidth=2, label='Temperature')
ax7.set_xlabel('Width (μm)', fontsize=11, fontweight='bold')
ax7.set_ylabel('Normalized Value', fontsize=11, fontweight='bold')
ax7.set_title('Metrics vs Width (Normalized)', fontsize=12, fontweight='bold')
ax7.legend()
ax7.grid(True, alpha=0.3)

# 8. Spider/Radar chart for multi-objective comparison
ax8 = plt.subplot(3, 3, 8, projection='polar')
categories = ['Power\n(lower better)', 'Delay\n(lower better)', 
              'Noise\n(higher better)', 'Temp\n(lower better)']
N = len(categories)

angles = [n / float(N) * 2 * np.pi for n in range(N)]
angles += angles[:1]

for i, scenario in enumerate(scenario_names):
    values = [
        1 - (results[scenario]['P_total']*1e6) / max([results[s]['P_total']*1e6 for s in scenario_names]),
        1 - (results[scenario]['delay']*1e12) / max([results[s]['delay']*1e12 for s in scenario_names]),
        (results[scenario]['noise_margin']*1e3) / max([results[s]['noise_margin']*1e3 for s in scenario_names]),
        1 - (results[scenario]['temperature'] - params.Tamb) / max([results[s]['temperature'] - params.Tamb for s in scenario_names])
    ]
    values += values[:1]
    ax8.plot(angles, values, 'o-', linewidth=2, label=scenario, color=colors[i])
    ax8.fill(angles, values, alpha=0.15, color=colors[i])

ax8.set_xticks(angles[:-1])
ax8.set_xticklabels(categories, size=9)
ax8.set_ylim(0, 1)
ax8.set_title('Multi-Objective Performance\n(Normalized)', 
              fontsize=12, fontweight='bold', pad=20)
ax8.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1))
ax8.grid(True)

# 9. 3D Surface plot: Power and Delay vs Width
ax9 = plt.subplot(3, 3, 9, projection='3d')
W_3d = np.linspace(params.W_min, params.W_max, 50)
W_um_3d = W_3d / 1e-6
power_3d = np.array([(calculate_dynamic_power(W, params) + 
                      calculate_leakage_power(W, params))*1e6 for W in W_3d])
delay_3d = np.array([calculate_delay(W, params)*1e12 for W in W_3d])

ax9.plot(W_um_3d, power_3d, delay_3d, 'b-', linewidth=3, label='Design Space')
for i, scenario in enumerate(scenario_names):
    ax9.scatter(results[scenario]['W_um'],
               results[scenario]['P_total']*1e6,
               results[scenario]['delay']*1e12,
               s=200, color=colors[i], marker='o', 
               edgecolor='black', linewidth=2, label=scenario)

ax9.set_xlabel('Width (μm)', fontsize=10, fontweight='bold')
ax9.set_ylabel('Power (μW)', fontsize=10, fontweight='bold')
ax9.set_zlabel('Delay (ps)', fontsize=10, fontweight='bold')
ax9.set_title('3D Design Space', fontsize=12, fontweight='bold')
ax9.legend(fontsize=8)

plt.tight_layout()
plt.savefig('circuit_optimization_analysis.png', dpi=150, bbox_inches='tight')
print("\nVisualization saved as 'circuit_optimization_analysis.png'")
plt.show()

# Summary table
print("\n" + "="*60)
print("OPTIMIZATION RESULTS SUMMARY")
print("="*60)
print(f"{'Scenario':<20} {'Width(μm)':<12} {'Power(μW)':<12} {'Delay(ps)':<12} {'NM(mV)':<12} {'Temp(°C)':<10}")
print("-"*88)
for scenario in scenario_names:
    r = results[scenario]
    print(f"{scenario:<20} {r['W_um']:>10.2f}  {r['P_total']*1e6:>10.2f}  "
          f"{r['delay']*1e12:>10.2f}  {r['noise_margin']*1e3:>10.2f}  {r['temperature']:>8.2f}")

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, differential_evolution
from scipy.integrate import odeint
import time

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

# ============================================================================
# SECTION 1: Define the Climate Model
# ============================================================================

def climate_model(T, t, params, forcing_func):
    """
    Simplified energy balance climate model
    
    Parameters:
    -----------
    T : float
        Temperature anomaly (°C)
    t : float
        Time (years)
    params : dict
        Model parameters:
        - lambda_climate: climate feedback parameter (W/m²/K)
        - C: heat capacity (J/m²/K)
        - kappa: ocean heat uptake efficiency (W/m²/K)
    forcing_func : function
        Radiative forcing as a function of time
    
    Returns:
    --------
    dT_dt : float
        Rate of temperature change
    """
    lambda_climate = params['lambda_climate']
    C = params['C']
    kappa = params['kappa']
    
    # Radiative forcing at time t
    F = forcing_func(t)
    
    # Energy balance equation (linearized)
    # C * dT/dt = F - lambda * T - kappa * T
    dT_dt = (F - lambda_climate * T - kappa * T) / C
    
    return dT_dt


def forcing_function(t, F0=0, F_rate=0.04):
    """
    Radiative forcing function (linearly increasing)
    Represents CO2 and other greenhouse gas forcing
    
    Parameters:
    -----------
    t : float or array
        Time in years
    F0 : float
        Initial forcing (W/m²)
    F_rate : float
        Rate of forcing increase (W/m²/year)
    """
    return F0 + F_rate * t


def simulate_temperature(params, time_points, forcing_func, T0=0.0):
    """
    Simulate temperature evolution given parameters
    
    Parameters:
    -----------
    params : dict
        Model parameters
    time_points : array
        Time points for simulation
    forcing_func : function
        Forcing function
    T0 : float
        Initial temperature anomaly
    
    Returns:
    --------
    T : array
        Temperature anomaly time series
    """
    T = odeint(climate_model, T0, time_points, args=(params, forcing_func))
    return T.flatten()


# ============================================================================
# SECTION 2: Generate Synthetic Observational Data
# ============================================================================

def generate_observations(time_points, true_params, forcing_func, noise_level=0.05):
    """
    Generate synthetic observational data with noise
    
    Parameters:
    -----------
    time_points : array
        Time points
    true_params : dict
        True parameter values
    forcing_func : function
        Forcing function
    noise_level : float
        Standard deviation of observation noise (°C)
    
    Returns:
    --------
    observations : array
        Synthetic temperature observations
    """
    # Generate true signal
    true_signal = simulate_temperature(true_params, time_points, forcing_func)
    
    # Add observational noise
    noise = np.random.normal(0, noise_level, size=len(time_points))
    observations = true_signal + noise
    
    return observations, true_signal


# ============================================================================
# SECTION 3: Define Objective Function for Optimization
# ============================================================================

def objective_function(param_array, time_points, observations, forcing_func):
    """
    Objective function to minimize: Root Mean Square Error (RMSE)
    
    Parameters:
    -----------
    param_array : array
        Array of parameters [lambda_climate, C, kappa]
    time_points : array
        Time points
    observations : array
        Observed temperature data
    forcing_func : function
        Forcing function
    
    Returns:
    --------
    rmse : float
        Root mean square error
    """
    # Unpack parameters
    lambda_climate, C, kappa = param_array
    
    # Create parameter dictionary
    params = {
        'lambda_climate': lambda_climate,
        'C': C,
        'kappa': kappa
    }
    
    # Simulate temperature
    try:
        T_sim = simulate_temperature(params, time_points, forcing_func)
        
        # Calculate RMSE
        rmse = np.sqrt(np.mean((T_sim - observations)**2))
        
        return rmse
    except:
        # Return large error if simulation fails
        return 1e10


# ============================================================================
# SECTION 4: Optimization Functions
# ============================================================================

def optimize_parameters_local(time_points, observations, forcing_func, 
                              initial_guess, bounds):
    """
    Optimize parameters using local optimization (L-BFGS-B)
    
    Parameters:
    -----------
    time_points : array
        Time points
    observations : array
        Observed data
    forcing_func : function
        Forcing function
    initial_guess : array
        Initial parameter guess
    bounds : list of tuples
        Parameter bounds
    
    Returns:
    --------
    result : OptimizeResult
        Optimization result
    """
    print("Starting local optimization (L-BFGS-B)...")
    start_time = time.time()
    
    result = minimize(
        objective_function,
        initial_guess,
        args=(time_points, observations, forcing_func),
        method='L-BFGS-B',
        bounds=bounds,
        options={'disp': True, 'maxiter': 1000}
    )
    
    elapsed_time = time.time() - start_time
    print(f"Local optimization completed in {elapsed_time:.2f} seconds")
    
    return result


def optimize_parameters_global(time_points, observations, forcing_func, bounds):
    """
    Optimize parameters using global optimization (Differential Evolution)
    
    Parameters:
    -----------
    time_points : array
        Time points
    observations : array
        Observed data
    forcing_func : function
        Forcing function
    bounds : list of tuples
        Parameter bounds
    
    Returns:
    --------
    result : OptimizeResult
        Optimization result
    """
    print("Starting global optimization (Differential Evolution)...")
    start_time = time.time()
    
    result = differential_evolution(
        objective_function,
        bounds,
        args=(time_points, observations, forcing_func),
        strategy='best1bin',
        maxiter=100,
        popsize=15,
        tol=1e-7,
        mutation=(0.5, 1.5),
        recombination=0.7,
        seed=42,
        disp=True,
        workers=1
    )
    
    elapsed_time = time.time() - start_time
    print(f"Global optimization completed in {elapsed_time:.2f} seconds")
    
    return result


# ============================================================================
# SECTION 5: Main Execution and Visualization
# ============================================================================

def main():
    """Main execution function"""
    
    # Define time points (years)
    time_points = np.linspace(0, 100, 101)
    
    # True parameters (to be estimated)
    true_params = {
        'lambda_climate': 1.5,  # W/m²/K (climate feedback)
        'C': 10.0,              # J/m²/K (heat capacity, scaled)
        'kappa': 0.5            # W/m²/K (ocean heat uptake)
    }
    
    print("=" * 70)
    print("CLIMATE MODEL PARAMETER OPTIMIZATION")
    print("=" * 70)
    print("\nTrue Parameters:")
    print(f"  λ (climate feedback): {true_params['lambda_climate']:.3f} W/m²/K")
    print(f"  C (heat capacity):    {true_params['C']:.3f} J/m²/K")
    print(f"  κ (ocean uptake):     {true_params['kappa']:.3f} W/m²/K")
    print()
    
    # Generate synthetic observations
    print("Generating synthetic observational data...")
    observations, true_signal = generate_observations(
        time_points, true_params, forcing_function, noise_level=0.05
    )
    print(f"Generated {len(observations)} observations with noise\n")
    
    # Parameter bounds for optimization
    # [lambda_climate, C, kappa]
    bounds = [
        (0.5, 3.0),   # lambda_climate
        (5.0, 20.0),  # C
        (0.1, 2.0)    # kappa
    ]
    
    # Initial guess (intentionally different from true values)
    initial_guess = np.array([2.0, 15.0, 1.0])
    print("Initial guess:")
    print(f"  λ: {initial_guess[0]:.3f} W/m²/K")
    print(f"  C: {initial_guess[1]:.3f} J/m²/K")
    print(f"  κ: {initial_guess[2]:.3f} W/m²/K")
    print()
    
    # Perform local optimization
    result_local = optimize_parameters_local(
        time_points, observations, forcing_function, initial_guess, bounds
    )
    
    print("\n" + "=" * 70)
    
    # Perform global optimization
    result_global = optimize_parameters_global(
        time_points, observations, forcing_function, bounds
    )
    
    # Extract optimized parameters
    opt_params_local = {
        'lambda_climate': result_local.x[0],
        'C': result_local.x[1],
        'kappa': result_local.x[2]
    }
    
    opt_params_global = {
        'lambda_climate': result_global.x[0],
        'C': result_global.x[1],
        'kappa': result_global.x[2]
    }
    
    # Simulate with optimized parameters
    T_initial = simulate_temperature(
        {'lambda_climate': initial_guess[0], 
         'C': initial_guess[1], 
         'kappa': initial_guess[2]},
        time_points, forcing_function
    )
    T_opt_local = simulate_temperature(opt_params_local, time_points, forcing_function)
    T_opt_global = simulate_temperature(opt_params_global, time_points, forcing_function)
    
    # Calculate errors
    rmse_initial = np.sqrt(np.mean((T_initial - observations)**2))
    rmse_local = np.sqrt(np.mean((T_opt_local - observations)**2))
    rmse_global = np.sqrt(np.mean((T_opt_global - observations)**2))
    
    # Print results
    print("\n" + "=" * 70)
    print("OPTIMIZATION RESULTS")
    print("=" * 70)
    
    print("\nInitial Guess RMSE: {:.6f} °C".format(rmse_initial))
    print("\nLocal Optimization (L-BFGS-B):")
    print(f"  λ: {opt_params_local['lambda_climate']:.6f} W/m²/K (true: {true_params['lambda_climate']:.6f})")
    print(f"  C: {opt_params_local['C']:.6f} J/m²/K (true: {true_params['C']:.6f})")
    print(f"  κ: {opt_params_local['kappa']:.6f} W/m²/K (true: {true_params['kappa']:.6f})")
    print(f"  RMSE: {rmse_local:.6f} °C")
    
    print("\nGlobal Optimization (Differential Evolution):")
    print(f"  λ: {opt_params_global['lambda_climate']:.6f} W/m²/K (true: {true_params['lambda_climate']:.6f})")
    print(f"  C: {opt_params_global['C']:.6f} J/m²/K (true: {true_params['C']:.6f})")
    print(f"  κ: {opt_params_global['kappa']:.6f} W/m²/K (true: {true_params['kappa']:.6f})")
    print(f"  RMSE: {rmse_global:.6f} °C")
    
    # Calculate parameter errors
    print("\nParameter Estimation Errors (Global Optimization):")
    error_lambda = abs(opt_params_global['lambda_climate'] - true_params['lambda_climate'])
    error_C = abs(opt_params_global['C'] - true_params['C'])
    error_kappa = abs(opt_params_global['kappa'] - true_params['kappa'])
    print(f"  Δλ: {error_lambda:.6f} W/m²/K ({100*error_lambda/true_params['lambda_climate']:.2f}%)")
    print(f"  ΔC: {error_C:.6f} J/m²/K ({100*error_C/true_params['C']:.2f}%)")
    print(f"  Δκ: {error_kappa:.6f} W/m²/K ({100*error_kappa/true_params['kappa']:.2f}%)")
    
    # ========================================================================
    # VISUALIZATION
    # ========================================================================
    
    fig = plt.figure(figsize=(16, 12))
    
    # Plot 1: Temperature Time Series
    ax1 = plt.subplot(2, 3, 1)
    ax1.plot(time_points, observations, 'o', alpha=0.5, markersize=4, 
             label='Observations (with noise)', color='black')
    ax1.plot(time_points, true_signal, '-', linewidth=2.5, 
             label='True Signal', color='green')
    ax1.plot(time_points, T_initial, '--', linewidth=2, 
             label='Initial Guess', color='red', alpha=0.7)
    ax1.plot(time_points, T_opt_global, '-', linewidth=2, 
             label='Optimized (Global)', color='blue')
    ax1.set_xlabel('Time (years)', fontsize=11)
    ax1.set_ylabel('Temperature Anomaly (°C)', fontsize=11)
    ax1.set_title('Temperature Evolution: Observations vs Model', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper left', fontsize=9)
    ax1.grid(True, alpha=0.3)
    
    # Plot 2: Residuals
    ax2 = plt.subplot(2, 3, 2)
    residuals_initial = T_initial - observations
    residuals_global = T_opt_global - observations
    ax2.plot(time_points, residuals_initial, 'o-', alpha=0.6, markersize=3,
             label=f'Initial (RMSE={rmse_initial:.4f}°C)', color='red')
    ax2.plot(time_points, residuals_global, 's-', alpha=0.6, markersize=3,
             label=f'Optimized (RMSE={rmse_global:.4f}°C)', color='blue')
    ax2.axhline(y=0, color='black', linestyle='--', linewidth=1)
    ax2.set_xlabel('Time (years)', fontsize=11)
    ax2.set_ylabel('Residual (°C)', fontsize=11)
    ax2.set_title('Model Residuals', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    # Plot 3: Parameter Comparison
    ax3 = plt.subplot(2, 3, 3)
    params_names = ['λ\n(W/m²/K)', 'C\n(J/m²/K)', 'κ\n(W/m²/K)']
    true_vals = [true_params['lambda_climate'], true_params['C'], true_params['kappa']]
    initial_vals = list(initial_guess)
    local_vals = [opt_params_local['lambda_climate'], opt_params_local['C'], opt_params_local['kappa']]
    global_vals = [opt_params_global['lambda_climate'], opt_params_global['C'], opt_params_global['kappa']]
    
    x_pos = np.arange(len(params_names))
    width = 0.2
    
    ax3.bar(x_pos - 1.5*width, true_vals, width, label='True', color='green', alpha=0.8)
    ax3.bar(x_pos - 0.5*width, initial_vals, width, label='Initial', color='red', alpha=0.6)
    ax3.bar(x_pos + 0.5*width, local_vals, width, label='Local Opt', color='orange', alpha=0.8)
    ax3.bar(x_pos + 1.5*width, global_vals, width, label='Global Opt', color='blue', alpha=0.8)
    
    ax3.set_ylabel('Parameter Value', fontsize=11)
    ax3.set_title('Parameter Estimation Comparison', fontsize=12, fontweight='bold')
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels(params_names, fontsize=10)
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3, axis='y')
    
    # Plot 4: RMSE Comparison
    ax4 = plt.subplot(2, 3, 4)
    methods = ['Initial\nGuess', 'Local\nOptimization', 'Global\nOptimization']
    rmse_values = [rmse_initial, rmse_local, rmse_global]
    colors_rmse = ['red', 'orange', 'blue']
    
    bars = ax4.bar(methods, rmse_values, color=colors_rmse, alpha=0.7, edgecolor='black')
    ax4.set_ylabel('RMSE (°C)', fontsize=11)
    ax4.set_title('Root Mean Square Error Comparison', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3, axis='y')
    
    # Add value labels on bars
    for bar, val in zip(bars, rmse_values):
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.5f}', ha='center', va='bottom', fontsize=9)
    
    # Plot 5: Forcing Function
    ax5 = plt.subplot(2, 3, 5)
    forcing_values = forcing_function(time_points)
    ax5.plot(time_points, forcing_values, '-', linewidth=2, color='purple')
    ax5.set_xlabel('Time (years)', fontsize=11)
    ax5.set_ylabel('Radiative Forcing (W/m²)', fontsize=11)
    ax5.set_title('External Radiative Forcing', fontsize=12, fontweight='bold')
    ax5.grid(True, alpha=0.3)
    ax5.fill_between(time_points, 0, forcing_values, alpha=0.3, color='purple')
    
    # Plot 6: Parameter Error Percentages
    ax6 = plt.subplot(2, 3, 6)
    error_percentages = [
        100 * error_lambda / true_params['lambda_climate'],
        100 * error_C / true_params['C'],
        100 * error_kappa / true_params['kappa']
    ]
    
    bars_err = ax6.bar(params_names, error_percentages, color='skyblue', 
                       alpha=0.7, edgecolor='black')
    ax6.set_ylabel('Relative Error (%)', fontsize=11)
    ax6.set_title('Parameter Estimation Error (Global Opt)', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3, axis='y')
    
    # Add value labels
    for bar, val in zip(bars_err, error_percentages):
        height = bar.get_height()
        ax6.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.2f}%', ha='center', va='bottom', fontsize=9)
    
    plt.tight_layout()
    plt.savefig('climate_model_optimization.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n" + "=" * 70)
    print("Visualization saved as 'climate_model_optimization.png'")
    print("=" * 70)


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