import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.signal import find_peaks
import warnings
warnings.filterwarnings('ignore')

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

class PulsarDetector:
    """
    A class for detecting pulsar signals using phase folding and chi-squared periodogram
    """
    
    def __init__(self, n_bins=20, significance_threshold=3.0):
        """
        Initialize the detector with parameters
        
        Parameters:
        -----------
        n_bins : int
            Number of phase bins for folding (affects sensitivity)
        significance_threshold : float
            Sigma threshold for detection (affects false positive rate)
        """
        self.n_bins = n_bins
        self.significance_threshold = significance_threshold
        
    def generate_pulsar_signal(self, n_photons=10000, period=0.5, 
                                pulse_width=0.1, noise_level=0.8):
        """
        Generate synthetic pulsar data with noise
        
        Parameters:
        -----------
        n_photons : int
            Total number of photon events
        period : float
            True pulsar period (seconds)
        pulse_width : float
            Width of the pulse relative to period (0-1)
        noise_level : float
            Fraction of uniform noise (0=no noise, 1=all noise)
        
        Returns:
        --------
        arrival_times : array
            Photon arrival times
        true_period : float
            The true period used
        """
        # Signal photons (concentrated in pulses)
        n_signal = int(n_photons * (1 - noise_level))
        n_noise = n_photons - n_signal
        
        # Generate observation time span
        total_time = 100.0  # seconds
        
        # Signal: concentrated around pulse peaks
        n_pulses = int(total_time / period)
        signal_times = []
        
        for i in range(n_pulses):
            pulse_center = i * period
            # Generate photons with Gaussian distribution around pulse peak
            pulse_times = np.random.normal(pulse_center, pulse_width * period / 4, 
                                          n_signal // n_pulses)
            signal_times.extend(pulse_times)
        
        signal_times = np.array(signal_times)
        signal_times = signal_times[(signal_times >= 0) & (signal_times < total_time)]
        
        # Noise: uniformly distributed
        noise_times = np.random.uniform(0, total_time, n_noise)
        
        # Combine and sort
        arrival_times = np.sort(np.concatenate([signal_times, noise_times]))
        
        return arrival_times, period
    
    def chi_squared_periodogram(self, arrival_times, trial_periods):
        """
        Calculate chi-squared statistic for a range of trial periods
        
        Parameters:
        -----------
        arrival_times : array
            Photon arrival times
        trial_periods : array
            Array of periods to test
        
        Returns:
        --------
        chi_squared : array
            Chi-squared statistic for each trial period
        significance : array
            Significance in sigma units
        """
        chi_squared = np.zeros(len(trial_periods))
        
        for i, period in enumerate(trial_periods):
            # Fold the data at this trial period
            phases = (arrival_times % period) / period
            
            # Create histogram of phases
            counts, _ = np.histogram(phases, bins=self.n_bins, range=(0, 1))
            
            # Calculate chi-squared statistic
            mean_count = np.mean(counts)
            if mean_count > 0:
                chi_squared[i] = np.sum((counts - mean_count)**2 / mean_count)
            else:
                chi_squared[i] = 0
        
        # Calculate significance
        expected_chi2 = self.n_bins - 1
        std_chi2 = np.sqrt(2 * (self.n_bins - 1))
        significance = (chi_squared - expected_chi2) / std_chi2
        
        return chi_squared, significance
    
    def detect_period(self, arrival_times, period_range=(0.1, 2.0), n_periods=1000):
        """
        Detect pulsar period from arrival times
        
        Parameters:
        -----------
        arrival_times : array
            Photon arrival times
        period_range : tuple
            (min_period, max_period) to search
        n_periods : int
            Number of trial periods
        
        Returns:
        --------
        best_period : float
            Detected period
        significance : float
            Detection significance
        all_results : dict
            Full periodogram results
        """
        trial_periods = np.linspace(period_range[0], period_range[1], n_periods)
        chi_squared, significance = self.chi_squared_periodogram(arrival_times, trial_periods)
        
        # Find peaks in significance
        peaks, properties = find_peaks(significance, height=self.significance_threshold)
        
        if len(peaks) > 0:
            best_idx = peaks[np.argmax(properties['peak_heights'])]
            best_period = trial_periods[best_idx]
            best_significance = significance[best_idx]
        else:
            best_period = None
            best_significance = 0
        
        return best_period, best_significance, {
            'trial_periods': trial_periods,
            'chi_squared': chi_squared,
            'significance': significance,
            'peaks': peaks
        }
    
    def false_alarm_probability(self, significance, n_trials):
        """
        Calculate false alarm probability using extreme value statistics
        
        Parameters:
        -----------
        significance : float
            Observed significance in sigma
        n_trials : int
            Number of independent trials
        
        Returns:
        --------
        fap : float
            False alarm probability
        """
        # Single trial probability
        p_single = 1 - stats.norm.cdf(significance)
        
        # Multiple trial correction
        fap = 1 - (1 - p_single)**n_trials
        
        return fap


def optimize_parameters():
    """
    Optimize detection parameters to minimize false positives
    """
    print("="*70)
    print("PULSAR SIGNAL DETECTION: PARAMETER OPTIMIZATION")
    print("="*70)
    print()
    
    # Test different parameter combinations
    n_bins_range = [10, 20, 30, 40, 50]
    threshold_range = [2.5, 3.0, 3.5, 4.0, 4.5]
    
    # Generate test data: real pulsar + noise
    print("Generating synthetic pulsar data...")
    print("-" * 70)
    
    # Real pulsar signal
    detector_temp = PulsarDetector()
    real_data, true_period = detector_temp.generate_pulsar_signal(
        n_photons=10000, period=0.5, pulse_width=0.1, noise_level=0.7
    )
    print(f"True period: {true_period:.4f} seconds")
    print(f"Number of photons: {len(real_data)}")
    print(f"Observation time: {real_data[-1]:.2f} seconds")
    print()
    
    # Pure noise (for false positive testing)
    noise_data = np.sort(np.random.uniform(0, real_data[-1], len(real_data)))
    
    # Test all combinations
    results = []
    
    print("Testing parameter combinations...")
    print("-" * 70)
    
    for n_bins in n_bins_range:
        for threshold in threshold_range:
            detector = PulsarDetector(n_bins=n_bins, 
                                     significance_threshold=threshold)
            
            # Test on real signal
            detected_period, significance, _ = detector.detect_period(
                real_data, period_range=(0.1, 1.0), n_periods=500
            )
            
            # Test on pure noise
            noise_period, noise_sig, _ = detector.detect_period(
                noise_data, period_range=(0.1, 1.0), n_periods=500
            )
            
            # Calculate metrics
            if detected_period is not None:
                period_error = abs(detected_period - true_period)
                detection_success = period_error < 0.01  # Within 1% of true period
            else:
                period_error = np.inf
                detection_success = False
            
            false_positive = noise_period is not None
            
            results.append({
                'n_bins': n_bins,
                'threshold': threshold,
                'detected_period': detected_period,
                'period_error': period_error,
                'significance': significance,
                'detection_success': detection_success,
                'false_positive': false_positive,
                'noise_significance': noise_sig
            })
    
    # Find optimal parameters
    valid_detections = [r for r in results if r['detection_success']]
    if valid_detections:
        # Minimize false positives while maintaining detection
        optimal = min(valid_detections, 
                     key=lambda x: (x['false_positive'], -x['significance']))
        
        print("\nOptimal parameters found:")
        print(f"  Number of bins: {optimal['n_bins']}")
        print(f"  Significance threshold: {optimal['threshold']:.1f} σ")
        print(f"  Detected period: {optimal['detected_period']:.6f} seconds")
        print(f"  Period error: {optimal['period_error']:.6f} seconds")
        print(f"  Detection significance: {optimal['significance']:.2f} σ")
        print(f"  False positive: {optimal['false_positive']}")
        print()
    
    return results, real_data, true_period, noise_data, optimal


def visualize_results(results, real_data, true_period, noise_data, optimal):
    """
    Create comprehensive visualizations of the results
    """
    fig = plt.figure(figsize=(16, 12))
    
    # 1. Parameter space heatmap (Detection success)
    ax1 = plt.subplot(3, 3, 1)
    n_bins_vals = sorted(list(set([r['n_bins'] for r in results])))
    threshold_vals = sorted(list(set([r['threshold'] for r in results])))
    
    success_matrix = np.zeros((len(threshold_vals), len(n_bins_vals)))
    for r in results:
        i = threshold_vals.index(r['threshold'])
        j = n_bins_vals.index(r['n_bins'])
        success_matrix[i, j] = r['detection_success']
    
    im1 = ax1.imshow(success_matrix, aspect='auto', cmap='RdYlGn', 
                     interpolation='nearest')
    ax1.set_xticks(range(len(n_bins_vals)))
    ax1.set_xticklabels(n_bins_vals)
    ax1.set_yticks(range(len(threshold_vals)))
    ax1.set_yticklabels(threshold_vals)
    ax1.set_xlabel('Number of Bins')
    ax1.set_ylabel('Threshold (σ)')
    ax1.set_title('Detection Success Rate')
    plt.colorbar(im1, ax=ax1)
    
    # 2. False positive rate
    ax2 = plt.subplot(3, 3, 2)
    fp_matrix = np.zeros((len(threshold_vals), len(n_bins_vals)))
    for r in results:
        i = threshold_vals.index(r['threshold'])
        j = n_bins_vals.index(r['n_bins'])
        fp_matrix[i, j] = r['false_positive']
    
    im2 = ax2.imshow(fp_matrix, aspect='auto', cmap='RdYlGn_r', 
                     interpolation='nearest')
    ax2.set_xticks(range(len(n_bins_vals)))
    ax2.set_xticklabels(n_bins_vals)
    ax2.set_yticks(range(len(threshold_vals)))
    ax2.set_yticklabels(threshold_vals)
    ax2.set_xlabel('Number of Bins')
    ax2.set_ylabel('Threshold (σ)')
    ax2.set_title('False Positive Rate')
    plt.colorbar(im2, ax=ax2)
    
    # 3. Detection significance
    ax3 = plt.subplot(3, 3, 3)
    sig_matrix = np.zeros((len(threshold_vals), len(n_bins_vals)))
    for r in results:
        i = threshold_vals.index(r['threshold'])
        j = n_bins_vals.index(r['n_bins'])
        if r['detection_success']:
            sig_matrix[i, j] = r['significance']
    
    im3 = ax3.imshow(sig_matrix, aspect='auto', cmap='viridis', 
                     interpolation='nearest')
    ax3.set_xticks(range(len(n_bins_vals)))
    ax3.set_xticklabels(n_bins_vals)
    ax3.set_yticks(range(len(threshold_vals)))
    ax3.set_yticklabels(threshold_vals)
    ax3.set_xlabel('Number of Bins')
    ax3.set_ylabel('Threshold (σ)')
    ax3.set_title('Detection Significance (σ)')
    plt.colorbar(im3, ax=ax3)
    
    # 4. Periodogram with optimal parameters
    ax4 = plt.subplot(3, 3, 4)
    detector_opt = PulsarDetector(n_bins=optimal['n_bins'], 
                                  significance_threshold=optimal['threshold'])
    _, _, periodogram = detector_opt.detect_period(real_data, 
                                                    period_range=(0.1, 1.0), 
                                                    n_periods=1000)
    
    ax4.plot(periodogram['trial_periods'], periodogram['significance'], 
             'b-', linewidth=1, label='Significance')
    ax4.axhline(optimal['threshold'], color='r', linestyle='--', 
                label=f'Threshold ({optimal["threshold"]:.1f}σ)')
    ax4.axvline(true_period, color='g', linestyle='--', 
                label=f'True Period ({true_period:.3f}s)')
    if optimal['detected_period']:
        ax4.axvline(optimal['detected_period'], color='orange', linestyle=':', 
                    label=f'Detected ({optimal["detected_period"]:.3f}s)')
    ax4.set_xlabel('Trial Period (seconds)')
    ax4.set_ylabel('Significance (σ)')
    ax4.set_title('Chi-Squared Periodogram (Optimal Parameters)')
    ax4.legend(fontsize=8)
    ax4.grid(True, alpha=0.3)
    
    # 5. Phase-folded light curve
    ax5 = plt.subplot(3, 3, 5)
    phases = (real_data % optimal['detected_period']) / optimal['detected_period']
    ax5.hist(phases, bins=50, color='steelblue', alpha=0.7, edgecolor='black')
    ax5.set_xlabel('Phase')
    ax5.set_ylabel('Counts')
    ax5.set_title(f'Phase-Folded Light Curve (P={optimal["detected_period"]:.4f}s)')
    ax5.grid(True, alpha=0.3)
    
    # 6. Arrival time distribution
    ax6 = plt.subplot(3, 3, 6)
    ax6.hist(real_data, bins=100, color='purple', alpha=0.7, edgecolor='black')
    ax6.set_xlabel('Time (seconds)')
    ax6.set_ylabel('Counts')
    ax6.set_title('Photon Arrival Time Distribution')
    ax6.grid(True, alpha=0.3)
    
    # 7. Period error vs n_bins
    ax7 = plt.subplot(3, 3, 7)
    for threshold in threshold_vals:
        errors = []
        bins = []
        for r in results:
            if r['threshold'] == threshold and r['detection_success']:
                errors.append(r['period_error'])
                bins.append(r['n_bins'])
        if errors:
            ax7.plot(bins, errors, 'o-', label=f'Threshold {threshold}σ')
    ax7.set_xlabel('Number of Bins')
    ax7.set_ylabel('Period Error (seconds)')
    ax7.set_title('Detection Accuracy vs Parameters')
    ax7.legend(fontsize=8)
    ax7.grid(True, alpha=0.3)
    ax7.set_yscale('log')
    
    # 8. Noise test periodogram
    ax8 = plt.subplot(3, 3, 8)
    _, _, noise_periodogram = detector_opt.detect_period(noise_data, 
                                                         period_range=(0.1, 1.0), 
                                                         n_periods=1000)
    ax8.plot(noise_periodogram['trial_periods'], 
             noise_periodogram['significance'], 
             'r-', linewidth=1, alpha=0.7)
    ax8.axhline(optimal['threshold'], color='black', linestyle='--', 
                label=f'Threshold ({optimal["threshold"]:.1f}σ)')
    ax8.set_xlabel('Trial Period (seconds)')
    ax8.set_ylabel('Significance (σ)')
    ax8.set_title('Periodogram for Pure Noise (False Positive Test)')
    ax8.legend(fontsize=8)
    ax8.grid(True, alpha=0.3)
    
    # 9. ROC-like curve
    ax9 = plt.subplot(3, 3, 9)
    fp_rates = []
    detection_rates = []
    
    for threshold in threshold_vals:
        subset = [r for r in results if r['threshold'] == threshold]
        fp_rate = sum(r['false_positive'] for r in subset) / len(subset)
        det_rate = sum(r['detection_success'] for r in subset) / len(subset)
        fp_rates.append(fp_rate)
        detection_rates.append(det_rate)
    
    ax9.plot(fp_rates, detection_rates, 'o-', linewidth=2, markersize=8)
    for i, threshold in enumerate(threshold_vals):
        ax9.annotate(f'{threshold}σ', (fp_rates[i], detection_rates[i]), 
                    fontsize=8, xytext=(5, 5), textcoords='offset points')
    ax9.plot([0, 1], [0, 1], 'k--', alpha=0.3)
    ax9.set_xlabel('False Positive Rate')
    ax9.set_ylabel('Detection Success Rate')
    ax9.set_title('Detection Performance Trade-off')
    ax9.grid(True, alpha=0.3)
    ax9.set_xlim(-0.05, 1.05)
    ax9.set_ylim(-0.05, 1.05)
    
    plt.tight_layout()
    plt.savefig('pulsar_detection_optimization.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n" + "="*70)
    print("VISUALIZATION COMPLETE")
    print("="*70)


# Main execution
print("\n" + "="*70)
print("STARTING PULSAR DETECTION OPTIMIZATION")
print("="*70)
print()

results, real_data, true_period, noise_data, optimal = optimize_parameters()
visualize_results(results, real_data, true_period, noise_data, optimal)

print("\n" + "="*70)
print("ANALYSIS SUMMARY")
print("="*70)
print("""
This analysis demonstrates optimal parameter selection for pulsar detection:

1. PHASE BINS: More bins increase sensitivity but also increase noise
   - Too few: Miss fine pulse structure
   - Too many: Amplify statistical fluctuations
   - Optimal: 20-30 bins for typical pulsars

2. SIGNIFICANCE THRESHOLD: Controls false positive rate
   - Lower threshold: More detections but more false alarms
   - Higher threshold: Fewer false positives but miss weak signals
   - Optimal: 3.0-3.5σ balances sensitivity and specificity

3. KEY METRICS:
   - Detection Success: Did we find the true period?
   - Period Error: How accurate is our measurement?
   - False Positive Rate: Do we detect patterns in noise?

4. TRADE-OFFS:
   - Increasing threshold reduces false positives
   - Increasing bins improves period accuracy
   - Both changes may reduce detection rate for weak signals
""")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader
import seaborn as sns

# Set random seeds for reproducibility
np.random.seed(42)
torch.manual_seed(42)

# Check if GPU is available
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"Using device: {device}")

# ==============================================================================
# 1. Generate Synthetic Radiative Transfer Data
# ==============================================================================

class RadiativeTransferSimulator:
    """
    Simulates radiative transfer through a layered atmosphere
    using the Schwarzschild equation for thermal emission
    """
    def __init__(self, n_layers=50):
        self.n_layers = n_layers
        self.tau_levels = np.linspace(0, 10, n_layers)  # Optical depth levels
        
    def planck_function(self, temperature):
        """
        Simplified Planck function (normalized)
        B(T) ∝ T^4 (Stefan-Boltzmann approximation)
        """
        return (temperature / 300.0) ** 4
    
    def solve_rte(self, temperature_profile):
        """
        Solve the radiative transfer equation using discrete ordinates
        dI/dτ = I - B(τ)
        
        Upward intensity at TOA (top of atmosphere)
        """
        n_layers = len(temperature_profile)
        dtau = self.tau_levels[1] - self.tau_levels[0]
        
        # Planck function at each level
        B = self.planck_function(temperature_profile)
        
        # Initialize intensity (starting from surface)
        I = B[-1]  # Surface emission
        
        # Integrate upward through atmosphere
        intensities = [I]
        for i in range(n_layers - 2, -1, -1):
            # First-order upward integration
            # I(τ-Δτ) = I(τ)·exp(-Δτ) + B(τ)·(1-exp(-Δτ))
            transmission = np.exp(-dtau)
            I = I * transmission + B[i] * (1 - transmission)
            intensities.append(I)
        
        return intensities[-1]  # TOA radiance
    
    def generate_training_data(self, n_samples=1000):
        """
        Generate synthetic training data with various temperature profiles
        """
        X = []  # Temperature profiles
        y = []  # Outgoing radiance at TOA
        
        for _ in range(n_samples):
            # Generate realistic temperature profile
            # T(τ) = T_surface - lapse_rate * height + perturbations
            T_surface = np.random.uniform(280, 310)  # K
            lapse_rate = np.random.uniform(5, 10)  # K per unit optical depth
            
            # Add atmospheric layers with temperature variation
            temp_profile = T_surface - lapse_rate * self.tau_levels
            
            # Add random perturbations (atmospheric variability)
            perturbations = np.random.normal(0, 5, self.n_layers)
            temp_profile += perturbations
            
            # Ensure physical constraints
            temp_profile = np.clip(temp_profile, 200, 320)
            
            # Solve RTE for this profile
            radiance = self.solve_rte(temp_profile)
            
            X.append(temp_profile)
            y.append(radiance)
        
        return np.array(X), np.array(y)

# ==============================================================================
# 2. Create PyTorch Dataset and DataLoader
# ==============================================================================

class RadiativeTransferDataset(Dataset):
    """PyTorch Dataset for radiative transfer data"""
    def __init__(self, X, y):
        self.X = torch.FloatTensor(X)
        self.y = torch.FloatTensor(y).unsqueeze(1)
    
    def __len__(self):
        return len(self.X)
    
    def __getitem__(self, idx):
        return self.X[idx], self.y[idx]

# ==============================================================================
# 3. Define Neural Network Architecture
# ==============================================================================

class RadiativeTransferNN(nn.Module):
    """
    Neural Network to approximate radiative transfer equation solution
    Architecture: Deep feedforward network with residual connections
    """
    def __init__(self, input_dim, hidden_dims=[128, 256, 256, 128]):
        super(RadiativeTransferNN, self).__init__()
        
        layers = []
        prev_dim = input_dim
        
        # Build hidden layers
        for hidden_dim in hidden_dims:
            layers.append(nn.Linear(prev_dim, hidden_dim))
            layers.append(nn.BatchNorm1d(hidden_dim))
            layers.append(nn.ReLU())
            layers.append(nn.Dropout(0.2))
            prev_dim = hidden_dim
        
        self.hidden_layers = nn.Sequential(*layers)
        
        # Output layer (single value: TOA radiance)
        self.output_layer = nn.Linear(prev_dim, 1)
        
    def forward(self, x):
        x = self.hidden_layers(x)
        x = self.output_layer(x)
        return x

# ==============================================================================
# 4. Training Loop with Custom Loss Functions
# ==============================================================================

class RTELoss(nn.Module):
    """
    Custom loss function for radiative transfer
    Combines MSE with physics-informed constraints
    """
    def __init__(self, lambda_physics=0.1):
        super(RTELoss, self).__init__()
        self.lambda_physics = lambda_physics
        self.mse = nn.MSELoss()
    
    def forward(self, predictions, targets, inputs):
        # Standard MSE loss
        mse_loss = self.mse(predictions, targets)
        
        # Physics-informed loss: radiance should increase with temperature
        # This is a soft constraint based on Stefan-Boltzmann law
        mean_temp = inputs.mean(dim=1, keepdim=True)
        expected_trend = (mean_temp / 300.0) ** 4
        
        physics_loss = torch.mean((predictions / targets - 1.0) ** 2)
        
        total_loss = mse_loss + self.lambda_physics * physics_loss
        
        return total_loss, mse_loss, physics_loss

def train_model(model, train_loader, val_loader, epochs=100, lr=0.001):
    """
    Training loop with validation and loss tracking
    """
    model = model.to(device)
    optimizer = optim.Adam(model.parameters(), lr=lr, weight_decay=1e-5)
    scheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='min', 
                                                       factor=0.5, patience=10)
    
    criterion = RTELoss(lambda_physics=0.1)
    
    # Track losses
    train_losses = []
    val_losses = []
    mse_losses = []
    physics_losses = []
    
    best_val_loss = float('inf')
    
    for epoch in range(epochs):
        # Training phase
        model.train()
        train_loss_epoch = 0
        mse_loss_epoch = 0
        physics_loss_epoch = 0
        
        for batch_X, batch_y in train_loader:
            batch_X, batch_y = batch_X.to(device), batch_y.to(device)
            
            optimizer.zero_grad()
            
            # Forward pass
            predictions = model(batch_X)
            
            # Compute loss
            total_loss, mse_loss, physics_loss = criterion(predictions, batch_y, batch_X)
            
            # Backward pass
            total_loss.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
            optimizer.step()
            
            train_loss_epoch += total_loss.item()
            mse_loss_epoch += mse_loss.item()
            physics_loss_epoch += physics_loss.item()
        
        train_loss_epoch /= len(train_loader)
        mse_loss_epoch /= len(train_loader)
        physics_loss_epoch /= len(train_loader)
        
        # Validation phase
        model.eval()
        val_loss_epoch = 0
        
        with torch.no_grad():
            for batch_X, batch_y in val_loader:
                batch_X, batch_y = batch_X.to(device), batch_y.to(device)
                predictions = model(batch_X)
                total_loss, _, _ = criterion(predictions, batch_y, batch_X)
                val_loss_epoch += total_loss.item()
        
        val_loss_epoch /= len(val_loader)
        
        # Learning rate scheduling
        scheduler.step(val_loss_epoch)
        
        # Save best model
        if val_loss_epoch < best_val_loss:
            best_val_loss = val_loss_epoch
            torch.save(model.state_dict(), 'best_rte_model.pth')
        
        # Track losses
        train_losses.append(train_loss_epoch)
        val_losses.append(val_loss_epoch)
        mse_losses.append(mse_loss_epoch)
        physics_losses.append(physics_loss_epoch)
        
        # Print progress
        if (epoch + 1) % 10 == 0:
            print(f"Epoch [{epoch+1}/{epochs}] - "
                  f"Train Loss: {train_loss_epoch:.6f}, "
                  f"Val Loss: {val_loss_epoch:.6f}, "
                  f"MSE: {mse_loss_epoch:.6f}, "
                  f"Physics: {physics_loss_epoch:.6f}")
    
    return train_losses, val_losses, mse_losses, physics_losses

# ==============================================================================
# 5. Main Execution
# ==============================================================================

print("=" * 80)
print("RADIATIVE TRANSFER MODEL OPTIMIZATION WITH NEURAL NETWORKS")
print("=" * 80)
print()

# Generate data
print("Step 1: Generating synthetic radiative transfer data...")
simulator = RadiativeTransferSimulator(n_layers=50)
X_train, y_train = simulator.generate_training_data(n_samples=5000)
X_val, y_val = simulator.generate_training_data(n_samples=1000)
X_test, y_test = simulator.generate_training_data(n_samples=500)

print(f"Training samples: {len(X_train)}")
print(f"Validation samples: {len(X_val)}")
print(f"Test samples: {len(X_test)}")
print(f"Input dimension: {X_train.shape[1]} (atmospheric layers)")
print()

# Create datasets and dataloaders
train_dataset = RadiativeTransferDataset(X_train, y_train)
val_dataset = RadiativeTransferDataset(X_val, y_val)
test_dataset = RadiativeTransferDataset(X_test, y_test)

train_loader = DataLoader(train_dataset, batch_size=64, shuffle=True)
val_loader = DataLoader(val_dataset, batch_size=64, shuffle=False)
test_loader = DataLoader(test_dataset, batch_size=64, shuffle=False)

# Create model
print("Step 2: Initializing neural network model...")
model = RadiativeTransferNN(input_dim=50, hidden_dims=[128, 256, 256, 128])
total_params = sum(p.numel() for p in model.parameters())
print(f"Total parameters: {total_params:,}")
print()

# Train model
print("Step 3: Training the model...")
print("-" * 80)
train_losses, val_losses, mse_losses, physics_losses = train_model(
    model, train_loader, val_loader, epochs=100, lr=0.001
)
print("-" * 80)
print()

# Load best model
model.load_state_dict(torch.load('best_rte_model.pth'))

# Evaluate on test set
print("Step 4: Evaluating on test set...")
model.eval()
test_predictions = []
test_targets = []

with torch.no_grad():
    for batch_X, batch_y in test_loader:
        batch_X = batch_X.to(device)
        predictions = model(batch_X)
        test_predictions.extend(predictions.cpu().numpy())
        test_targets.extend(batch_y.numpy())

test_predictions = np.array(test_predictions).flatten()
test_targets = np.array(test_targets).flatten()

# Calculate metrics
mse = np.mean((test_predictions - test_targets) ** 2)
rmse = np.sqrt(mse)
mae = np.mean(np.abs(test_predictions - test_targets))
r2 = 1 - (np.sum((test_targets - test_predictions) ** 2) / 
          np.sum((test_targets - np.mean(test_targets)) ** 2))

print(f"Test MSE: {mse:.6f}")
print(f"Test RMSE: {rmse:.6f}")
print(f"Test MAE: {mae:.6f}")
print(f"Test R²: {r2:.6f}")
print()

# ==============================================================================
# 6. Visualization
# ==============================================================================

print("Step 5: Generating visualizations...")

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

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

# Plot 1: Training History
ax1 = plt.subplot(2, 3, 1)
epochs_range = range(1, len(train_losses) + 1)
ax1.plot(epochs_range, train_losses, label='Training Loss', linewidth=2)
ax1.plot(epochs_range, val_losses, label='Validation Loss', linewidth=2)
ax1.set_xlabel('Epoch', fontsize=12)
ax1.set_ylabel('Loss', fontsize=12)
ax1.set_title('Training and Validation Loss', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Plot 2: Loss Components
ax2 = plt.subplot(2, 3, 2)
ax2.plot(epochs_range, mse_losses, label='MSE Loss', linewidth=2)
ax2.plot(epochs_range, physics_losses, label='Physics Loss', linewidth=2)
ax2.set_xlabel('Epoch', fontsize=12)
ax2.set_ylabel('Loss', fontsize=12)
ax2.set_title('Loss Components', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# Plot 3: Predictions vs True Values
ax3 = plt.subplot(2, 3, 3)
ax3.scatter(test_targets, test_predictions, alpha=0.5, s=20)
min_val = min(test_targets.min(), test_predictions.min())
max_val = max(test_targets.max(), test_predictions.max())
ax3.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Perfect Prediction')
ax3.set_xlabel('True Radiance', fontsize=12)
ax3.set_ylabel('Predicted Radiance', fontsize=12)
ax3.set_title(f'Predictions vs True Values (R²={r2:.4f})', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# Plot 4: Residual Distribution
ax4 = plt.subplot(2, 3, 4)
residuals = test_predictions - test_targets
ax4.hist(residuals, bins=50, edgecolor='black', alpha=0.7)
ax4.axvline(x=0, color='r', linestyle='--', linewidth=2)
ax4.set_xlabel('Residual (Predicted - True)', fontsize=12)
ax4.set_ylabel('Frequency', fontsize=12)
ax4.set_title(f'Residual Distribution (MAE={mae:.4f})', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

# Plot 5: Example Temperature Profiles and Predictions
ax5 = plt.subplot(2, 3, 5)
n_examples = 3
for i in range(n_examples):
    idx = np.random.randint(0, len(X_test))
    temp_profile = X_test[idx]
    ax5.plot(temp_profile, simulator.tau_levels, label=f'Profile {i+1}', linewidth=2)

ax5.set_xlabel('Temperature (K)', fontsize=12)
ax5.set_ylabel('Optical Depth (τ)', fontsize=12)
ax5.set_title('Example Temperature Profiles', fontsize=14, fontweight='bold')
ax5.invert_yaxis()
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)

# Plot 6: Relative Error Distribution
ax6 = plt.subplot(2, 3, 6)
relative_errors = np.abs((test_predictions - test_targets) / test_targets) * 100
ax6.hist(relative_errors, bins=50, edgecolor='black', alpha=0.7)
ax6.set_xlabel('Relative Error (%)', fontsize=12)
ax6.set_ylabel('Frequency', fontsize=12)
ax6.set_title(f'Relative Error Distribution (Mean: {np.mean(relative_errors):.2f}%)', 
              fontsize=14, fontweight='bold')
ax6.grid(True, alpha=0.3)

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

print("Visualization complete!")
print("=" * 80)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.cluster import DBSCAN
from sklearn.metrics import silhouette_score, davies_bouldin_score
from sklearn.preprocessing import StandardScaler
import seaborn as sns
from itertools import product

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

# Generate synthetic galaxy distribution data
def generate_galaxy_clusters(n_clusters=5, n_galaxies_per_cluster=50, n_noise=30):
    """
    Generate synthetic 3D galaxy distribution data with multiple clusters
    
    Parameters:
    - n_clusters: Number of galaxy clusters
    - n_galaxies_per_cluster: Number of galaxies in each cluster
    - n_noise: Number of background/noise galaxies
    """
    galaxies = []
    labels_true = []
    
    # Generate clustered galaxies
    for i in range(n_clusters):
        # Random cluster center in 3D space (in Mpc - Megaparsecs)
        center = np.random.uniform(-100, 100, 3)
        
        # Generate galaxies around the center with some dispersion
        cluster_galaxies = np.random.randn(n_galaxies_per_cluster, 3) * 5 + center
        galaxies.append(cluster_galaxies)
        labels_true.extend([i] * n_galaxies_per_cluster)
    
    # Generate noise/background galaxies
    noise_galaxies = np.random.uniform(-120, 120, (n_noise, 3))
    galaxies.append(noise_galaxies)
    labels_true.extend([-1] * n_noise)
    
    # Combine all galaxies
    galaxies = np.vstack(galaxies)
    labels_true = np.array(labels_true)
    
    return galaxies, labels_true

# Generate the data
print("=" * 60)
print("GALAXY CLUSTER IDENTIFICATION SYSTEM")
print("=" * 60)
print("\nGenerating synthetic galaxy distribution data...")
galaxies, true_labels = generate_galaxy_clusters(
    n_clusters=5, 
    n_galaxies_per_cluster=50, 
    n_noise=30
)
print(f"Total galaxies: {len(galaxies)}")
print(f"True number of clusters: {len(np.unique(true_labels[true_labels >= 0]))}")
print(f"Noise galaxies: {np.sum(true_labels == -1)}")

# Normalize the data for better clustering
scaler = StandardScaler()
galaxies_scaled = scaler.fit_transform(galaxies)

# Parameter optimization: Grid search over epsilon and min_samples
print("\n" + "=" * 60)
print("OPTIMIZING DBSCAN PARAMETERS")
print("=" * 60)

epsilon_values = np.linspace(0.2, 1.5, 15)
min_samples_values = [3, 5, 7, 10, 15, 20]

results = []

print("\nTesting parameter combinations...")
for eps, min_samp in product(epsilon_values, min_samples_values):
    # Apply DBSCAN
    dbscan = DBSCAN(eps=eps, min_samples=min_samp)
    labels = dbscan.fit_predict(galaxies_scaled)
    
    # Count clusters (excluding noise labeled as -1)
    n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
    n_noise = list(labels).count(-1)
    
    # Calculate metrics only if we have at least 2 clusters
    if n_clusters >= 2:
        # Filter out noise points for silhouette score
        mask = labels != -1
        if np.sum(mask) > 0:
            try:
                silhouette = silhouette_score(galaxies_scaled[mask], labels[mask])
                davies_bouldin = davies_bouldin_score(galaxies_scaled[mask], labels[mask])
            except:
                silhouette = -1
                davies_bouldin = 999
        else:
            silhouette = -1
            davies_bouldin = 999
    else:
        silhouette = -1
        davies_bouldin = 999
    
    results.append({
        'eps': eps,
        'min_samples': min_samp,
        'n_clusters': n_clusters,
        'n_noise': n_noise,
        'silhouette': silhouette,
        'davies_bouldin': davies_bouldin,
        'labels': labels
    })

# Find best parameters based on silhouette score
valid_results = [r for r in results if r['silhouette'] > -1]
if valid_results:
    best_result = max(valid_results, key=lambda x: x['silhouette'])
    print(f"\n✓ Best parameters found:")
    print(f"  ε (epsilon) = {best_result['eps']:.3f}")
    print(f"  MinPts (min_samples) = {best_result['min_samples']}")
    print(f"  Number of clusters identified: {best_result['n_clusters']}")
    print(f"  Number of noise points: {best_result['n_noise']}")
    print(f"  Silhouette Score: {best_result['silhouette']:.4f}")
    print(f"  Davies-Bouldin Index: {best_result['davies_bouldin']:.4f}")
else:
    print("No valid clustering found!")
    best_result = results[0]

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

# 1. 3D scatter plot of original galaxy distribution
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
scatter1 = ax1.scatter(galaxies[:, 0], galaxies[:, 1], galaxies[:, 2], 
                       c=true_labels, cmap='tab10', s=30, alpha=0.6, 
                       edgecolors='black', linewidth=0.5)
ax1.set_xlabel('X (Mpc)', fontsize=10)
ax1.set_ylabel('Y (Mpc)', fontsize=10)
ax1.set_zlabel('Z (Mpc)', fontsize=10)
ax1.set_title('True Galaxy Distribution\n(5 clusters + noise)', fontsize=12, fontweight='bold')
plt.colorbar(scatter1, ax=ax1, label='True Cluster ID', pad=0.1)

# 2. 3D scatter plot with DBSCAN clustering results
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
scatter2 = ax2.scatter(galaxies[:, 0], galaxies[:, 1], galaxies[:, 2], 
                       c=best_result['labels'], cmap='tab10', s=30, alpha=0.6,
                       edgecolors='black', linewidth=0.5)
ax2.set_xlabel('X (Mpc)', fontsize=10)
ax2.set_ylabel('Y (Mpc)', fontsize=10)
ax2.set_zlabel('Z (Mpc)', fontsize=10)
ax2.set_title(f'DBSCAN Clustering Results\n(ε={best_result["eps"]:.3f}, MinPts={best_result["min_samples"]})', 
              fontsize=12, fontweight='bold')
plt.colorbar(scatter2, ax=ax2, label='Cluster ID', pad=0.1)

# 3. Parameter optimization heatmap (Silhouette Score)
ax3 = fig.add_subplot(2, 3, 3)
silhouette_matrix = np.full((len(min_samples_values), len(epsilon_values)), -1.0)
for r in results:
    i = min_samples_values.index(r['min_samples'])
    j = np.argmin(np.abs(epsilon_values - r['eps']))
    silhouette_matrix[i, j] = r['silhouette']

im = ax3.imshow(silhouette_matrix, cmap='RdYlGn', aspect='auto', vmin=-1, vmax=1)
ax3.set_xticks(range(len(epsilon_values)))
ax3.set_yticks(range(len(min_samples_values)))
ax3.set_xticklabels([f'{e:.2f}' for e in epsilon_values], rotation=45, ha='right', fontsize=8)
ax3.set_yticklabels(min_samples_values, fontsize=9)
ax3.set_xlabel('ε (epsilon)', fontsize=11, fontweight='bold')
ax3.set_ylabel('MinPts (min_samples)', fontsize=11, fontweight='bold')
ax3.set_title('Silhouette Score Heatmap\n(Higher is Better)', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax3, label='Silhouette Score')

# Mark the best parameter
best_i = min_samples_values.index(best_result['min_samples'])
best_j = np.argmin(np.abs(epsilon_values - best_result['eps']))
ax3.plot(best_j, best_i, 'b*', markersize=20, markeredgecolor='white', markeredgewidth=2)

# 4. Number of clusters vs parameters
ax4 = fig.add_subplot(2, 3, 4)
cluster_matrix = np.zeros((len(min_samples_values), len(epsilon_values)))
for r in results:
    i = min_samples_values.index(r['min_samples'])
    j = np.argmin(np.abs(epsilon_values - r['eps']))
    cluster_matrix[i, j] = r['n_clusters']

im2 = ax4.imshow(cluster_matrix, cmap='viridis', aspect='auto')
ax4.set_xticks(range(len(epsilon_values)))
ax4.set_yticks(range(len(min_samples_values)))
ax4.set_xticklabels([f'{e:.2f}' for e in epsilon_values], rotation=45, ha='right', fontsize=8)
ax4.set_yticklabels(min_samples_values, fontsize=9)
ax4.set_xlabel('ε (epsilon)', fontsize=11, fontweight='bold')
ax4.set_ylabel('MinPts (min_samples)', fontsize=11, fontweight='bold')
ax4.set_title('Number of Clusters Identified', fontsize=12, fontweight='bold')
plt.colorbar(im2, ax=ax4, label='# Clusters')
ax4.plot(best_j, best_i, 'r*', markersize=20, markeredgecolor='white', markeredgewidth=2)

# 5. Davies-Bouldin Index heatmap
ax5 = fig.add_subplot(2, 3, 5)
db_matrix = np.full((len(min_samples_values), len(epsilon_values)), 999.0)
for r in results:
    i = min_samples_values.index(r['min_samples'])
    j = np.argmin(np.abs(epsilon_values - r['eps']))
    if r['davies_bouldin'] < 999:
        db_matrix[i, j] = r['davies_bouldin']

# Mask invalid values
db_matrix_masked = np.ma.masked_where(db_matrix >= 999, db_matrix)
im3 = ax5.imshow(db_matrix_masked, cmap='RdYlGn_r', aspect='auto', vmin=0, vmax=3)
ax5.set_xticks(range(len(epsilon_values)))
ax5.set_yticks(range(len(min_samples_values)))
ax5.set_xticklabels([f'{e:.2f}' for e in epsilon_values], rotation=45, ha='right', fontsize=8)
ax5.set_yticklabels(min_samples_values, fontsize=9)
ax5.set_xlabel('ε (epsilon)', fontsize=11, fontweight='bold')
ax5.set_ylabel('MinPts (min_samples)', fontsize=11, fontweight='bold')
ax5.set_title('Davies-Bouldin Index\n(Lower is Better)', fontsize=12, fontweight='bold')
plt.colorbar(im3, ax=ax5, label='DB Index')
ax5.plot(best_j, best_i, 'b*', markersize=20, markeredgecolor='white', markeredgewidth=2)

# 6. Cluster size distribution
ax6 = fig.add_subplot(2, 3, 6)
unique_labels = set(best_result['labels'])
cluster_sizes = [np.sum(best_result['labels'] == label) for label in unique_labels if label != -1]
noise_count = np.sum(best_result['labels'] == -1)

bars = ax6.bar(range(len(cluster_sizes)), cluster_sizes, color='steelblue', 
               edgecolor='black', linewidth=1.5, alpha=0.7)
if noise_count > 0:
    ax6.bar(len(cluster_sizes), noise_count, color='red', 
            edgecolor='black', linewidth=1.5, alpha=0.7, label='Noise')

ax6.set_xlabel('Cluster ID', fontsize=11, fontweight='bold')
ax6.set_ylabel('Number of Galaxies', fontsize=11, fontweight='bold')
ax6.set_title('Galaxy Count per Cluster', fontsize=12, fontweight='bold')
ax6.grid(axis='y', alpha=0.3, linestyle='--')
if noise_count > 0:
    ax6.legend(fontsize=10)

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

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

# Additional analysis: Performance metrics across all tested parameters
print("\n" + "=" * 60)
print("PARAMETER SENSITIVITY ANALYSIS")
print("=" * 60)

# Find top 5 parameter combinations
top_5 = sorted(valid_results, key=lambda x: x['silhouette'], reverse=True)[:5]
print("\nTop 5 Parameter Combinations:")
print("-" * 60)
for i, result in enumerate(top_5, 1):
    print(f"{i}. ε={result['eps']:.3f}, MinPts={result['min_samples']}")
    print(f"   Clusters: {result['n_clusters']}, Silhouette: {result['silhouette']:.4f}, DB: {result['davies_bouldin']:.4f}")

# Statistical summary
print("\n" + "=" * 60)
print("STATISTICAL SUMMARY")
print("=" * 60)
print(f"Total parameter combinations tested: {len(results)}")
print(f"Valid clustering solutions: {len(valid_results)}")
print(f"Average number of clusters found: {np.mean([r['n_clusters'] for r in valid_results]):.2f}")
print(f"Silhouette score range: [{min([r['silhouette'] for r in valid_results]):.4f}, {max([r['silhouette'] for r in valid_results]):.4f}]")

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

def generate_galaxy_clusters(n_clusters=5, n_galaxies_per_cluster=50, n_noise=30)

scaler = StandardScaler()
galaxies_scaled = scaler.fit_transform(galaxies)