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)

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)