import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
from sklearn.cluster import KMeans
from matplotlib.patches import Rectangle
import seaborn as sns

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

# ========================================
# 1. Define Survey Area Parameters
# ========================================
area_size = 10.0  # 10 km x 10 km
n_samples = 15    # Maximum number of drill holes
grid_resolution = 100  # For visualization

# ========================================
# 2. Generate Synthetic Geological Map
# ========================================
print("Generating synthetic geological environment...")

# Create coordinate grids
x_grid = np.linspace(0, area_size, grid_resolution)
y_grid = np.linspace(0, area_size, grid_resolution)
X_mesh, Y_mesh = np.meshgrid(x_grid, y_grid)

# Define multiple geological zones (5 different types)
# Using sinusoidal patterns to create realistic geological boundaries
geology_map = np.zeros((grid_resolution, grid_resolution))

# Zone 1: Sedimentary basin (northwest)
zone1 = ((X_mesh < 5) & (Y_mesh > 5)).astype(float)

# Zone 2: Igneous intrusion (center)
zone2 = (np.sqrt((X_mesh - 5)**2 + (Y_mesh - 5)**2) < 2.5).astype(float) * 2

# Zone 3: Metamorphic belt (diagonal)
zone3 = (np.abs(X_mesh - Y_mesh) < 1.5).astype(float) * 3

# Zone 4: Volcanic (southeast)
zone4 = ((X_mesh > 6) & (Y_mesh < 4)).astype(float) * 4

# Zone 5: Alluvial deposits (remaining areas)
geology_map = zone1 + zone2 + zone3 + zone4

# ========================================
# 3. Generate Chemical Gradient
# ========================================
# Simulate metal concentration gradient (e.g., copper)
# Higher in center, decreasing outward with some noise
distance_from_center = np.sqrt((X_mesh - 5)**2 + (Y_mesh - 5)**2)
chemical_gradient = 100 * np.exp(-distance_from_center / 3.0)
chemical_gradient += np.random.randn(grid_resolution, grid_resolution) * 5

# Simulate pH gradient (east-west trend)
ph_gradient = 6.5 + 0.3 * X_mesh + np.random.randn(grid_resolution, grid_resolution) * 0.2

# ========================================
# 4. Latin Hypercube Sampling (LHS)
# ========================================
print(f"Generating {n_samples} sampling locations using Latin Hypercube Sampling...")

def latin_hypercube_sampling(n_samples, dimensions, bounds):
    """
    Generate Latin Hypercube Samples
    
    Parameters:
    - n_samples: number of samples
    - dimensions: number of dimensions (2 for x, y)
    - bounds: [(min, max), ...] for each dimension
    """
    samples = np.zeros((n_samples, dimensions))
    
    for d in range(dimensions):
        # Divide range into n_samples intervals
        intervals = np.linspace(bounds[d][0], bounds[d][1], n_samples + 1)
        
        # Randomly sample within each interval
        for i in range(n_samples):
            samples[i, d] = np.random.uniform(intervals[i], intervals[i + 1])
    
    # Shuffle each dimension independently
    for d in range(dimensions):
        np.random.shuffle(samples[:, d])
    
    return samples

# Generate initial LHS samples
lhs_samples = latin_hypercube_sampling(n_samples, 2, [(0, area_size), (0, area_size)])

# ========================================
# 5. Extract Features at Sample Locations
# ========================================
def get_features_at_location(x, y):
    """Extract geological and chemical features at given location"""
    # Convert to grid indices
    ix = int(np.clip(x / area_size * grid_resolution, 0, grid_resolution - 1))
    iy = int(np.clip(y / area_size * grid_resolution, 0, grid_resolution - 1))
    
    return {
        'geology': geology_map[iy, ix],
        'metal_conc': chemical_gradient[iy, ix],
        'ph': ph_gradient[iy, ix]
    }

# Extract features for all samples
sample_features = []
for i in range(n_samples):
    features = get_features_at_location(lhs_samples[i, 0], lhs_samples[i, 1])
    sample_features.append([
        lhs_samples[i, 0], 
        lhs_samples[i, 1],
        features['geology'],
        features['metal_conc'],
        features['ph']
    ])

sample_features = np.array(sample_features)

# ========================================
# 6. Diversity Metrics Calculation
# ========================================
print("Calculating diversity metrics...")

# Geological diversity: number of unique geology types covered
unique_geologies = len(np.unique(sample_features[:, 2]))
geology_diversity = unique_geologies / 5.0  # Normalized (5 total zones)

# Chemical diversity: coverage of concentration range
conc_range = np.max(sample_features[:, 3]) - np.min(sample_features[:, 3])
full_conc_range = np.max(chemical_gradient) - np.min(chemical_gradient)
chemical_diversity = conc_range / full_conc_range

# pH diversity
ph_range = np.max(sample_features[:, 4]) - np.min(sample_features[:, 4])
full_ph_range = np.max(ph_gradient) - np.min(ph_gradient)
ph_diversity = ph_range / full_ph_range

# Spatial diversity: measure how well-distributed samples are
# Using minimum spanning tree concept
spatial_distances = cdist(lhs_samples, lhs_samples)
np.fill_diagonal(spatial_distances, np.inf)
min_distances = np.min(spatial_distances, axis=1)
spatial_diversity = np.mean(min_distances) / (area_size / np.sqrt(n_samples))

print(f"\nDiversity Metrics:")
print(f"  Geological Diversity: {geology_diversity:.2%}")
print(f"  Chemical Diversity: {chemical_diversity:.2%}")
print(f"  pH Diversity: {ph_diversity:.2%}")
print(f"  Spatial Diversity Index: {spatial_diversity:.2f}")

# ========================================
# 7. K-means Clustering for Zone Analysis
# ========================================
print("\nPerforming k-means clustering on features...")

# Normalize features for clustering
feature_matrix = sample_features[:, 2:5]
feature_matrix_normalized = (feature_matrix - feature_matrix.mean(axis=0)) / feature_matrix.std(axis=0)

# Apply k-means to group similar samples
n_clusters = 5
kmeans = KMeans(n_clusters=n_clusters, random_state=42)
clusters = kmeans.fit_predict(feature_matrix_normalized)

# ========================================
# 8. Visualization
# ========================================
print("\nGenerating visualizations...")

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

# Plot 1: Geological Map with Sample Locations
ax1 = plt.subplot(2, 3, 1)
im1 = ax1.contourf(X_mesh, Y_mesh, geology_map, levels=20, cmap='terrain')
scatter1 = ax1.scatter(lhs_samples[:, 0], lhs_samples[:, 1], 
                       c='red', s=200, marker='*', 
                       edgecolors='black', linewidths=2, 
                       label='Drill Locations', zorder=5)
ax1.set_xlabel('X (km)', fontsize=12)
ax1.set_ylabel('Y (km)', fontsize=12)
ax1.set_title('Geological Map with Sampling Locations', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
plt.colorbar(im1, ax=ax1, label='Geology Type')

# Plot 2: Metal Concentration Map
ax2 = plt.subplot(2, 3, 2)
im2 = ax2.contourf(X_mesh, Y_mesh, chemical_gradient, levels=20, cmap='YlOrRd')
ax2.scatter(lhs_samples[:, 0], lhs_samples[:, 1], 
           c='blue', s=200, marker='*', 
           edgecolors='white', linewidths=2, zorder=5)
ax2.set_xlabel('X (km)', fontsize=12)
ax2.set_ylabel('Y (km)', fontsize=12)
ax2.set_title('Metal Concentration (ppm)', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
plt.colorbar(im2, ax=ax2, label='Concentration (ppm)')

# Plot 3: pH Gradient Map
ax3 = plt.subplot(2, 3, 3)
im3 = ax3.contourf(X_mesh, Y_mesh, ph_gradient, levels=20, cmap='RdYlBu_r')
ax3.scatter(lhs_samples[:, 0], lhs_samples[:, 1], 
           c='green', s=200, marker='*', 
           edgecolors='black', linewidths=2, zorder=5)
ax3.set_xlabel('X (km)', fontsize=12)
ax3.set_ylabel('Y (km)', fontsize=12)
ax3.set_title('pH Distribution', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3)
plt.colorbar(im3, ax=ax3, label='pH')

# Plot 4: Sample Clustering
ax4 = plt.subplot(2, 3, 4)
scatter4 = ax4.scatter(lhs_samples[:, 0], lhs_samples[:, 1], 
                       c=clusters, s=300, marker='o', 
                       cmap='Set3', edgecolors='black', linewidths=2)
for i in range(n_samples):
    ax4.annotate(f'{i+1}', (lhs_samples[i, 0], lhs_samples[i, 1]), 
                ha='center', va='center', fontweight='bold')
ax4.set_xlabel('X (km)', fontsize=12)
ax4.set_ylabel('Y (km)', fontsize=12)
ax4.set_title('Sample Clusters (Similar Characteristics)', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)
plt.colorbar(scatter4, ax=ax4, label='Cluster ID')

# Plot 5: Feature Distribution
ax5 = plt.subplot(2, 3, 5)
feature_names = ['Geology Type', 'Metal Conc.', 'pH']
x_pos = np.arange(len(feature_names))

# Normalize all features to 0-1 scale for comparison
geology_norm = sample_features[:, 2] / 4.0  # Max geology type is 4
metal_norm = (sample_features[:, 3] - sample_features[:, 3].min()) / (sample_features[:, 3].max() - sample_features[:, 3].min())
ph_norm = (sample_features[:, 4] - sample_features[:, 4].min()) / (sample_features[:, 4].max() - sample_features[:, 4].min())

bp = ax5.boxplot([geology_norm, metal_norm, ph_norm], 
                  labels=feature_names, patch_artist=True)
for patch, color in zip(bp['boxes'], ['lightblue', 'lightcoral', 'lightgreen']):
    patch.set_facecolor(color)
ax5.set_ylabel('Normalized Value (0-1)', fontsize=12)
ax5.set_title('Feature Distribution Across Samples', fontsize=14, fontweight='bold')
ax5.grid(True, alpha=0.3, axis='y')

# Plot 6: Diversity Metrics Bar Chart
ax6 = plt.subplot(2, 3, 6)
metrics = ['Geological\nDiversity', 'Chemical\nDiversity', 'pH\nDiversity', 'Spatial\nDiversity']
values = [geology_diversity, chemical_diversity, ph_diversity, min(spatial_diversity, 1.0)]
colors_bar = ['#2ecc71', '#3498db', '#e74c3c', '#f39c12']

bars = ax6.bar(metrics, values, color=colors_bar, edgecolor='black', linewidth=2)
ax6.set_ylabel('Diversity Score', fontsize=12)
ax6.set_title('Optimization Performance Metrics', fontsize=14, fontweight='bold')
ax6.set_ylim([0, 1.2])
ax6.axhline(y=1.0, color='gray', linestyle='--', linewidth=1, label='Target')
ax6.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar, value in zip(bars, values):
    height = bar.get_height()
    ax6.text(bar.get_x() + bar.get_width()/2., height,
            f'{value:.2%}',
            ha='center', va='bottom', fontweight='bold', fontsize=11)

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

# ========================================
# 9. Generate Summary Report
# ========================================
print("\n" + "="*60)
print("SAMPLING STRATEGY OPTIMIZATION REPORT")
print("="*60)
print(f"\nSurvey Area: {area_size} km × {area_size} km")
print(f"Total Samples: {n_samples} drill holes")
print(f"Sampling Method: Latin Hypercube Sampling (LHS)")
print(f"\nCoverage Summary:")
print(f"  - Geological zones covered: {unique_geologies}/5 ({geology_diversity:.1%})")
print(f"  - Metal concentration range: {conc_range:.1f} ppm ({chemical_diversity:.1%} of total)")
print(f"  - pH range covered: {ph_range:.2f} units ({ph_diversity:.1%} of total)")
print(f"  - Average minimum distance between samples: {np.mean(min_distances):.2f} km")
print(f"\nCluster Analysis:")
for cluster_id in range(n_clusters):
    cluster_samples = np.where(clusters == cluster_id)[0] + 1
    print(f"  Cluster {cluster_id + 1}: Samples {list(cluster_samples)}")

print("\n" + "="*60)
print("RECOMMENDATIONS:")
print("="*60)
print("1. The sampling strategy achieves excellent spatial coverage")
print("2. All major geological zones are represented")
print("3. Chemical gradients are well-captured across the survey area")
print("4. Consider prioritizing samples in Clusters 1-2 for initial analysis")
print("5. If budget allows, add 3-5 infill samples in undersampled zones")
print("="*60)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution
from scipy.stats import norm
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, ConstantKernel as C
import warnings
warnings.filterwarnings('ignore')

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

# ============================================================================
# TRUE RESPONSE FUNCTION (Unknown in real experiments)
# ============================================================================
def true_response_function(params):
    """
    Simulates the true signal intensity as a function of instrument parameters.
    
    Parameters:
    - params[0]: Collision Energy (CE) [10-50 eV]
    - params[1]: Capillary Voltage [2000-4000 V]
    - params[2]: Mobile Phase pH [2.0-7.0]
    - params[3]: Column Temperature [25-50 °C]
    
    Returns:
    - Signal intensity (arbitrary units)
    """
    ce, voltage, ph, temp = params
    
    # Optimal conditions (unknown to optimizer)
    ce_opt = 28.5
    voltage_opt = 3200
    ph_opt = 4.2
    temp_opt = 38.0
    
    # Response surface with different sensitivities for each parameter
    signal = 100000 * np.exp(
        -0.008 * (ce - ce_opt)**2
        -0.0000015 * (voltage - voltage_opt)**2
        -0.25 * (ph - ph_opt)**2
        -0.015 * (temp - temp_opt)**2
    )
    
    # Add interaction effect between pH and voltage
    signal *= (1 + 0.15 * np.sin((ph - 4) * (voltage - 3000) / 500))
    
    # Add realistic noise (±5%)
    noise = np.random.normal(0, 0.05 * signal)
    
    return signal + noise

# ============================================================================
# PARAMETER BOUNDS
# ============================================================================
bounds = np.array([
    [10, 50],      # Collision Energy (eV)
    [2000, 4000],  # Capillary Voltage (V)
    [2.0, 7.0],    # pH
    [25, 50]       # Temperature (°C)
])

param_names = ['Collision Energy (eV)', 'Capillary Voltage (V)', 
               'Mobile Phase pH', 'Column Temperature (°C)']

# ============================================================================
# BAYESIAN OPTIMIZATION IMPLEMENTATION
# ============================================================================
class BayesianOptimizer:
    def __init__(self, bounds, n_initial=10):
        self.bounds = bounds
        self.n_params = len(bounds)
        self.X_observed = []
        self.y_observed = []
        
        # Gaussian Process with Matern kernel
        kernel = C(1.0, (1e-3, 1e3)) * Matern(length_scale=1.0, nu=2.5)
        self.gp = GaussianProcessRegressor(
            kernel=kernel,
            n_restarts_optimizer=10,
            alpha=1e-6,
            normalize_y=True
        )
        
        # Initial random sampling
        self.initial_sampling(n_initial)
    
    def initial_sampling(self, n):
        """Perform initial random sampling"""
        print("=== Initial Random Sampling ===")
        for i in range(n):
            # Latin Hypercube Sampling for better coverage
            params = np.array([
                np.random.uniform(b[0], b[1]) for b in self.bounds
            ])
            signal = true_response_function(params)
            
            self.X_observed.append(params)
            self.y_observed.append(signal)
            
            print(f"Experiment {i+1}: Signal = {signal:.1f}")
    
    def expected_improvement(self, X, xi=0.01):
        """
        Calculate Expected Improvement acquisition function
        
        EI(x) = E[max(f(x) - f(x+), 0)]
        """
        X = np.atleast_2d(X)
        
        mu, sigma = self.gp.predict(X, return_std=True)
        mu = mu.reshape(-1, 1)
        sigma = sigma.reshape(-1, 1)
        
        y_max = np.max(self.y_observed)
        
        with np.errstate(divide='warn'):
            imp = mu - y_max - xi
            Z = imp / sigma
            ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
            ei[sigma == 0.0] = 0.0
        
        return -ei.flatten()  # Negative for minimization
    
    def propose_location(self):
        """Find the point with maximum Expected Improvement"""
        # Normalize bounds for optimization
        bounds_norm = [(0, 1)] * self.n_params
        
        # Use differential evolution for global optimization
        result = differential_evolution(
            lambda x: self.expected_improvement(self.denormalize(x)),
            bounds_norm,
            maxiter=100,
            seed=42
        )
        
        return self.denormalize(result.x)
    
    def denormalize(self, x_norm):
        """Convert normalized [0,1] to actual parameter ranges"""
        return np.array([
            x_norm[i] * (self.bounds[i][1] - self.bounds[i][0]) + self.bounds[i][0]
            for i in range(self.n_params)
        ])
    
    def optimize(self, n_iter=20):
        """Run Bayesian Optimization"""
        print("\n=== Bayesian Optimization ===")
        
        for iteration in range(n_iter):
            # Fit GP to current data
            X_array = np.array(self.X_observed)
            y_array = np.array(self.y_observed)
            self.gp.fit(X_array, y_array)
            
            # Propose next point
            next_params = self.propose_location()
            next_signal = true_response_function(next_params)
            
            # Update observations
            self.X_observed.append(next_params)
            self.y_observed.append(next_signal)
            
            best_signal = np.max(self.y_observed)
            print(f"Iteration {iteration+1}: Signal = {next_signal:.1f}, Best = {best_signal:.1f}")
        
        # Final results
        best_idx = np.argmax(self.y_observed)
        best_params = self.X_observed[best_idx]
        best_signal = self.y_observed[best_idx]
        
        return best_params, best_signal

# ============================================================================
# RUN OPTIMIZATION
# ============================================================================
print("=" * 60)
print("LC-MS PARAMETER OPTIMIZATION")
print("=" * 60)

optimizer = BayesianOptimizer(bounds, n_initial=10)
best_params, best_signal = optimizer.optimize(n_iter=20)

print("\n" + "=" * 60)
print("OPTIMIZATION RESULTS")
print("=" * 60)
for i, name in enumerate(param_names):
    print(f"{name}: {best_params[i]:.2f}")
print(f"\nMaximum Signal Intensity: {best_signal:.1f}")
print("=" * 60)

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

# Plot 1: Optimization Progress
ax1 = plt.subplot(3, 3, 1)
iterations = np.arange(1, len(optimizer.y_observed) + 1)
cumulative_best = np.maximum.accumulate(optimizer.y_observed)

ax1.plot(iterations, optimizer.y_observed, 'o-', alpha=0.6, label='Observed Signal')
ax1.plot(iterations, cumulative_best, 'r-', linewidth=2, label='Best Signal')
ax1.axvline(x=10, color='gray', linestyle='--', alpha=0.5, label='Random → Bayesian')
ax1.set_xlabel('Experiment Number', fontsize=11)
ax1.set_ylabel('Signal Intensity', fontsize=11)
ax1.set_title('Optimization Progress', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2-5: Response surfaces for each parameter
X_array = np.array(optimizer.X_observed)
y_array = np.array(optimizer.y_observed)

for idx in range(4):
    ax = plt.subplot(3, 3, idx + 2)
    
    # Create grid for this parameter
    param_range = np.linspace(bounds[idx][0], bounds[idx][1], 100)
    signals = []
    
    for p in param_range:
        test_params = best_params.copy()
        test_params[idx] = p
        signals.append(true_response_function(test_params))
    
    ax.plot(param_range, signals, 'b-', linewidth=2, label='Response Surface')
    ax.scatter(X_array[:, idx], y_array, alpha=0.5, s=50, c='orange', 
               edgecolors='black', label='Experiments')
    ax.axvline(x=best_params[idx], color='red', linestyle='--', 
               linewidth=2, label='Optimum')
    
    ax.set_xlabel(param_names[idx], fontsize=10)
    ax.set_ylabel('Signal Intensity', fontsize=10)
    ax.set_title(f'Effect of {param_names[idx]}', fontsize=11, fontweight='bold')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

# Plot 6: 2D contour (pH vs Temperature)
ax6 = plt.subplot(3, 3, 6)
ph_range = np.linspace(bounds[2][0], bounds[2][1], 50)
temp_range = np.linspace(bounds[3][0], bounds[3][1], 50)
PH, TEMP = np.meshgrid(ph_range, temp_range)

Z = np.zeros_like(PH)
for i in range(len(ph_range)):
    for j in range(len(temp_range)):
        test_params = best_params.copy()
        test_params[2] = ph_range[i]
        test_params[3] = temp_range[j]
        Z[j, i] = true_response_function(test_params)

contour = ax6.contourf(PH, TEMP, Z, levels=20, cmap='viridis')
ax6.scatter(X_array[:, 2], X_array[:, 3], c='white', s=50, 
            edgecolors='black', alpha=0.7)
ax6.plot(best_params[2], best_params[3], 'r*', markersize=20, 
         label='Optimum')
ax6.set_xlabel('Mobile Phase pH', fontsize=10)
ax6.set_ylabel('Column Temperature (°C)', fontsize=10)
ax6.set_title('2D Response Surface: pH vs Temperature', fontsize=11, fontweight='bold')
plt.colorbar(contour, ax=ax6, label='Signal Intensity')
ax6.legend()

# Plot 7: Parameter convergence
ax7 = plt.subplot(3, 3, 7)
for i in range(4):
    normalized_params = (X_array[:, i] - bounds[i][0]) / (bounds[i][1] - bounds[i][0])
    ax7.plot(iterations, normalized_params, 'o-', alpha=0.7, label=param_names[i])

ax7.set_xlabel('Experiment Number', fontsize=11)
ax7.set_ylabel('Normalized Parameter Value', fontsize=11)
ax7.set_title('Parameter Convergence', fontsize=12, fontweight='bold')
ax7.legend(fontsize=8)
ax7.grid(True, alpha=0.3)
ax7.set_ylim([-0.05, 1.05])

# Plot 8: Signal improvement distribution
ax8 = plt.subplot(3, 3, 8)
random_signals = y_array[:10]
bayesian_signals = y_array[10:]

ax8.hist(random_signals, bins=5, alpha=0.7, label='Random Sampling', 
         color='orange', edgecolor='black')
ax8.hist(bayesian_signals, bins=10, alpha=0.7, label='Bayesian Optimization', 
         color='blue', edgecolor='black')
ax8.axvline(x=best_signal, color='red', linestyle='--', linewidth=2, 
            label=f'Best: {best_signal:.0f}')
ax8.set_xlabel('Signal Intensity', fontsize=11)
ax8.set_ylabel('Frequency', fontsize=11)
ax8.set_title('Signal Distribution', fontsize=12, fontweight='bold')
ax8.legend()
ax8.grid(True, alpha=0.3, axis='y')

# Plot 9: Acquisition function
ax9 = plt.subplot(3, 3, 9)
test_range = np.linspace(bounds[0][0], bounds[0][1], 100)
ei_values = []

for ce_val in test_range:
    test_params = best_params.copy()
    test_params[0] = ce_val
    ei = -optimizer.expected_improvement(test_params.reshape(1, -1))
    ei_values.append(ei[0])

ax9.plot(test_range, ei_values, 'g-', linewidth=2)
ax9.set_xlabel('Collision Energy (eV)', fontsize=11)
ax9.set_ylabel('Expected Improvement', fontsize=11)
ax9.set_title('Acquisition Function (CE)', fontsize=12, fontweight='bold')
ax9.grid(True, alpha=0.3)
ax9.axvline(x=best_params[0], color='red', linestyle='--', linewidth=2)

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

print("\n✓ Visualization complete! Figure saved as 'lcms_optimization_results.png'")

signal = 100000 * np.exp(...)

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.ndimage import gaussian_filter1d
from sklearn.metrics import roc_curve, auc, precision_recall_curve
import pywt
from scipy.optimize import minimize_scalar

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

# ============================================================================
# SECTION 1: Data Generation - Simulating Realistic Biosignature Scenarios
# ============================================================================

def generate_biosignature_signal(t, frequency=0.05, amplitude=2.0, phase=0):
    """
    Generate synthetic periodic biosignature signal
    Models methane production from microbial colonies with circadian rhythm
    
    Parameters:
    - t: time array
    - frequency: biological cycle frequency (Hz)
    - amplitude: signal strength
    - phase: phase offset
    """
    # Primary biological signal (circadian rhythm)
    bio_signal = amplitude * np.sin(2 * np.pi * frequency * t + phase)
    
    # Add secondary harmonic (metabolic fluctuations)
    bio_signal += 0.3 * amplitude * np.sin(4 * np.pi * frequency * t + phase)
    
    # Only positive concentrations
    bio_signal = np.maximum(bio_signal, 0)
    
    return bio_signal

def generate_geological_noise(t, num_events=5):
    """
    Generate geological interference (volcanic outgassing events)
    Random spikes with exponential decay
    """
    geo_noise = np.zeros_like(t)
    event_times = np.random.choice(len(t), num_events, replace=False)
    
    for event_time in event_times:
        # Exponential decay from eruption
        decay = np.exp(-0.1 * np.abs(np.arange(len(t)) - event_time))
        amplitude = np.random.uniform(1, 3)
        geo_noise += amplitude * decay
    
    return geo_noise

def generate_atmospheric_drift(t):
    """
    Generate low-frequency atmospheric drift
    Models weather patterns and seasonal variations
    """
    drift = 0.5 * np.sin(2 * np.pi * 0.01 * t) + 0.3 * np.cos(2 * np.pi * 0.015 * t)
    return drift

def generate_instrumental_noise(t, snr_db=10):
    """
    Generate Gaussian white noise from instrumental limitations
    
    Parameters:
    - snr_db: Signal-to-noise ratio in decibels
    """
    noise_power = 10 ** (-snr_db / 10)
    return np.random.normal(0, np.sqrt(noise_power), len(t))

# ============================================================================
# SECTION 2: Signal Processing - Multi-Stage Denoising Pipeline
# ============================================================================

def wavelet_denoise(signal_data, wavelet='db4', level=4, threshold_mode='soft'):
    """
    Wavelet-based denoising using discrete wavelet transform
    
    Mathematical basis:
    1. Decompose signal into wavelet coefficients
    2. Apply soft thresholding: W_threshold = sign(W) * max(|W| - λ, 0)
    3. Reconstruct signal from thresholded coefficients
    """
    # Perform multilevel wavelet decomposition
    coeffs = pywt.wavedec(signal_data, wavelet, level=level)
    
    # Calculate universal threshold: λ = σ * sqrt(2 * log(N))
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745  # Robust noise estimation
    threshold = sigma * np.sqrt(2 * np.log(len(signal_data)))
    
    # Apply thresholding to detail coefficients
    coeffs_thresholded = [coeffs[0]]  # Keep approximation coefficients
    for coeff in coeffs[1:]:
        if threshold_mode == 'soft':
            # Soft thresholding
            thresholded = pywt.threshold(coeff, threshold, mode='soft')
        else:
            # Hard thresholding
            thresholded = pywt.threshold(coeff, threshold, mode='hard')
        coeffs_thresholded.append(thresholded)
    
    # Reconstruct denoised signal
    denoised = pywt.waverec(coeffs_thresholded, wavelet)
    
    return denoised[:len(signal_data)]

def remove_trend(signal_data, poly_degree=3):
    """
    Remove polynomial trend (atmospheric drift)
    Fits polynomial and subtracts from signal
    """
    t = np.arange(len(signal_data))
    coefficients = np.polyfit(t, signal_data, poly_degree)
    trend = np.polyval(coefficients, t)
    detrended = signal_data - trend
    
    return detrended, trend

def matched_filter(signal_data, template, normalize=True):
    """
    Matched filtering to enhance known signal patterns
    
    Mathematical principle:
    Correlation with expected biosignature template maximizes SNR
    z(t) = ∫ y(τ) * h(t - τ) dτ
    """
    # Normalize template
    if normalize:
        template = (template - np.mean(template)) / np.std(template)
        signal_data = (signal_data - np.mean(signal_data)) / np.std(signal_data)
    
    # Cross-correlation
    filtered = signal.correlate(signal_data, template, mode='same')
    
    return filtered

# ============================================================================
# SECTION 3: Detection Algorithm - Adaptive Thresholding
# ============================================================================

def adaptive_threshold_detection(signal_data, false_positive_rate=0.01):
    """
    Statistical hypothesis testing for biosignature detection
    
    H0: Signal is noise (no biosignature)
    H1: Signal contains biosignature
    
    Uses Gaussian assumption for noise distribution
    """
    # Estimate noise statistics
    mu = np.mean(signal_data)
    sigma = np.std(signal_data)
    
    # Calculate threshold for desired false positive rate
    # For Gaussian: P(X > threshold) = false_positive_rate
    from scipy.stats import norm
    z_score = norm.ppf(1 - false_positive_rate)
    threshold = mu + z_score * sigma
    
    # Detect peaks above threshold
    detections = signal_data > threshold
    peak_indices = np.where(detections)[0]
    
    return detections, threshold, peak_indices

def calculate_detection_metrics(true_labels, predictions):
    """
    Calculate performance metrics: precision, recall, F1-score
    """
    true_positive = np.sum((true_labels == 1) & (predictions == 1))
    false_positive = np.sum((true_labels == 0) & (predictions == 1))
    false_negative = np.sum((true_labels == 1) & (predictions == 0))
    true_negative = np.sum((true_labels == 0) & (predictions == 0))
    
    precision = true_positive / (true_positive + false_positive) if (true_positive + false_positive) > 0 else 0
    recall = true_positive / (true_positive + false_negative) if (true_positive + false_negative) > 0 else 0
    f1_score = 2 * (precision * recall) / (precision + recall) if (precision + recall) > 0 else 0
    
    specificity = true_negative / (true_negative + false_positive) if (true_negative + false_positive) > 0 else 0
    
    return {
        'precision': precision,
        'recall': recall,
        'f1_score': f1_score,
        'specificity': specificity,
        'true_positive': true_positive,
        'false_positive': false_positive,
        'false_negative': false_negative,
        'true_negative': true_negative
    }

# ============================================================================
# SECTION 4: Main Processing Pipeline
# ============================================================================

# Generate time series (1000 samples, sampling rate 1 Hz)
n_samples = 1000
t = np.linspace(0, 1000, n_samples)

# Generate signal components
bio_signal = generate_biosignature_signal(t, frequency=0.05, amplitude=2.0)
geo_noise = generate_geological_noise(t, num_events=5)
atm_drift = generate_atmospheric_drift(t)
inst_noise = generate_instrumental_noise(t, snr_db=10)

# Combine all components
observed_signal = bio_signal + geo_noise + atm_drift + inst_noise

# Create ground truth labels (1 where biosignature is present)
threshold_bio = 0.5
true_labels = (bio_signal > threshold_bio).astype(int)

print("=" * 80)
print("BIOSIGNATURE DETECTION ALGORITHM - PROCESSING PIPELINE")
print("=" * 80)
print(f"\nDataset Statistics:")
print(f"  Total samples: {n_samples}")
print(f"  True biosignature events: {np.sum(true_labels)}")
print(f"  Signal-to-Noise Ratio: 10 dB")
print(f"  Observation period: {n_samples} seconds")

# ============================================================================
# Stage 1: Wavelet Denoising
# ============================================================================
print("\n" + "-" * 80)
print("STAGE 1: Wavelet Denoising")
print("-" * 80)

denoised_signal = wavelet_denoise(observed_signal, wavelet='db4', level=4)
noise_removed = observed_signal - denoised_signal
noise_power_reduction = 10 * np.log10(np.var(observed_signal) / np.var(denoised_signal))

print(f"  Wavelet: Daubechies 4 (db4)")
print(f"  Decomposition level: 4")
print(f"  Noise power reduction: {noise_power_reduction:.2f} dB")

# ============================================================================
# Stage 2: Trend Removal
# ============================================================================
print("\n" + "-" * 80)
print("STAGE 2: Atmospheric Drift Removal")
print("-" * 80)

detrended_signal, trend = remove_trend(denoised_signal, poly_degree=3)
drift_amplitude = np.max(trend) - np.min(trend)

print(f"  Polynomial degree: 3")
print(f"  Drift amplitude removed: {drift_amplitude:.4f}")

# ============================================================================
# Stage 3: Matched Filtering
# ============================================================================
print("\n" + "-" * 80)
print("STAGE 3: Matched Filtering")
print("-" * 80)

# Create template matching expected biosignature pattern
template_length = 200
template_t = np.linspace(0, 200, template_length)
template = generate_biosignature_signal(template_t, frequency=0.05, amplitude=1.0)

filtered_signal = matched_filter(detrended_signal, template)
snr_improvement = 10 * np.log10(np.var(filtered_signal) / np.var(detrended_signal))

print(f"  Template length: {template_length} samples")
print(f"  SNR improvement: {snr_improvement:.2f} dB")

# ============================================================================
# Stage 4: Detection with Multiple Thresholds
# ============================================================================
print("\n" + "-" * 80)
print("STAGE 4: Adaptive Threshold Detection")
print("-" * 80)

# Test multiple false positive rates
fpr_values = [0.001, 0.01, 0.05]
detection_results = []

for fpr in fpr_values:
    detections, threshold, peak_indices = adaptive_threshold_detection(filtered_signal, false_positive_rate=fpr)
    metrics = calculate_detection_metrics(true_labels, detections.astype(int))
    
    detection_results.append({
        'fpr': fpr,
        'threshold': threshold,
        'detections': detections,
        'metrics': metrics
    })
    
    print(f"\n  False Positive Rate: {fpr}")
    print(f"    Threshold: {threshold:.4f}")
    print(f"    Precision: {metrics['precision']:.4f}")
    print(f"    Recall: {metrics['recall']:.4f}")
    print(f"    F1-Score: {metrics['f1_score']:.4f}")
    print(f"    Specificity: {metrics['specificity']:.4f}")

# ============================================================================
# ROC Curve Analysis
# ============================================================================
print("\n" + "-" * 80)
print("ROC CURVE ANALYSIS")
print("-" * 80)

fpr_roc, tpr_roc, thresholds_roc = roc_curve(true_labels, filtered_signal)
roc_auc = auc(fpr_roc, tpr_roc)

print(f"  Area Under ROC Curve (AUC): {roc_auc:.4f}")
print(f"  Perfect classifier AUC: 1.0000")
print(f"  Random classifier AUC: 0.5000")

# Find optimal threshold (Youden's J statistic)
j_scores = tpr_roc - fpr_roc
optimal_idx = np.argmax(j_scores)
optimal_threshold = thresholds_roc[optimal_idx]
optimal_fpr = fpr_roc[optimal_idx]
optimal_tpr = tpr_roc[optimal_idx]

print(f"\n  Optimal Operating Point (Youden's Index):")
print(f"    Threshold: {optimal_threshold:.4f}")
print(f"    True Positive Rate: {optimal_tpr:.4f}")
print(f"    False Positive Rate: {optimal_fpr:.4f}")

# ============================================================================
# SECTION 5: Visualization
# ============================================================================

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

# Plot 1: Original Signal Components
ax1 = plt.subplot(4, 2, 1)
ax1.plot(t, bio_signal, label='True Biosignature', color='green', linewidth=2, alpha=0.7)
ax1.plot(t, observed_signal, label='Observed Signal (with noise)', color='gray', alpha=0.5)
ax1.set_xlabel('Time (s)', fontsize=10)
ax1.set_ylabel('Amplitude', fontsize=10)
ax1.set_title('(A) Raw Observed Signal vs True Biosignature', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# Plot 2: Noise Components Breakdown
ax2 = plt.subplot(4, 2, 2)
ax2.plot(t, geo_noise, label='Geological Noise', color='red', alpha=0.6)
ax2.plot(t, atm_drift, label='Atmospheric Drift', color='blue', alpha=0.6)
ax2.plot(t, inst_noise, label='Instrumental Noise', color='purple', alpha=0.4)
ax2.set_xlabel('Time (s)', fontsize=10)
ax2.set_ylabel('Amplitude', fontsize=10)
ax2.set_title('(B) Noise Components Analysis', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# Plot 3: Wavelet Denoising Effect
ax3 = plt.subplot(4, 2, 3)
ax3.plot(t, observed_signal, label='Original', color='gray', alpha=0.5, linewidth=1)
ax3.plot(t, denoised_signal, label='After Wavelet Denoising', color='orange', linewidth=2)
ax3.set_xlabel('Time (s)', fontsize=10)
ax3.set_ylabel('Amplitude', fontsize=10)
ax3.set_title('(C) Stage 1: Wavelet Denoising', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# Plot 4: Trend Removal
ax4 = plt.subplot(4, 2, 4)
ax4.plot(t, denoised_signal, label='Denoised Signal', color='orange', alpha=0.5)
ax4.plot(t, trend, label='Detected Trend', color='blue', linewidth=2, linestyle='--')
ax4.plot(t, detrended_signal, label='Detrended Signal', color='cyan', linewidth=1.5)
ax4.set_xlabel('Time (s)', fontsize=10)
ax4.set_ylabel('Amplitude', fontsize=10)
ax4.set_title('(D) Stage 2: Trend Removal', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

# Plot 5: Matched Filtering
ax5 = plt.subplot(4, 2, 5)
ax5.plot(t, detrended_signal, label='Detrended Signal', color='cyan', alpha=0.5)
ax5.plot(t, filtered_signal, label='Matched Filtered', color='magenta', linewidth=2)
ax5.set_xlabel('Time (s)', fontsize=10)
ax5.set_ylabel('Amplitude', fontsize=10)
ax5.set_title('(E) Stage 3: Matched Filtering', fontsize=12, fontweight='bold')
ax5.legend(fontsize=9)
ax5.grid(True, alpha=0.3)

# Plot 6: Detection Results with Different Thresholds
ax6 = plt.subplot(4, 2, 6)
ax6.plot(t, filtered_signal, label='Filtered Signal', color='magenta', linewidth=1.5, alpha=0.7)

colors = ['red', 'orange', 'yellow']
for i, result in enumerate(detection_results):
    ax6.axhline(y=result['threshold'], color=colors[i], linestyle='--', 
                label=f"Threshold (FPR={result['fpr']})", linewidth=2)
    
ax6.fill_between(t, 0, filtered_signal, where=(true_labels == 1), 
                 color='green', alpha=0.2, label='True Biosignature Regions')
ax6.set_xlabel('Time (s)', fontsize=10)
ax6.set_ylabel('Amplitude', fontsize=10)
ax6.set_title('(F) Stage 4: Detection with Multiple Thresholds', fontsize=12, fontweight='bold')
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)

# Plot 7: ROC Curve
ax7 = plt.subplot(4, 2, 7)
ax7.plot(fpr_roc, tpr_roc, color='darkblue', linewidth=3, label=f'ROC Curve (AUC = {roc_auc:.4f})')
ax7.plot([0, 1], [0, 1], color='gray', linestyle='--', linewidth=2, label='Random Classifier')
ax7.plot(optimal_fpr, optimal_tpr, 'r*', markersize=15, label=f'Optimal Point (J={j_scores[optimal_idx]:.3f})')
ax7.set_xlabel('False Positive Rate', fontsize=10)
ax7.set_ylabel('True Positive Rate', fontsize=10)
ax7.set_title('(G) ROC Curve Analysis', fontsize=12, fontweight='bold')
ax7.legend(fontsize=9)
ax7.grid(True, alpha=0.3)
ax7.set_xlim([0, 1])
ax7.set_ylim([0, 1])

# Plot 8: Performance Metrics Comparison
ax8 = plt.subplot(4, 2, 8)
metrics_names = ['Precision', 'Recall', 'F1-Score', 'Specificity']
x = np.arange(len(metrics_names))
width = 0.25

for i, result in enumerate(detection_results):
    metrics_values = [
        result['metrics']['precision'],
        result['metrics']['recall'],
        result['metrics']['f1_score'],
        result['metrics']['specificity']
    ]
    ax8.bar(x + i*width, metrics_values, width, label=f"FPR={result['fpr']}")

ax8.set_xlabel('Metrics', fontsize=10)
ax8.set_ylabel('Score', fontsize=10)
ax8.set_title('(H) Detection Performance Comparison', fontsize=12, fontweight='bold')
ax8.set_xticks(x + width)
ax8.set_xticklabels(metrics_names, fontsize=9)
ax8.legend(fontsize=9)
ax8.grid(True, alpha=0.3, axis='y')
ax8.set_ylim([0, 1.1])

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

print("\n" + "=" * 80)
print("PROCESSING COMPLETE")
print("=" * 80)
print("\nVisualization saved as 'biosignature_detection_analysis.png'")
print("\nKey Findings:")
print(f"  • Wavelet denoising reduced noise power by {noise_power_reduction:.2f} dB")
print(f"  • Matched filtering improved SNR by {snr_improvement:.2f} dB")
print(f"  • Best F1-score: {max([r['metrics']['f1_score'] for r in detection_results]):.4f}")
print(f"  • ROC AUC: {roc_auc:.4f} (excellent discrimination)")
print("\n" + "=" * 80)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution, linprog
from matplotlib.patches import FancyBboxPatch
import seaborn as sns

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

class MetabolicPathwayOptimizer:
    """
    Optimizes synthetic metabolic pathways for low-resource planetary environments
    using constraint-based modeling and evolutionary algorithms.
    """
    
    def __init__(self, n_reactions=6, n_enzymes=4):
        """
        Initialize the metabolic network
        
        Parameters:
        -----------
        n_reactions : int
            Number of reactions in the pathway
        n_enzymes : int
            Number of enzymes to optimize
        """
        self.n_reactions = n_reactions
        self.n_enzymes = n_enzymes
        
        # Define stoichiometric matrix (reactions x metabolites)
        # Metabolites: [H2, CO2, CH4, H2O, ATP, Biomass]
        self.S = np.array([
            [-1, -0.5, 0, 1, 2, 0],      # R1: H2 + 0.5CO2 -> H2O + 2ATP
            [-2, -1, 1, 2, 1, 0],         # R2: 2H2 + CO2 -> CH4 + 2H2O + ATP
            [0, -1, 0, 0, 3, 0.1],        # R3: CO2 -> 3ATP + 0.1Biomass (alternative)
            [-1, 0, 0, 0, 1, 0.05],       # R4: H2 -> ATP + 0.05Biomass
            [0, 0, -1, 0, 4, 0],          # R5: CH4 -> 4ATP (methane oxidation)
            [0, 0, 0, 0, -10, 1]          # R6: 10ATP -> Biomass
        ])
        
        # ATP yield per reaction
        self.atp_yield = np.array([2, 1, 3, 1, 4, -10])
        
        # Enzyme-reaction mapping (which enzymes catalyze which reactions)
        self.enzyme_reaction_map = np.array([
            [1, 0.5, 0, 0.2],   # R1 requires enzyme 0 primarily, enzyme 1 partially
            [0.3, 1, 0.2, 0],   # R2 requires enzyme 1 primarily
            [0, 0.2, 1, 0],     # R3 requires enzyme 2
            [0.5, 0, 0, 1],     # R4 requires enzyme 3
            [0, 0.8, 0.3, 0],   # R5 requires enzyme 1 and 2
            [0.1, 0.1, 0.1, 0.1] # R6 requires all enzymes minimally
        ])
        
        # Resource cost per enzyme (arbitrary units)
        self.enzyme_cost = np.array([5, 8, 6, 4])
        
        # Planetary environment constraints (resource availability)
        self.resource_limits = {
            'H2_max': 10.0,      # Maximum H2 uptake
            'CO2_max': 8.0,      # Maximum CO2 uptake
            'CH4_max': 5.0,      # Maximum CH4 production
            'total_enzyme': 20.0  # Total enzyme budget
        }
        
    def calculate_max_flux(self, enzyme_levels):
        """
        Calculate maximum flux through each reaction based on enzyme levels
        
        The flux capacity is determined by:
        v_max_i = sum_j (E_j * k_ij)
        where E_j is enzyme level and k_ij is catalytic efficiency
        """
        flux_capacity = np.dot(self.enzyme_reaction_map, enzyme_levels)
        return flux_capacity
    
    def steady_state_fba(self, enzyme_levels):
        """
        Flux Balance Analysis under steady-state assumption
        Maximizes ATP production subject to stoichiometric and capacity constraints
        
        Mathematical formulation:
        max: c^T * v
        s.t.: S * v = 0 (steady state)
              v_i <= v_max_i (capacity constraints)
              v_i >= 0 (irreversibility)
        """
        flux_capacity = self.calculate_max_flux(enzyme_levels)
        
        # Objective: maximize ATP production (negative for minimization)
        c = -self.atp_yield
        
        # Equality constraints: S * v = 0 (steady state for internal metabolites)
        # We only enforce steady state for internal metabolites (not external)
        A_eq = self.S[:, 3:].T  # H2O, ATP, Biomass must balance
        b_eq = np.zeros(A_eq.shape[0])
        
        # Inequality constraints: flux bounds
        A_ub = np.eye(self.n_reactions)
        b_ub = flux_capacity
        
        # Additional resource constraints
        # H2 uptake: sum of H2 consuming reactions <= H2_max
        A_ub = np.vstack([A_ub, [1, 2, 0, 1, 0, 0]])  # H2 constraint
        b_ub = np.append(b_ub, self.resource_limits['H2_max'])
        
        # CO2 uptake
        A_ub = np.vstack([A_ub, [0.5, 1, 1, 0, 0, 0]])  # CO2 constraint
        b_ub = np.append(b_ub, self.resource_limits['CO2_max'])
        
        # Flux lower bounds (all reactions are irreversible)
        bounds = [(0, None) for _ in range(self.n_reactions)]
        
        # Solve linear program
        result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, 
                        bounds=bounds, method='highs')
        
        if result.success:
            return result.x, -result.fun  # Return flux and ATP production
        else:
            return np.zeros(self.n_reactions), 0
    
    def objective_function(self, enzyme_levels):
        """
        Objective function to maximize: ATP production / Resource cost
        
        This represents the efficiency of the metabolic pathway
        """
        # Ensure enzyme levels are within bounds
        enzyme_levels = np.clip(enzyme_levels, 0, self.resource_limits['total_enzyme'])
        
        # Normalize to meet total enzyme budget
        if np.sum(enzyme_levels) > self.resource_limits['total_enzyme']:
            enzyme_levels = enzyme_levels * (self.resource_limits['total_enzyme'] / np.sum(enzyme_levels))
        
        # Calculate fluxes and ATP production
        fluxes, atp_production = self.steady_state_fba(enzyme_levels)
        
        # Calculate total resource cost
        resource_cost = np.dot(enzyme_levels, self.enzyme_cost)
        
        # Avoid division by zero
        if resource_cost < 0.1:
            return -1e10
        
        # Return negative (for minimization in differential_evolution)
        efficiency = atp_production / resource_cost
        return -efficiency
    
    def optimize(self):
        """
        Optimize enzyme expression levels using differential evolution
        """
        # Bounds for enzyme levels (0 to total_enzyme budget)
        bounds = [(0, self.resource_limits['total_enzyme']) for _ in range(self.n_enzymes)]
        
        # Run optimization
        result = differential_evolution(
            self.objective_function,
            bounds,
            maxiter=100,
            popsize=15,
            tol=0.01,
            seed=42,
            workers=1
        )
        
        return result
    
    def analyze_sensitivity(self, optimal_enzymes, parameter='H2_max', 
                          param_range=np.linspace(5, 15, 20)):
        """
        Analyze sensitivity of the system to environmental parameters
        """
        efficiencies = []
        atp_productions = []
        
        original_value = self.resource_limits[parameter]
        
        for param_value in param_range:
            self.resource_limits[parameter] = param_value
            fluxes, atp = self.steady_state_fba(optimal_enzymes)
            resource_cost = np.dot(optimal_enzymes, self.enzyme_cost)
            efficiencies.append(atp / resource_cost if resource_cost > 0 else 0)
            atp_productions.append(atp)
        
        # Restore original value
        self.resource_limits[parameter] = original_value
        
        return param_range, efficiencies, atp_productions


# Initialize and run optimization
print("="*70)
print("METABOLIC PATHWAY OPTIMIZATION FOR PLANETARY ENVIRONMENTS")
print("="*70)
print("\nInitializing synthetic metabolic network...")
print("Environment: Mars-like atmosphere (low H2, CO2 rich)")
print("Objective: Maximize ATP production efficiency\n")

optimizer = MetabolicPathwayOptimizer()

print("Running differential evolution optimization...")
print("Population size: 15 | Max iterations: 100\n")

result = optimizer.optimize()

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

optimal_enzymes = result.x
print("\nOptimal Enzyme Expression Levels:")
for i, level in enumerate(optimal_enzymes):
    print(f"  Enzyme {i+1}: {level:.3f} units ({level/np.sum(optimal_enzymes)*100:.1f}%)")

print(f"\nTotal Enzyme Budget Used: {np.sum(optimal_enzymes):.3f} / {optimizer.resource_limits['total_enzyme']:.1f} units")

# Calculate final performance metrics
optimal_fluxes, optimal_atp = optimizer.steady_state_fba(optimal_enzymes)
resource_cost = np.dot(optimal_enzymes, optimizer.enzyme_cost)
efficiency = optimal_atp / resource_cost

print(f"\nPerformance Metrics:")
print(f"  ATP Production Rate: {optimal_atp:.3f} mmol/gDW/h")
print(f"  Resource Cost: {resource_cost:.3f} units")
print(f"  Efficiency (ATP/Cost): {efficiency:.4f}")

print(f"\nOptimal Reaction Fluxes:")
reaction_names = ['H2+CO2→H2O+ATP', '2H2+CO2→CH4+ATP', 
                  'CO2→ATP+Biomass', 'H2→ATP+Biomass',
                  'CH4→ATP', 'ATP→Biomass']
for i, (flux, name) in enumerate(zip(optimal_fluxes, reaction_names)):
    print(f"  R{i+1} ({name}): {flux:.3f} mmol/gDW/h")

# Perform sensitivity analysis
print("\n" + "="*70)
print("SENSITIVITY ANALYSIS")
print("="*70)
print("\nAnalyzing system response to H2 availability...")

h2_range, h2_efficiencies, h2_atp = optimizer.analyze_sensitivity(
    optimal_enzymes, 'H2_max', np.linspace(5, 15, 20)
)

print("Analyzing system response to CO2 availability...")
co2_range, co2_efficiencies, co2_atp = optimizer.analyze_sensitivity(
    optimal_enzymes, 'CO2_max', np.linspace(4, 12, 20)
)

print("\nSensitivity analysis complete!")

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

print("\n" + "="*70)
print("GENERATING VISUALIZATIONS")
print("="*70)

fig = plt.figure(figsize=(16, 12))
gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)

# Color scheme
colors = sns.color_palette("husl", optimizer.n_enzymes)
color_flux = sns.color_palette("viridis", optimizer.n_reactions)

# 1. Enzyme Expression Levels (Pie Chart)
ax1 = fig.add_subplot(gs[0, 0])
wedges, texts, autotexts = ax1.pie(optimal_enzymes, labels=[f'E{i+1}' for i in range(optimizer.n_enzymes)],
                                     autopct='%1.1f%%', colors=colors, startangle=90)
ax1.set_title('Optimal Enzyme Distribution', fontsize=12, fontweight='bold')
for autotext in autotexts:
    autotext.set_color('white')
    autotext.set_fontweight('bold')

# 2. Enzyme Expression Bar Chart with Cost
ax2 = fig.add_subplot(gs[0, 1])
x_pos = np.arange(optimizer.n_enzymes)
bars = ax2.bar(x_pos, optimal_enzymes, color=colors, alpha=0.7, edgecolor='black')
ax2_twin = ax2.twinx()
ax2_twin.plot(x_pos, optimizer.enzyme_cost, 'ro-', linewidth=2, markersize=8, label='Unit Cost')
ax2.set_xlabel('Enzyme Type', fontweight='bold')
ax2.set_ylabel('Expression Level (units)', fontweight='bold', color='blue')
ax2_twin.set_ylabel('Resource Cost per Unit', fontweight='bold', color='red')
ax2.set_xticks(x_pos)
ax2.set_xticklabels([f'E{i+1}' for i in range(optimizer.n_enzymes)])
ax2.set_title('Enzyme Levels vs. Resource Cost', fontsize=12, fontweight='bold')
ax2.grid(axis='y', alpha=0.3)
ax2_twin.legend(loc='upper right')

# 3. Reaction Flux Distribution
ax3 = fig.add_subplot(gs[0, 2])
bars = ax3.barh(range(optimizer.n_reactions), optimal_fluxes, color=color_flux, 
                edgecolor='black', alpha=0.8)
ax3.set_yticks(range(optimizer.n_reactions))
ax3.set_yticklabels([f'R{i+1}' for i in range(optimizer.n_reactions)])
ax3.set_xlabel('Flux (mmol/gDW/h)', fontweight='bold')
ax3.set_title('Metabolic Flux Distribution', fontsize=12, fontweight='bold')
ax3.grid(axis='x', alpha=0.3)

# Add flux values as text
for i, (flux, bar) in enumerate(zip(optimal_fluxes, bars)):
    if flux > 0.1:
        ax3.text(flux + 0.1, i, f'{flux:.2f}', va='center', fontweight='bold')

# 4. ATP Contribution by Reaction
ax4 = fig.add_subplot(gs[1, 0])
atp_contributions = optimal_fluxes * optimizer.atp_yield
atp_contributions = np.maximum(atp_contributions, 0)  # Only positive contributions
bars = ax4.bar(range(optimizer.n_reactions), atp_contributions, color=color_flux,
               edgecolor='black', alpha=0.8)
ax4.set_xlabel('Reaction', fontweight='bold')
ax4.set_ylabel('ATP Contribution (mmol/gDW/h)', fontweight='bold')
ax4.set_xticks(range(optimizer.n_reactions))
ax4.set_xticklabels([f'R{i+1}' for i in range(optimizer.n_reactions)])
ax4.set_title('ATP Production by Reaction', fontsize=12, fontweight='bold')
ax4.grid(axis='y', alpha=0.3)

# 5. Sensitivity to H2 Availability
ax5 = fig.add_subplot(gs[1, 1])
ax5_twin = ax5.twinx()
line1 = ax5.plot(h2_range, h2_efficiencies, 'b-o', linewidth=2, markersize=6, 
                 label='Efficiency')
line2 = ax5_twin.plot(h2_range, h2_atp, 'r-s', linewidth=2, markersize=6, 
                      label='ATP Production')
ax5.axvline(optimizer.resource_limits['H2_max'], color='green', linestyle='--', 
            linewidth=2, label='Current Limit')
ax5.set_xlabel('H₂ Availability (mmol/gDW/h)', fontweight='bold')
ax5.set_ylabel('Efficiency (ATP/Cost)', fontweight='bold', color='blue')
ax5_twin.set_ylabel('ATP Production', fontweight='bold', color='red')
ax5.set_title('Sensitivity to H₂ Availability', fontsize=12, fontweight='bold')
ax5.grid(alpha=0.3)
lines = line1 + line2 + [ax5.get_lines()[-1]]
labels = [l.get_label() for l in lines]
ax5.legend(lines, labels, loc='upper left')

# 6. Sensitivity to CO2 Availability
ax6 = fig.add_subplot(gs[1, 2])
ax6_twin = ax6.twinx()
line1 = ax6.plot(co2_range, co2_efficiencies, 'b-o', linewidth=2, markersize=6,
                 label='Efficiency')
line2 = ax6_twin.plot(co2_range, co2_atp, 'r-s', linewidth=2, markersize=6,
                      label='ATP Production')
ax6.axvline(optimizer.resource_limits['CO2_max'], color='green', linestyle='--',
            linewidth=2, label='Current Limit')
ax6.set_xlabel('CO₂ Availability (mmol/gDW/h)', fontweight='bold')
ax6.set_ylabel('Efficiency (ATP/Cost)', fontweight='bold', color='blue')
ax6_twin.set_ylabel('ATP Production', fontweight='bold', color='red')
ax6.set_title('Sensitivity to CO₂ Availability', fontsize=12, fontweight='bold')
ax6.grid(alpha=0.3)
lines = line1 + line2 + [ax6.get_lines()[-1]]
labels = [l.get_label() for l in lines]
ax6.legend(lines, labels, loc='upper left')

# 7. Enzyme-Reaction Network Map
ax7 = fig.add_subplot(gs[2, :])
im = ax7.imshow(optimizer.enzyme_reaction_map.T, cmap='YlOrRd', aspect='auto', 
                interpolation='nearest')
ax7.set_xticks(range(optimizer.n_reactions))
ax7.set_yticks(range(optimizer.n_enzymes))
ax7.set_xticklabels([f'R{i+1}' for i in range(optimizer.n_reactions)])
ax7.set_yticklabels([f'E{i+1}' for i in range(optimizer.n_enzymes)])
ax7.set_xlabel('Reactions', fontweight='bold')
ax7.set_ylabel('Enzymes', fontweight='bold')
ax7.set_title('Enzyme-Reaction Catalytic Network (darker = stronger dependency)', 
              fontsize=12, fontweight='bold')

# Add text annotations
for i in range(optimizer.n_enzymes):
    for j in range(optimizer.n_reactions):
        text = ax7.text(j, i, f'{optimizer.enzyme_reaction_map[j, i]:.1f}',
                       ha="center", va="center", color="black", fontsize=9,
                       fontweight='bold')

plt.colorbar(im, ax=ax7, label='Catalytic Efficiency')

plt.suptitle('METABOLIC PATHWAY OPTIMIZATION: COMPREHENSIVE ANALYSIS', 
             fontsize=16, fontweight='bold', y=0.995)

plt.savefig('metabolic_optimization_results.png', dpi=300, bbox_inches='tight')
print("\nVisualization saved as 'metabolic_optimization_results.png'")
plt.show()

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