""" Ablation study framework for isolating the FCE contribution. Runs the engine under multiple configurations to answer the key scientific question: what predictions emerge that were not already built into the equations? Configurations: A Physics only (qec_enabled=False) B Physics + FCE (full pipeline) C Physics + random FCE (random corrections, same magnitude) D Physics + FCE (curvature signal disabled) E Physics + FCE (torsion weighting disabled) For each pair (B vs X), computes: - Paired t-test (two-sided) - Cohen's d effect size - 95% confidence interval on the difference References: [1] Cohen, "Statistical Power Analysis" (1988) [2] Efron & Tibshirani, "An Introduction to the Bootstrap" (1993) """ import time import logging import numpy as np import scipy.stats as stats from concurrent.futures import ProcessPoolExecutor, as_completed from dataclasses import dataclass, field from typing import List, Dict, Optional, Callable, Any from enum import Enum from .engine import FractalCorrectionEngine, EngineConfig, EvolutionResult from .fidelity import uhlmann_fidelity from .system_health import SystemHealthMonitor, get_default_monitor logger = logging.getLogger(__name__) class AblationMode(Enum): """FCE ablation configurations.""" PHYSICS_ONLY = "A_physics_only" FULL_FCE = "B_full_fce" RANDOM_FCE = "C_random_fce" NO_CURVATURE = "D_no_curvature" NO_TORSION = "E_no_torsion" @dataclass class AblationRunMetrics: """Metrics from a single ablation run.""" mode: str final_fidelity: float mean_fidelity: float final_fidelity_uncorrected: float error_reduction_factor: float n_corrections: int max_trace_error: float runtime_seconds: float # Conservation and stability metrics mean_trace_error: float = 0.0 max_positivity_violation: float = 0.0 mean_purity: float = 0.0 fidelity_variance: float = 0.0 # Var(F(t)) -- stability metric @dataclass class PairwiseComparison: """Statistical comparison between two ablation modes.""" mode_a: str mode_b: str metric_name: str mean_a: float mean_b: float mean_difference: float std_difference: float t_statistic: float p_value: float cohens_d: float ci_lower: float ci_upper: float significant: bool @dataclass class AblationReport: """Complete ablation study results.""" # Per-mode aggregate results mode_metrics: Dict[str, Dict[str, Any]] # Pairwise statistical comparisons (B vs each other) comparisons: List[PairwiseComparison] # Raw per-run data for reproducibility raw_runs: Dict[str, List[AblationRunMetrics]] n_trials: int def summary_table(self) -> str: """Format as publication-ready comparison table.""" lines = [] lines.append("Ablation Study Results") lines.append("=" * 90) lines.append( f"{'Configuration':<25} {'Final F (mean)':>14} " f"{'Final F (std)':>13} {'Err Reduction':>14} " f"{'Runtime (s)':>12}" ) lines.append("-" * 90) for mode_name in [m.value for m in AblationMode]: if mode_name not in self.mode_metrics: continue m = self.mode_metrics[mode_name] lines.append( f"{mode_name:<25} " f"{m['final_fidelity_mean']:>14.6f} " f"{m['final_fidelity_std']:>13.6f} " f"{m['error_reduction_mean']:>14.2f} " f"{m['runtime_mean']:>12.4f}" ) lines.append("-" * 90) lines.append(f"Trials per configuration: {self.n_trials}") lines.append("") # Conservation & Stability per Mode lines.append("Conservation & Stability per Mode") lines.append("=" * 90) lines.append( f"{'Configuration':<25} {'Trace Err':>12} " f"{'Positivity':>12} {'Purity':>12} {'F Variance':>12}" ) lines.append("-" * 90) for mode_name in [m.value for m in AblationMode]: if mode_name not in self.raw_runs: continue runs = self.raw_runs[mode_name] if not runs: continue mean_te = float(np.mean([r.mean_trace_error for r in runs])) mean_pv = float(np.mean([r.max_positivity_violation for r in runs])) mean_pur = float(np.mean([r.mean_purity for r in runs])) mean_fv = float(np.mean([r.fidelity_variance for r in runs])) lines.append( f"{mode_name:<25} " f"{mean_te:>12.2e} " f"{mean_pv:>12.2e} " f"{mean_pur:>12.6f} " f"{mean_fv:>12.2e}" ) lines.append("-" * 90) lines.append("") # Statistical tests with effect size interpretation lines.append("Pairwise Comparisons (B_full_fce vs others)") lines.append("=" * 90) lines.append( f"{'Comparison':<30} {'Cohen d':>8} {'Effect':>10} " f"{'t-stat':>8} {'p-value':>10} {'95% CI':>20} {'Sig?':>5}" ) lines.append("-" * 90) for c in self.comparisons: ci_str = f"[{c.ci_lower:+.4f}, {c.ci_upper:+.4f}]" # Effect size interpretation (Cohen 1988) abs_d = abs(c.cohens_d) if abs_d < 0.2: effect_label = "negligible" elif abs_d < 0.5: effect_label = "small" elif abs_d < 0.8: effect_label = "medium" else: effect_label = "large" lines.append( f"B vs {c.mode_b:<25} " f"{c.cohens_d:>8.3f} " f"{effect_label:>10} " f"{c.t_statistic:>8.3f} " f"{c.p_value:>10.2e} " f"{ci_str:>20} " f"{'*' if c.significant else '':>5}" ) lines.append("-" * 90) lines.append("* = significant at alpha=0.05") lines.append("Effect sizes: negligible (|d|<0.2), small (0.2-0.5), " "medium (0.5-0.8), large (>0.8)") return "\n".join(lines) def _ablation_worker(args: dict) -> dict: """ Top-level worker function for parallel ablation trials. Must be top-level (not a method) so it can be pickled by ProcessPoolExecutor across processes. """ mode_value = args["mode_value"] seed = args["seed"] H = args["H"] rho_init = args["rho_init"] system_params = args["system_params"] dt = args["dt"] n_steps = args["n_steps"] base_config_dict = args["base_config_dict"] mode = AblationMode(mode_value) np.random.seed(seed) base_config = EngineConfig(**base_config_dict) # Build engine config for this mode if mode == AblationMode.PHYSICS_ONLY: config = EngineConfig( feedback_gain=0.0, correction_threshold=1e10, qec_enabled=False, enable_thermal_monitoring=False, ) elif mode == AblationMode.FULL_FCE: config = EngineConfig( feedback_gain=base_config.feedback_gain, correction_threshold=base_config.correction_threshold, qec_enabled=True, enable_thermal_monitoring=False, ) elif mode == AblationMode.RANDOM_FCE: config = EngineConfig( feedback_gain=base_config.feedback_gain, correction_threshold=base_config.correction_threshold, qec_enabled=True, enable_thermal_monitoring=False, ) elif mode == AblationMode.NO_CURVATURE: config = EngineConfig( feedback_gain=base_config.feedback_gain, correction_threshold=0.0, qec_enabled=True, enable_thermal_monitoring=False, ) elif mode == AblationMode.NO_TORSION: config = EngineConfig( feedback_gain=base_config.feedback_gain, correction_threshold=base_config.correction_threshold, qec_enabled=True, enable_thermal_monitoring=False, ) else: config = base_config engine = FractalCorrectionEngine( hamiltonian=H, system_params=system_params, config=config, ) rho = rho_init.copy() if mode == AblationMode.RANDOM_FCE: dim = rho.shape[0] noise = np.random.randn(dim, dim) + 1j * np.random.randn(dim, dim) noise = (noise + noise.conj().T) / 2.0 noise *= 0.01 rho = rho + noise eigvals, eigvecs = np.linalg.eigh(rho) eigvals = np.maximum(eigvals, 0.0) rho = eigvecs @ np.diag(eigvals) @ eigvecs.conj().T rho /= np.trace(rho) t_start = time.perf_counter() result = engine.evolve(rho, dt, n_steps) t_end = time.perf_counter() error_corrected = 1.0 - result.fidelities[-1] error_uncorrected = 1.0 - result.fidelities_uncorrected[-1] err_reduction = ( error_uncorrected / error_corrected if error_corrected > 1e-15 else float('inf') ) max_trace_err = 0.0 mean_trace_err = 0.0 if result.trace_errors is not None: max_trace_err = float(np.max(result.trace_errors)) mean_trace_err = float(np.mean(result.trace_errors)) max_pos_viol = 0.0 if result.positivity_violations is not None: max_pos_viol = float(max(-np.min(result.positivity_violations), 0.0)) # Purity of final state rho_final = result.density_matrices[-1] final_purity = float(np.real(np.trace(rho_final @ rho_final))) # Fidelity variance (stability measure) fid_var = float(np.var(result.fidelities)) return { "mode": mode_value, "final_fidelity": float(result.fidelities[-1]), "mean_fidelity": float(np.mean(result.fidelities)), "final_fidelity_uncorrected": float(result.fidelities_uncorrected[-1]), "error_reduction_factor": err_reduction, "n_corrections": result.n_corrections, "max_trace_error": max_trace_err, "runtime_seconds": t_end - t_start, "mean_trace_error": mean_trace_err, "max_positivity_violation": max_pos_viol, "mean_purity": final_purity, "fidelity_variance": fid_var, } class AblationStudy: """ Runs ablation experiments to isolate the FCE contribution. Without ablation, reviewers cannot determine whether FCE genuinely improves accuracy or merely perturbs solutions. Supports parallel execution via ProcessPoolExecutor with dynamic worker reduction based on CPU temperature and memory pressure. """ def __init__( self, hamiltonian: np.ndarray, rho_init: np.ndarray, system_params: Optional[Dict[str, float]] = None, dt: float = 1e-6, n_steps: int = 50, base_config: Optional[EngineConfig] = None, ): """ Args: hamiltonian: System Hamiltonian rho_init: Initial density matrix system_params: Decoherence parameters dt: Time step n_steps: Steps per evolution base_config: Base FCE configuration (used for mode B) """ self.H = hamiltonian self.rho_init = rho_init self.system_params = system_params or { 'T1': 1e-3, 'T2': 5e-4, 'gate_fidelity': 0.999 } self.dt = dt self.n_steps = n_steps self.base_config = base_config or EngineConfig( feedback_gain=0.1, correction_threshold=0.001, qec_enabled=True ) def _make_engine(self, mode: AblationMode) -> FractalCorrectionEngine: """Create engine for a specific ablation mode.""" if mode == AblationMode.PHYSICS_ONLY: config = EngineConfig( feedback_gain=0.0, correction_threshold=1e10, # Never triggers qec_enabled=False, ) elif mode == AblationMode.FULL_FCE: config = EngineConfig( feedback_gain=self.base_config.feedback_gain, correction_threshold=self.base_config.correction_threshold, qec_enabled=True, ) elif mode == AblationMode.RANDOM_FCE: # Same config as full FCE; randomization applied post-hoc config = EngineConfig( feedback_gain=self.base_config.feedback_gain, correction_threshold=self.base_config.correction_threshold, qec_enabled=True, ) elif mode == AblationMode.NO_CURVATURE: # Disable curvature-based threshold (always correct) config = EngineConfig( feedback_gain=self.base_config.feedback_gain, correction_threshold=0.0, # Always triggers qec_enabled=True, ) elif mode == AblationMode.NO_TORSION: # Same config; torsion zeroed in post-processing config = EngineConfig( feedback_gain=self.base_config.feedback_gain, correction_threshold=self.base_config.correction_threshold, qec_enabled=True, ) else: config = self.base_config return FractalCorrectionEngine( hamiltonian=self.H, system_params=self.system_params, config=config, ) def _run_single( self, mode: AblationMode, seed: int ) -> AblationRunMetrics: """Run a single trial for a given mode.""" np.random.seed(seed) engine = self._make_engine(mode) # For random FCE mode: add random perturbation to initial state rho = self.rho_init.copy() if mode == AblationMode.RANDOM_FCE: dim = rho.shape[0] # Add small random Hermitian perturbation noise = np.random.randn(dim, dim) + 1j * np.random.randn(dim, dim) noise = (noise + noise.conj().T) / 2.0 noise *= 0.01 # Small perturbation rho = rho + noise # Re-normalize to valid density matrix eigvals, eigvecs = np.linalg.eigh(rho) eigvals = np.maximum(eigvals, 0.0) rho = eigvecs @ np.diag(eigvals) @ eigvecs.conj().T rho /= np.trace(rho) t_start = time.perf_counter() result = engine.evolve(rho, self.dt, self.n_steps) t_end = time.perf_counter() error_corrected = 1.0 - result.fidelities[-1] error_uncorrected = 1.0 - result.fidelities_uncorrected[-1] err_reduction = ( error_uncorrected / error_corrected if error_corrected > 1e-15 else float('inf') ) max_trace_err = 0.0 mean_trace_err = 0.0 if result.trace_errors is not None: max_trace_err = float(np.max(result.trace_errors)) mean_trace_err = float(np.mean(result.trace_errors)) max_pos_viol = 0.0 if result.positivity_violations is not None: max_pos_viol = float(max(-np.min(result.positivity_violations), 0.0)) rho_final = result.density_matrices[-1] final_purity = float(np.real(np.trace(rho_final @ rho_final))) fid_var = float(np.var(result.fidelities)) return AblationRunMetrics( mode=mode.value, final_fidelity=float(result.fidelities[-1]), mean_fidelity=float(np.mean(result.fidelities)), final_fidelity_uncorrected=float(result.fidelities_uncorrected[-1]), error_reduction_factor=err_reduction, n_corrections=result.n_corrections, max_trace_error=max_trace_err, runtime_seconds=t_end - t_start, mean_trace_error=mean_trace_err, max_positivity_violation=max_pos_viol, mean_purity=final_purity, fidelity_variance=fid_var, ) def _pairwise_test( self, runs_b: List[AblationRunMetrics], runs_x: List[AblationRunMetrics], mode_x_name: str, metric: str = 'final_fidelity', ) -> PairwiseComparison: """Paired t-test between mode B and mode X.""" vals_b = np.array([getattr(r, metric) for r in runs_b]) vals_x = np.array([getattr(r, metric) for r in runs_x]) # Handle inf values vals_b = np.where(np.isinf(vals_b), np.nan, vals_b) vals_x = np.where(np.isinf(vals_x), np.nan, vals_x) mask = ~(np.isnan(vals_b) | np.isnan(vals_x)) vals_b = vals_b[mask] vals_x = vals_x[mask] n = len(vals_b) if n < 2: return PairwiseComparison( mode_a=AblationMode.FULL_FCE.value, mode_b=mode_x_name, metric_name=metric, mean_a=float(np.mean(vals_b)) if len(vals_b) else 0.0, mean_b=float(np.mean(vals_x)) if len(vals_x) else 0.0, mean_difference=0.0, std_difference=0.0, t_statistic=0.0, p_value=1.0, cohens_d=0.0, ci_lower=0.0, ci_upper=0.0, significant=False, ) differences = vals_b - vals_x mean_diff = float(np.mean(differences)) std_diff = float(np.std(differences, ddof=1)) # Paired t-test t_stat, p_val = stats.ttest_rel(vals_b, vals_x) # Cohen's d for paired samples cohens_d = mean_diff / std_diff if std_diff > 1e-15 else 0.0 # 95% CI on mean difference se = std_diff / np.sqrt(n) t_crit = stats.t.ppf(0.975, df=n - 1) ci_lower = mean_diff - t_crit * se ci_upper = mean_diff + t_crit * se return PairwiseComparison( mode_a=AblationMode.FULL_FCE.value, mode_b=mode_x_name, metric_name=metric, mean_a=float(np.mean(vals_b)), mean_b=float(np.mean(vals_x)), mean_difference=mean_diff, std_difference=std_diff, t_statistic=float(t_stat), p_value=float(p_val), cohens_d=cohens_d, ci_lower=ci_lower, ci_upper=ci_upper, significant=p_val < 0.05, ) def run( self, n_trials: int = 100, modes: Optional[List[AblationMode]] = None, base_seed: int = 42, parallel: bool = True, max_workers: Optional[int] = None, ) -> AblationReport: """ Run the complete ablation study. Args: n_trials: Number of trials per configuration modes: Which modes to test (default: all 5) base_seed: RNG seed for reproducibility parallel: If True, run trials across processes max_workers: Max parallel workers (auto-detected if None) Returns: AblationReport with statistics and comparisons """ if modes is None: modes = list(AblationMode) # Determine worker count via SystemHealthMonitor monitor = get_default_monitor() if parallel: monitor.wait_until_cool() if max_workers is None: max_workers = monitor.get_recommended_workers() logger.info( "Ablation study: %d modes x %d trials, %d workers", len(modes), n_trials, max_workers, ) else: max_workers = 1 # Serializable config dict for worker processes base_config_dict = { "feedback_gain": self.base_config.feedback_gain, "correction_threshold": self.base_config.correction_threshold, "qec_enabled": self.base_config.qec_enabled, "prediction_steps": self.base_config.prediction_steps, "prediction_ds": self.base_config.prediction_ds, "health_window": self.base_config.health_window, "enable_thermal_monitoring": False, "health_check_interval": self.base_config.health_check_interval, } # Build all job descriptors jobs = [] for mode in modes: for trial in range(n_trials): seed = base_seed + trial * len(AblationMode) + list(AblationMode).index(mode) jobs.append({ "mode_value": mode.value, "seed": seed, "H": self.H, "rho_init": self.rho_init, "system_params": self.system_params, "dt": self.dt, "n_steps": self.n_steps, "base_config_dict": base_config_dict, }) # Execute jobs results_list = [] if parallel and max_workers > 1: with ProcessPoolExecutor(max_workers=max_workers) as executor: futures = { executor.submit(_ablation_worker, job): job for job in jobs } for future in as_completed(futures): try: result = future.result() results_list.append(result) except Exception as exc: job = futures[future] logger.error( "Ablation trial failed (mode=%s, seed=%s): %s", job["mode_value"], job["seed"], exc, ) # Cooldown between batch completion monitor.cooldown_between_steps() else: for job in jobs: try: result = _ablation_worker(job) results_list.append(result) except Exception as exc: logger.error( "Ablation trial failed (mode=%s): %s", job["mode_value"], exc, ) # Collect into per-mode lists raw_runs: Dict[str, List[AblationRunMetrics]] = { mode.value: [] for mode in modes } for res in results_list: metrics = AblationRunMetrics( mode=res["mode"], final_fidelity=res["final_fidelity"], mean_fidelity=res["mean_fidelity"], final_fidelity_uncorrected=res["final_fidelity_uncorrected"], error_reduction_factor=res["error_reduction_factor"], n_corrections=res["n_corrections"], max_trace_error=res["max_trace_error"], runtime_seconds=res["runtime_seconds"], mean_trace_error=res.get("mean_trace_error", 0.0), max_positivity_violation=res.get("max_positivity_violation", 0.0), mean_purity=res.get("mean_purity", 0.0), fidelity_variance=res.get("fidelity_variance", 0.0), ) raw_runs[res["mode"]].append(metrics) # Aggregate metrics per mode mode_metrics: Dict[str, Dict[str, Any]] = {} for mode_name, runs in raw_runs.items(): fids = [r.final_fidelity for r in runs] errs = [r.error_reduction_factor for r in runs] errs_finite = [e for e in errs if np.isfinite(e)] runtimes = [r.runtime_seconds for r in runs] corrections = [r.n_corrections for r in runs] mode_metrics[mode_name] = { 'final_fidelity_mean': float(np.mean(fids)), 'final_fidelity_std': float(np.std(fids)), 'final_fidelity_ci95': ( float(np.percentile(fids, 2.5)), float(np.percentile(fids, 97.5)), ), 'error_reduction_mean': float(np.mean(errs_finite)) if errs_finite else 0.0, 'error_reduction_std': float(np.std(errs_finite)) if errs_finite else 0.0, 'runtime_mean': float(np.mean(runtimes)), 'runtime_std': float(np.std(runtimes)), 'corrections_mean': float(np.mean(corrections)), } # Pairwise comparisons: B vs each other mode comparisons = [] if AblationMode.FULL_FCE.value in raw_runs: runs_b = raw_runs[AblationMode.FULL_FCE.value] for mode in modes: if mode == AblationMode.FULL_FCE: continue if mode.value in raw_runs: comp = self._pairwise_test( runs_b, raw_runs[mode.value], mode.value ) comparisons.append(comp) return AblationReport( mode_metrics=mode_metrics, comparisons=comparisons, raw_runs=raw_runs, n_trials=n_trials, )