""" ============================================================= MMA-CCT: Multi-Messenger Astronomy - Conditional Collapse Theory Neutrino Detection System based on CCT-ODE Framework Architecture: Transformer with Cross-Channel Sensitivity Tracking Training: Particle Statistics (Fermi-Dirac) + Sensitivity Amplification Detection: Rare Event Hypersensitivity + Entropic Collapse Author: CCT-ODE Framework Version: 1.0 ============================================================= """ import torch import torch.nn as nn import torch.nn.functional as F import torch.optim as optim from torch.utils.data import Dataset, DataLoader import math import numpy as np from typing import Dict, List, Tuple, Optional from dataclasses import dataclass import warnings # ============================================================================= # SECTION 1: CORE MATHEMATICAL CONSTANTS # ============================================================================= # Physical constants (normalized for computational efficiency) C_LIGHT = 299792458.0 # m/s (speed of light) EV_TO_JOULES = 1.60218e-19 MEV_TO_EV = 1.0e6 GEV_TO_EV = 1.0e9 # Neutrino detection parameters (IceCube-style) NEUTRINO_ENERGY_RANGE = (0.1, 100.0) # GeV NEUTRINO_FLUX_TOTAL = 1.0e3 # per km² per year (atmospheric) BACKGROUND_RATE = 1.0e9 # background events per km² per year SIGNAL_BRANCHING_RATIO = 1e-5 # rare signal events # Detection geometry (simplified cylinder detector) DETECTOR_VOLUME = 1.0 # km³ (normalized) PHOTON_SPEED = 0.8 * C_LIGHT # effective speed in medium @dataclass class PhysicsConstants: """Container for physics simulation constants""" c: float = C_LIGHT mev_to_ev: float = MEV_TO_EV gev_to_ev: float = GEV_TO_EV ice_density: float = 0.92 # g/cm³ refractive_index: float = 1.33 # ice cherenkov_angle: float = math.cos(41 * math.pi / 180) # ~41 degrees absorption_length: float = 100.0 # meters scattering_length: float = 25.0 # meters PHYSICS = PhysicsConstants() # ============================================================================= # SECTION 2: NEUTRINO DATA SIMULATION (REALISTIC) # ============================================================================= class NeutrinoSimulator: """ Realistic neutrino detection simulator based on IceCube-like detector. Simulates: - Cherenkov light emission from charged particles - Photon propagation and arrival times - Background noise (cosmic rays, atmospheric muons) - Flavor oscillation effects - Energy-dependent event signatures """ # Neutrino flavors FLAVOR_ELECTRON = 0 FLAVOR_MUON = 1 FLAVOR_TAU = 2 # Interaction types INTERACTION_CC = 0 # Charged current INTERACTION_NC = 1 # Neutral current def __init__(self, device='cpu'): self.device = device self.physics = PHYSICS def simulate_event(self, is_signal: bool = False, flavor: Optional[int] = None, energy: Optional[float] = None, direction: Optional[torch.Tensor] = None) -> Dict[str, torch.Tensor]: """ Simulate a single neutrino or background event. Returns: Dictionary with: - 'hits': [num_hits, 4] tensor (x, y, z, t) - 'features': [num_features] tensor (energy, direction, time, etc.) - 'label': int (0=background, 1=signal) - 'metadata': dict with event properties """ # Generate signal or background event if is_signal: event = self._simulate_neutrino(flavor, energy, direction) else: event = self._simulate_background() return event def _simulate_neutrino(self, flavor: Optional[int], energy: Optional[float], direction: Optional[torch.Tensor]) -> Dict[str, torch.Tensor]: """Simulate a real neutrino interaction""" # Randomize flavor if not specified if flavor is None: flavor = torch.randint(0, 3, (1,)).item() # Randomize energy (power law spectrum: dN/dE ∝ E^-2) if energy is None: energy = torch.rand(1).item() energy = 10 ** (energy * 3) # 1 GeV to 10 TeV range energy = min(energy, 100.0) # Cap at 100 GeV # Randomize direction (uniform on sphere) if direction is None: theta = torch.rand(1).item() * math.pi phi = torch.rand(1).item() * 2 * math.pi direction = torch.tensor([ math.sin(theta) * math.cos(phi), math.sin(theta) * math.sin(phi), math.cos(theta) ]) # Determine interaction type interaction_type = torch.randint(0, 2, (1,)).item() # Generate Cherenkov light hits num_hits = self._calculate_hits(energy, flavor, interaction_type) hits = self._generate_cherenkov_hits(num_hits, energy, direction, flavor) # Extract features features = self._extract_features(hits, energy, direction, flavor, interaction_type) return { 'hits': hits, 'features': features, 'label': torch.tensor(1), 'metadata': { 'flavor': flavor, 'energy': energy, 'direction': direction, 'interaction_type': interaction_type, 'num_hits': num_hits } } def _simulate_background(self) -> Dict[str, torch.Tensor]: """Simulate background event (muons, noise)""" # Background is mostly atmospheric muons # Has different hit pattern than neutrinos num_hits = torch.randint(5, 50, (1,)).item() # Random positions in detector volume x = (torch.rand(num_hits) - 0.5) * 1000 # meters y = (torch.rand(num_hits) - 0.5) * 1000 z = (torch.rand(num_hits) - 0.5) * 1000 # Times spread along a track (muons are track-like) t = torch.rand(num_hits) * 1000 # nanoseconds hits = torch.stack([x, y, z, t], dim=1) # Low reconstructed energy (background) energy = torch.rand(1).item() * 0.5 # Random direction theta = torch.rand(1).item() * math.pi phi = torch.rand(1).item() * 2 * math.pi direction = torch.tensor([ math.sin(theta) * math.cos(phi), math.sin(theta) * math.sin(phi), math.cos(theta) ]) features = torch.zeros(12) features[0] = energy # reconstructed energy features[1:4] = direction # direction features[4] = num_hits / 100.0 # normalized hit count return { 'hits': hits, 'features': features, 'label': torch.tensor(0), 'metadata': { 'flavor': -1, 'energy': energy, 'type': 'background' } } def _calculate_hits(self, energy: float, flavor: int, interaction_type: int) -> int: """Calculate number of Cherenkov photons detected""" # Base number of Cherenkov photons # E_lepton ~ 0.75 * E_nu for CC interactions E_lepton = 0.75 * energy if interaction_type == 0 else 0.0 if E_lepton > 0.5: # GeV threshold for Cherenkov # Cherenkov photon yield dNdx = 350 # photons per meter per MeV (at ~41 deg) track_length = min(E_lepton * 5, 1000) # meters (rough estimate) num_photons = dNdx * E_lepton * 1000 * track_length * 1e-6 # per MeV else: num_photons = 0 # Detection efficiency (~1%) efficiency = 0.01 # Poisson sampling detected = np.random.poisson(num_photons * efficiency) return max(detected, 1) # At least 1 hit def _generate_cherenkov_hits(self, num_hits: int, energy: float, direction: torch.Tensor, flavor: int) -> torch.Tensor: """Generate realistic Cherenkov ring hits""" hits = [] # Determine track characteristics based on flavor if flavor == self.FLAVOR_MUON: # Long track, diffuse pattern track_length = min(energy * 10, 1000) ring_radius = 50 + energy * 5 elif flavor == self.FLAVOR_ELECTRON: # Electromagnetic shower, more compact track_length = min(energy * 2, 200) ring_radius = 30 + energy * 2 else: # TAU # Short track with decay track_length = min(energy * 3, 300) ring_radius = 40 + energy * 3 # Generate hits along Cherenkov cone for i in range(num_hits): # Position along track t_track = torch.rand(1).item() * track_length # Position on cone surface phi = torch.rand(1).item() * 2 * math.pi r_offset = ring_radius * (0.5 + torch.rand(1).item()) # Create position relative to track base_pos = direction * t_track perp_x = torch.tensor([direction[1], -direction[0], 0.0]) perp_x = perp_x / (torch.norm(perp_x) + 1e-6) perp_y = torch.linalg.cross(direction, perp_x) pos = base_pos + r_offset * (perp_x * math.cos(phi) + perp_y * math.sin(phi)) # Time of arrival (propagation delay) dist = torch.norm(pos) time = dist / PHYSICS.cherenkov_angle / (PHYSICS.c * 1e-9) # nanoseconds time += torch.rand(1).item() * 50 # timing jitter hits.append([pos[0], pos[1], pos[2], time]) hits_tensor = torch.tensor(hits, dtype=torch.float32) # Normalize positions to [-1, 1] range if hits_tensor.shape[0] > 0: hits_tensor[:, :3] = hits_tensor[:, :3] / 500.0 # 500m normalization return hits_tensor def _extract_features(self, hits: torch.Tensor, energy: float, direction: torch.Tensor, flavor: int, interaction_type: int) -> torch.Tensor: """Extract physics features from hit pattern""" features = torch.zeros(12) # Basic features features[0] = math.log10(energy + 1e-6) / 3.0 # log-normalized energy features[1:4] = direction # direction cosines features[4] = hits.shape[0] / 100.0 # normalized hit count if hits.shape[0] > 3: # Spatial statistics centroid = hits[:, :3].mean(dim=0) features[5:8] = centroid # Temporal spread time_spread = hits[:, 3].std() features[8] = time_spread / 100.0 # PCA-based track length estimate centered = hits[:, :3] - centroid cov = centered.T @ centered / max(centered.shape[0], 1) eigenvalues = torch.linalg.eigvalsh(cov) track_score = eigenvalues[-1] / (eigenvalues.sum() + 1e-6) features[9] = track_score # Angular estimate from hit distribution features[10] = energy / 100.0 # energy proxy for track features[11] = interaction_type / 1.0 # interaction type else: # Not enough hits features[5:8] = torch.zeros(3) features[8:12] = torch.zeros(4) return features class NeutrinoDataset(Dataset): """Dataset for neutrino detection training""" def __init__(self, num_samples=10000, signal_ratio=0.1, device='cpu'): self.num_samples = num_samples self.signal_ratio = signal_ratio self.device = device self.simulator = NeutrinoSimulator(device) # Pre-generate all data self.data = [] self._generate_data() def _generate_data(self): """Generate dataset with realistic event distribution""" num_signals = int(self.num_samples * self.signal_ratio) num_background = self.num_samples - num_signals # Generate background events for _ in range(num_background): event = self.simulator.simulate_event(is_signal=False) self.data.append(event) # Generate signal events with varied properties for _ in range(num_signals): flavor = torch.randint(0, 3, (1,)).item() energy = 10 ** (torch.rand(1).item() * 3) # 1 GeV to 10 TeV energy = min(energy, 100.0) theta = torch.rand(1).item() * math.pi phi = torch.rand(1).item() * 2 * math.pi direction = torch.tensor([ math.sin(theta) * math.cos(phi), math.sin(theta) * math.sin(phi), math.cos(theta) ]) event = self.simulator.simulate_event( is_signal=True, flavor=flavor, #energy=energy.item(), energy=energy, direction=direction ) self.data.append(event) # Shuffle indices = torch.randperm(len(self.data)).tolist() self.data = [self.data[i] for i in indices] def __len__(self): return self.num_samples def __getitem__(self, idx): event = self.data[idx] # Pad hits to max length max_hits = 200 hits = event['hits'] if hits.shape[0] < max_hits: padding = torch.zeros(max_hits - hits.shape[0], 4) hits = torch.cat([hits, padding], dim=0) else: hits = hits[:max_hits] return { 'hits': hits, 'features': event['features'], 'label': event['label'] } # ============================================================================= # SECTION 3: CCT-ODE TRANSFORMER ENCODER # ============================================================================= class PositionalEncoding(nn.Module): """Positional encoding for transformer""" def __init__(self, d_model: int, max_len: int = 500): super().__init__() pe = torch.zeros(max_len, d_model) position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1) div_term = torch.exp(torch.arange(0, d_model, 2).float() * (-math.log(10000.0) / d_model)) pe[:, 0::2] = torch.sin(position * div_term) pe[:, 1::2] = torch.cos(position * div_term) self.register_buffer('pe', pe.unsqueeze(0)) def forward(self, x): return x + self.pe[:, :x.size(1)] class MultiHeadAttention(nn.Module): """ Multi-head attention with sensitivity tracking. Implements the CCT-ODE collapse operator. """ def __init__(self, d_model: int, n_heads: int, dropout: float = 0.1): super().__init__() assert d_model % n_heads == 0 self.d_model = d_model self.n_heads = n_heads self.d_k = d_model // n_heads self.W_q = nn.Linear(d_model, d_model) self.W_k = nn.Linear(d_model, d_model) self.W_v = nn.Linear(d_model, d_model) self.W_o = nn.Linear(d_model, d_model) self.dropout = nn.Dropout(dropout) self.layer_norm = nn.LayerNorm(d_model) # Sensitivity tracking (CCT-ODE) self.sensitivity_history = [] self.collapse_history = [] def forward(self, query, key, value, mask=None, return_sensitivity=False): """ Forward pass with optional sensitivity tracking. Returns: - output: [batch, seq, d_model] - sensitivity: [batch, n_heads, seq, seq] (optional) """ batch_size = query.size(0) seq_len = query.size(1) # Linear projections Q = self.W_q(query) K = self.W_k(key) V = self.W_v(value) # Reshape for multi-head Q = Q.view(batch_size, seq_len, self.n_heads, self.d_k).transpose(1, 2) K = K.view(batch_size, -1, self.n_heads, self.d_k).transpose(1, 2) V = V.view(batch_size, -1, self.n_heads, self.d_k).transpose(1, 2) # Attention scores scores = torch.matmul(Q, K.transpose(-2, -1)) / math.sqrt(self.d_k) if mask is not None: scores = scores.masked_fill(mask == 0, -1e9) # CCT-ODE: Softmax as collapse operator attention_weights = F.softmax(scores, dim=-1) attention_weights = self.dropout(attention_weights) # Compute sensitivity (collapse potential) if return_sensitivity: # Sensitivity = how much each token affects others sensitivity = attention_weights.clone() self.sensitivity_history.append(sensitivity.detach()) else: sensitivity = None # Apply attention context = torch.matmul(attention_weights, V) # Reshape and project context = context.transpose(1, 2).contiguous().view(batch_size, seq_len, self.d_model) output = self.W_o(context) # CCT-ODE: Track entropy reduction (collapse) if return_sensitivity and attention_weights is not None: entropy_before = -torch.sum(attention_weights * torch.log(attention_weights + 1e-10), dim=-1).mean() self.collapse_history.append(entropy_before.item()) return output, sensitivity class CCTODEFeedForward(nn.Module): """Feed-forward network with ODE dynamics""" def __init__(self, d_model: int, d_ff: int = 2048, dropout: float = 0.1): super().__init__() self.linear1 = nn.Linear(d_model, d_ff) self.linear2 = nn.Linear(d_ff, d_model) self.dropout = nn.Dropout(dropout) self.activation = nn.GELU() def forward(self, x): return self.linear2(self.dropout(self.activation(self.linear1(x)))) class CCTODELayer(nn.Module): """ Single CCT-ODE layer combining attention and FFN. Tracks sensitivity matrices per layer. """ def __init__(self, d_model: int, n_heads: int, d_ff: int = 2048, dropout: float = 0.1): super().__init__() self.attention = MultiHeadAttention(d_model, n_heads, dropout) self.feedforward = CCTODEFeedForward(d_model, d_ff, dropout) self.norm1 = nn.LayerNorm(d_model) self.norm2 = nn.LayerNorm(d_model) self.dropout1 = nn.Dropout(dropout) self.dropout2 = nn.Dropout(dropout) # Layer-specific sensitivity matrices (CCT-ODE) self.layer_sensitivity = { 'energy': None, 'time': None, 'direction': None, 'flavor': None } def forward(self, x, mask=None, track_sensitivity=False): """ Forward pass with optional sensitivity tracking. Args: x: [batch, seq, d_model] track_sensitivity: Whether to store sensitivity matrices Returns: - output: [batch, seq, d_model] - sensitivity_dict: dict of sensitivity matrices """ # Self-attention with residual attn_out, sensitivity = self.attention( self.norm1(x), self.norm1(x), self.norm1(x), mask, return_sensitivity=track_sensitivity ) x = x + self.dropout1(attn_out) # Store sensitivity if tracking if track_sensitivity: self.layer_sensitivity['attention'] = sensitivity # Feed-forward with residual ff_out = self.feedforward(self.norm2(x)) x = x + self.dropout2(ff_out) return x, self.layer_sensitivity # ============================================================================= # SECTION 4: NEUTRINO TRANSFORMER (MAIN ARCHITECTURE) # ============================================================================= class NeutrinoTransformer(nn.Module): """ CCT-ODE Transformer for Neutrino Detection. Implements: - Hit-based processing (Transformer on hit coordinates) - Feature-based processing (MLP on physics features) - Sensitivity matrix tracking - Entropic collapse detection - Hypersensitivity for rare events """ def __init__( self, d_model: int = 256, n_heads: int = 8, n_layers: int = 6, d_ff: int = 1024, max_hits: int = 200, n_features: int = 12, n_classes: int = 2, dropout: float = 0.1 ): super().__init__() self.d_model = d_model self.n_heads = n_heads self.n_layers = n_layers self.max_hits = max_hits self.n_classes = n_classes # ============================================ # HIT PROCESSING (Transformer on coordinates) # ============================================ # Input embedding for hits (x, y, z, t -> d_model) self.hit_embedding = nn.Sequential( nn.Linear(4, d_model // 2), nn.LayerNorm(d_model // 2), nn.GELU(), nn.Linear(d_model // 2, d_model), nn.LayerNorm(d_model), nn.GELU() ) self.hit_positional_encoding = PositionalEncoding(d_model, max_hits) self.hit_transformer_layers = nn.ModuleList([ CCTODELayer(d_model, n_heads, d_ff, dropout) for _ in range(n_layers) ]) # ============================================ # FEATURE PROCESSING (MLP on physics features) # ============================================ self.feature_encoder = nn.Sequential( nn.Linear(n_features, d_model), nn.LayerNorm(d_model), nn.GELU(), nn.Linear(d_model, d_model), nn.LayerNorm(d_model), nn.GELU() ) # ============================================ # FUSION LAYER # ============================================ self.fusion = nn.MultiheadAttention( d_model, n_heads, dropout, batch_first=True ) self.fusion_norm = nn.LayerNorm(d_model) # ============================================ # SENSITIVITY TRACKING (CCT-ODE) # ============================================ self.sensitivity_matrices = { 'hit_attention': [], 'feature_attention': [], 'cross_attention': [], 'layer_entropies': [] } # ============================================ # CLASSIFICATION HEAD # ============================================ self.classifier = nn.Sequential( nn.Linear(d_model * 2, d_model), nn.LayerNorm(d_model), nn.GELU(), nn.Dropout(dropout), nn.Linear(d_model, d_model // 2), nn.LayerNorm(d_model // 2), nn.GELU(), nn.Dropout(dropout), nn.Linear(d_model // 2, n_classes) ) # ============================================ # RARE EVENT DETECTION HEAD # ============================================ self.rare_event_detector = nn.Sequential( nn.Linear(d_model, d_model // 4), nn.LayerNorm(d_model // 4), nn.GELU(), nn.Linear(d_model // 4, 1), nn.Sigmoid() ) def forward(self, hits, features, track_sensitivity=True): """ Forward pass of the neutrino transformer. Args: hits: [batch, max_hits, 4] - Hit coordinates features: [batch, n_features] - Physics features track_sensitivity: Whether to compute sensitivity matrices Returns: - logits: [batch, n_classes] - sensitivity_dict: dict of sensitivity matrices - rare_event_prob: [batch, 1] """ batch_size = hits.shape[0] # ============================================ # STEP 1: Process Hits (Transformer) # ============================================ # Create mask for padded hits hit_mask = (hits.abs().sum(dim=-1) > 1e-6).float() # [batch, max_hits] # Embed hits hits_embedded = self.hit_embedding(hits) # [batch, max_hits, d_model] hits_embedded = self.hit_positional_encoding(hits_embedded) # Pass through transformer layers hit_representation = hits_embedded for layer_idx, layer in enumerate(self.hit_transformer_layers): hit_representation, layer_sensitivity = layer( hit_representation, mask=hit_mask.unsqueeze(1).unsqueeze(2), track_sensitivity=track_sensitivity ) if track_sensitivity and layer_sensitivity.get('attention') is not None: self.sensitivity_matrices['hit_attention'].append( layer_sensitivity['attention'].detach() ) # Pool hit representations (attention-weighted) hit_pooled = self._attention_pool(hit_representation, hit_mask) # ============================================ # STEP 2: Process Features (MLP) # ============================================ feature_representation = self.feature_encoder(features) # ============================================ # STEP 3: Fuse Representations # ============================================ # Cross-attention between hits and features fused, cross_attn = self.fusion( query=feature_representation.unsqueeze(1), key=hit_representation, value=hit_representation, key_padding_mask=(hit_mask == 0) ) fused = fused.squeeze(1) # [batch, d_model] fused = self.fusion_norm(fused + feature_representation) if track_sensitivity: self.sensitivity_matrices['cross_attention'].append(cross_attn.detach()) # ============================================ # STEP 4: Classification # ============================================ # Concatenate hit-pooled and fused representations joint_representation = torch.cat([hit_pooled, fused], dim=-1) logits = self.classifier(joint_representation) # ============================================ # STEP 5: Rare Event Detection # ============================================ rare_event_prob = self.rare_event_detector(fused) # ============================================ # STEP 6: Compute Entropy (CCT-ODE) # ============================================ if track_sensitivity: total_entropy = self._compute_entropy() self.sensitivity_matrices['layer_entropies'].append(total_entropy) return { 'logits': logits, 'hit_representation': hit_pooled, 'feature_representation': feature_representation, 'fused_representation': fused, 'rare_event_prob': rare_event_prob, 'sensitivity_matrices': self.sensitivity_matrices if track_sensitivity else None } def _attention_pool(self, x, mask): """ Attention-weighted pooling over sequence dimension. Args: x: [batch, seq, d_model] mask: [batch, seq] - 1 for valid, 0 for padding Returns: pooled: [batch, d_model] """ # Simple attention weights based on mask weights = mask.unsqueeze(-1).float() # [batch, seq, 1] weights = weights / (weights.sum(dim=1, keepdim=True) + 1e-6) pooled = (x * weights).sum(dim=1) # [batch, d_model] return pooled def _compute_entropy(self): """Compute entropy of attention distributions (CCT-ODE collapse metric)""" total_entropy = 0.0 for sensitivity_list in [self.sensitivity_matrices['hit_attention']]: for sensitivity in sensitivity_list: # sensitivity: [batch, heads, seq, seq] prob = F.softmax(sensitivity, dim=-1) entropy = -torch.sum(prob * torch.log(prob + 1e-10), dim=-1) total_entropy += entropy.mean().item() return total_entropy def reset_sensitivity_tracking(self): """Reset sensitivity matrices for new batch""" self.sensitivity_matrices = { 'hit_attention': [], 'feature_attention': [], 'cross_attention': [], 'layer_entropies': [] } # ============================================================================= # SECTION 5: CCT-ODE LOSS FUNCTIONS # ============================================================================= class CCTODELoss(nn.Module): """ Conditional Collapse Theory loss function. Combines: - Cross-entropy for classification - Sensitivity amplification for rare events - Fermi-Dirac regularization for particle statistics - Entropy collapse bonus """ def __init__(self, rare_event_weight=5.0, sensitivity_weight=0.1, fermi_temp=1.0): super().__init__() self.rare_event_weight = rare_event_weight self.sensitivity_weight = sensitivity_weight self.fermi_temp = fermi_temp self.ce_loss = nn.CrossEntropyLoss() def forward(self, outputs, targets, sensitivity_matrices=None): """ Compute CCT-ODE loss. Args: outputs: dict from neutrino transformer targets: [batch] tensor of class labels sensitivity_matrices: optional dict of sensitivity matrices Returns: total_loss: scalar loss_dict: dict of individual loss components """ # Standard cross-entropy ce = self.ce_loss(outputs['logits'], targets) # Rare event bonus (encourage high rare_event_prob for signal) signal_mask = (targets == 1).float() rare_bonus = -self.rare_event_weight * ( outputs['rare_event_prob'].squeeze() * signal_mask ).mean() # Sensitivity amplification loss sensitivity_loss = 0.0 if sensitivity_matrices is not None: sensitivity_loss = self._compute_sensitivity_loss(sensitivity_matrices) # Fermi-Dirac regularization (encourage weight distribution following physics) fermi_loss = self._compute_fermi_regularization(outputs) total_loss = ce + rare_bonus + self.sensitivity_weight * sensitivity_loss + 0.01 * fermi_loss loss_dict = { 'cross_entropy': ce.item(), 'rare_event_bonus': rare_bonus.item(), 'sensitivity_loss': sensitivity_loss.item() if isinstance(sensitivity_loss, torch.Tensor) else sensitivity_loss, 'fermi_regularization': fermi_loss.item(), 'total': total_loss.item() } return total_loss, loss_dict def _compute_sensitivity_loss(self, sensitivity_matrices): """Encourage high sensitivity to rare events""" if not sensitivity_matrices.get('hit_attention'): return 0.0 # Maximize attention concentration (low entropy = high sensitivity) total_concentration = 0.0 count = 0 for sensitivity in sensitivity_matrices['hit_attention']: # sensitivity: [batch, heads, seq, seq] prob = F.softmax(sensitivity, dim=-1) # Entropy of attention distribution entropy = -torch.sum(prob * torch.log(prob + 1e-10), dim=-1) # Low entropy = high concentration = good concentration = -entropy.mean() total_concentration += concentration count += 1 if count > 0: return -total_concentration / count # Negative because we want to minimize entropy return 0.0 def _compute_fermi_regularization(self, outputs): """ Fermi-Dirac inspired regularization. Encourages weight distribution to follow quantum statistics. """ fused_rep = outputs['fused_representation'] # Compute energy levels (as proxy for quantum states) energy = torch.norm(fused_rep, dim=-1) # Magnitude as energy proxy # Fermi-Dirac distribution: f(E) = 1 / (exp((E - mu) / kT) + 1) # Here we just encourage exponential distribution of energies energy_sorted, _ = torch.sort(energy) energy_diff = energy_sorted[1:] - energy_sorted[:-1] # Encourage decreasing occupation fermi_loss = F.relu(energy_diff[:-1] - energy_diff[1:]).mean() return fermi_loss # ============================================================================= # SECTION 6: CCT-ODE TRAINING UTILITIES # ============================================================================= class AdaptiveSensitivityCallback: """ Callback for adaptive sensitivity adjustment during training. Monitors entropy evolution and adjusts training parameters accordingly. """ def __init__(self, model, sensitivity_target=0.1): self.model = model self.sensitivity_target = sensitivity_target self.entropy_history = [] self.phase = 'normal' # normal, sensitive, periodic def on_batch_end(self, batch_idx, loss_dict, outputs): """Called after each training batch""" # Extract entropy from outputs if outputs and 'sensitivity_matrices' in outputs: entropies = outputs['sensitivity_matrices'].get('layer_entropies', []) if entropies: current_entropy = entropies[-1] self.entropy_history.append(current_entropy) # Adaptive sensitivity adjustment if len(self.entropy_history) > 10: recent_avg = np.mean(self.entropy_history[-10:]) if recent_avg > self.sensitivity_target * 2: # Too much entropy - increase sensitivity self.phase = 'sensitive' return {'sensitivity_scale': 2.0, 'learning_rate_factor': 1.5} elif recent_avg < self.sensitivity_target * 0.5: # Very low entropy - stable mode self.phase = 'periodic' return {'sensitivity_scale': 0.5, 'learning_rate_factor': 0.5} else: self.phase = 'normal' return {'sensitivity_scale': 1.0, 'learning_rate_factor': 1.0} return {'sensitivity_scale': 1.0, 'learning_rate_factor': 1.0} def train_neutrino_cct( model: nn.Module, train_loader: DataLoader, val_loader: DataLoader, epochs: int = 50, device: str = 'cuda', lr: float = 1e-4, weight_decay: float = 0.01, rare_event_weight: float = 5.0 ): """ Train the neutrino detection model using CCT-ODE framework. Args: model: NeutrinoTransformer model train_loader: Training data loader val_loader: Validation data loader epochs: Number of training epochs device: 'cuda' or 'cpu' lr: Learning rate weight_decay: Weight decay for AdamW rare_event_weight: Weight for rare event detection loss Returns: history: dict of training metrics """ model = model.to(device) criterion = CCTODELoss(rare_event_weight=rare_event_weight) optimizer = optim.AdamW(model.parameters(), lr=lr, weight_decay=weight_decay) # Learning rate scheduler scheduler = optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=epochs) # Sensitivity callback sensitivity_callback = AdaptiveSensitivityCallback(model) history = { 'train_loss': [], 'val_loss': [], 'train_acc': [], 'val_acc': [], 'sensitivity': [], 'rare_event_auc': [] } best_val_loss = float('inf') for epoch in range(epochs): model.train() epoch_losses = [] epoch_correct = 0 epoch_total = 0 for batch_idx, batch in enumerate(train_loader): # Move data to device hits = batch['hits'].to(device) features = batch['features'].to(device) labels = batch['label'].to(device) # Forward pass optimizer.zero_grad() outputs = model(hits, features, track_sensitivity=True) # Compute loss loss, loss_dict = criterion(outputs, labels, outputs['sensitivity_matrices']) # Backward pass loss.backward() # Gradient clipping for stability torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0) optimizer.step() # Adaptive sensitivity adjustment adjust = sensitivity_callback.on_batch_end(batch_idx, loss_dict, outputs) # Metrics epoch_losses.append(loss.item()) preds = outputs['logits'].argmax(dim=1) epoch_correct += (preds == labels).sum().item() epoch_total += labels.size(0) # Reset sensitivity tracking model.reset_sensitivity_tracking() scheduler.step() # Epoch metrics train_loss = np.mean(epoch_losses) train_acc = epoch_correct / epoch_total # Validation model.eval() val_losses = [] val_correct = 0 val_total = 0 with torch.no_grad(): for batch in val_loader: hits = batch['hits'].to(device) features = batch['features'].to(device) labels = batch['label'].to(device) outputs = model(hits, features, track_sensitivity=False) loss, _ = criterion(outputs, labels, None) val_losses.append(loss.item()) preds = outputs['logits'].argmax(dim=1) val_correct += (preds == labels).sum().item() val_total += labels.size(0) val_loss = np.mean(val_losses) val_acc = val_correct / val_total # Store history history['train_loss'].append(train_loss) history['val_loss'].append(val_loss) history['train_acc'].append(train_acc) history['val_acc'].append(val_acc) history['sensitivity'].append(sensitivity_callback.phase) # Print progress print(f"Epoch {epoch+1}/{epochs}") print(f" Train Loss: {train_loss:.4f} | Train Acc: {train_acc:.4f}") print(f" Val Loss: {val_loss:.4f} | Val Acc: {val_acc:.4f}") print(f" Sensitivity Phase: {sensitivity_callback.phase}") print(f" LR: {scheduler.get_last_lr()[0]:.6f}") print() # Save best model if val_loss < best_val_loss: best_val_loss = val_loss torch.save(model.state_dict(), 'best_neutrino_cct_model.pt') print(f" ✓ Saved best model (val_loss: {val_loss:.4f})") return history # ============================================================================= # SECTION 7: EVALUATION AND ANALYSIS # ============================================================================= def evaluate_neutrino_detector(model, test_loader, device='cuda'): """ Comprehensive evaluation of neutrino detector. Computes: - Accuracy, Precision, Recall, F1 - ROC-AUC (important for rare event detection) - Confusion matrix - Sensitivity analysis per event type """ model.eval() model = model.to(device) all_preds = [] all_labels = [] all_rare_probs = [] all_hits = [] all_features = [] with torch.no_grad(): for batch in test_loader: hits = batch['hits'].to(device) features = batch['features'].to(device) labels = batch['label'].to(device) outputs = model(hits, features, track_sensitivity=False) preds = outputs['logits'].argmax(dim=1) rare_probs = outputs['rare_event_prob'].squeeze() all_preds.extend(preds.cpu().tolist()) all_labels.extend(labels.cpu().tolist()) all_rare_probs.extend(rare_probs.cpu().tolist()) all_hits.extend(hits.cpu()) all_features.extend(features.cpu()) # Convert to numpy all_preds = np.array(all_preds) all_labels = np.array(all_labels) all_rare_probs = np.array(all_rare_probs) # Compute metrics accuracy = (all_preds == all_labels).mean() # Confusion matrix tp = ((all_preds == 1) & (all_labels == 1)).sum() tn = ((all_preds == 0) & (all_labels == 0)).sum() fp = ((all_preds == 1) & (all_labels == 0)).sum() fn = ((all_preds == 0) & (all_labels == 1)).sum() precision = tp / (tp + fp + 1e-6) recall = tp / (tp + fn + 1e-6) f1 = 2 * precision * recall / (precision + recall + 1e-6) # ROC-AUC for rare event detection try: from sklearn.metrics import roc_auc_score, roc_curve auc = roc_auc_score(all_labels, all_rare_probs) fpr, tpr, thresholds = roc_curve(all_labels, all_rare_probs) except ImportError: auc = 0.5 fpr, tpr, thresholds = None, None, None print("=" * 60) print("NEUTRINO DETECTOR EVALUATION RESULTS") print("=" * 60) print(f"\nOverall Metrics:") print(f" Accuracy: {accuracy:.4f}") print(f" Precision: {precision:.4f}") print(f" Recall: {recall:.4f}") print(f" F1 Score: {f1:.4f}") print(f" ROC-AUC (Rare Event): {auc:.4f}") print(f"\nConfusion Matrix:") print(f" True Positives: {tp} | False Negatives: {fn}") print(f" False Positives: {fp} | True Negatives: {tn}") print(f"\nRare Event Detection:") print(f" Signal events: {all_labels.sum()}") print(f" Background events: {(1 - all_labels).sum()}") print(f" Detected signals: {tp}") print(f" Missed signals: {fn}") # Per-energy analysis print(f"\nPer-Energy Analysis:") for energy_bin in ['low', 'medium', 'high']: if energy_bin == 'low': mask = np.array([f[0] < -1.0 for f in all_features]) elif energy_bin == 'medium': mask = np.array([f[0] >= -1.0 and f[0] < 0.5 for f in all_features]) else: mask = np.array([f[0] >= 0.5 for f in all_features]) if mask.sum() > 0: bin_acc = (all_preds[mask] == all_labels[mask]).mean() print(f" {energy_bin.capitalize()} energy: Accuracy = {bin_acc:.4f} (n={mask.sum()})") print("=" * 60) return { 'accuracy': accuracy, 'precision': precision, 'recall': recall, 'f1': f1, 'auc': auc, 'confusion_matrix': {'tp': tp, 'tn': tn, 'fp': fp, 'fn': fn}, 'predictions': all_preds, 'labels': all_labels, 'rare_probs': all_rare_probs, 'fpr': fpr, 'tpr': tpr } def analyze_sensitivity_matrices(model, data_loader, device='cuda'): """ Analyze the sensitivity matrices learned by the model. This reveals which aspects of the neutrino detection the model is most sensitive to, validating the CCT-ODE framework. """ model.eval() model = model.to(device) all_sensitivities = [] with torch.no_grad(): for batch in data_loader: hits = batch['hits'].to(device) features = batch['features'].to(device) labels = batch['label'].to(device) outputs = model(hits, features, track_sensitivity=True) # Collect sensitivity matrices if outputs['sensitivity_matrices']: all_sensitivities.append(outputs['sensitivity_matrices']) model.reset_sensitivity_tracking() print("=" * 60) print("SENSITIVITY MATRIX ANALYSIS") print("=" * 60) print("\nAnalyzing layer-by-layer sensitivity...") print("Higher sensitivity = Model is more responsive to that aspect") # Aggregate analysis if all_sensitivities: hit_attn = [s['hit_attention'] for s in all_sensitivities if 'hit_attention' in s] if hit_attn: print(f"\nHit Attention Sensitivities:") print(f" Number of layers analyzed: {len(hit_attn)}") print(f" Average attention entropy: {np.mean([e.mean().item() for es in hit_attn for e in es]):.4f}") print("=" * 60) # ============================================================================= # SECTION 8: MAIN EXECUTION # ============================================================================= def main(): """ Main function to train and evaluate the neutrino detector. """ # Configuration CONFIG = { 'd_model': 256, 'n_heads': 8, 'n_layers': 6, 'd_ff': 1024, 'max_hits': 200, 'n_features': 12, 'n_classes': 2, 'dropout': 0.1, 'train_samples': 50000, 'val_samples': 10000, 'test_samples': 10000, 'signal_ratio': 0.1, # 10% signal, 90% background 'epochs': 30, 'batch_size': 64, 'lr': 1e-4, 'weight_decay': 0.01, 'rare_event_weight': 5.0, 'device': 'cuda' if torch.cuda.is_available() else 'cpu' } print("=" * 60) print("MMA-CCT: Multi-Messenger Astronomy - Conditional Collapse Theory") print("Neutrino Detection System") print("=" * 60) print("\nConfiguration:") for key, value in CONFIG.items(): print(f" {key}: {value}") print(f"\nDevice: {CONFIG['device']}") # Create datasets print("\n[1/6] Generating datasets...") train_dataset = NeutrinoDataset( num_samples=CONFIG['train_samples'], signal_ratio=CONFIG['signal_ratio'], device=CONFIG['device'] ) val_dataset = NeutrinoDataset( num_samples=CONFIG['val_samples'], signal_ratio=CONFIG['signal_ratio'], device=CONFIG['device'] ) test_dataset = NeutrinoDataset( num_samples=CONFIG['test_samples'], signal_ratio=CONFIG['signal_ratio'], device=CONFIG['device'] ) train_loader = DataLoader(train_dataset, batch_size=CONFIG['batch_size'], shuffle=True, num_workers=0) val_loader = DataLoader(val_dataset, batch_size=CONFIG['batch_size'], shuffle=False, num_workers=0) test_loader = DataLoader(test_dataset, batch_size=CONFIG['batch_size'], shuffle=False, num_workers=0) print(f" Train: {len(train_dataset)} samples") print(f" Val: {len(val_dataset)} samples") print(f" Test: {len(test_dataset)} samples") print(f" Signal ratio: {CONFIG['signal_ratio']*100:.1f}%") # Create model print("\n[2/6] Creating model...") model = NeutrinoTransformer( d_model=CONFIG['d_model'], n_heads=CONFIG['n_heads'], n_layers=CONFIG['n_layers'], d_ff=CONFIG['d_ff'], max_hits=CONFIG['max_hits'], n_features=CONFIG['n_features'], n_classes=CONFIG['n_classes'], dropout=CONFIG['dropout'] ) num_params = sum(p.numel() for p in model.parameters() if p.requires_grad) print(f" Model parameters: {num_params:,}") # Train print("\n[3/6] Training model...") print("-" * 60) history = train_neutrino_cct( model=model, train_loader=train_loader, val_loader=val_loader, epochs=CONFIG['epochs'], device=CONFIG['device'], lr=CONFIG['lr'], weight_decay=CONFIG['weight_decay'], rare_event_weight=CONFIG['rare_event_weight'] ) print("-" * 60) # Load best model print("\n[4/6] Loading best model...") model.load_state_dict(torch.load('best_neutrino_cct_model.pt')) print(" Best model loaded.") # Evaluate print("\n[5/6] Evaluating model...") results = evaluate_neutrino_detector(model, test_loader, device=CONFIG['device']) # Analyze sensitivities print("\n[6/6] Analyzing sensitivity matrices...") analyze_sensitivity_matrices(model, test_loader, device=CONFIG['device']) print("\n" + "=" * 60) print("TRAINING COMPLETE") print("=" * 60) print(f"Best validation loss: {min(history['val_loss']):.4f}") print(f"Final test accuracy: {results['accuracy']:.4f}") print(f"Final test AUC: {results['auc']:.4f}") print("\nModel saved to: best_neutrino_cct_model.pt") return model, results, history if __name__ == "__main__": model, results, history = main()