""" Parameter sensitivity analysis for the FCE feedback gain (alpha). The FCE strength alpha appears throughout the framework. This module sweeps alpha across a range and measures output metrics to determine whether FCE genuinely improves accuracy or merely perturbs solutions. Sweep values: [0, 0.01, 0.02, 0.05, 0.10, 0.20, 0.50, 1.0] For each alpha, record: - Final fidelity - Error reduction factor - Number of corrections - Hausdorff dimension - Runtime Plot Error(alpha) to show the sensitivity landscape. References: [1] Saltelli et al., "Global Sensitivity Analysis" (2008) """ import time import logging import numpy as np from concurrent.futures import ProcessPoolExecutor, as_completed from dataclasses import dataclass from typing import List, Optional, Dict, Any from .engine import FractalCorrectionEngine, EngineConfig, EvolutionResult from .system_health import get_default_monitor logger = logging.getLogger(__name__) @dataclass class SensitivityPoint: """Metrics at a single alpha value.""" alpha: float final_fidelity: float mean_fidelity: float error_reduction: float n_corrections: int hausdorff_dimension: float runtime_seconds: float # Optional: std from repeated trials final_fidelity_std: float = 0.0 @dataclass class SensitivityResult: """Complete sensitivity sweep results.""" points: List[SensitivityPoint] optimal_alpha: float optimal_fidelity: float @property def alpha_values(self) -> np.ndarray: return np.array([p.alpha for p in self.points]) @property def fidelity_values(self) -> np.ndarray: return np.array([p.final_fidelity for p in self.points]) @property def error_values(self) -> np.ndarray: return 1.0 - self.fidelity_values def summary_table(self) -> str: """Format as publication-ready sensitivity table.""" lines = [] lines.append("Parameter Sensitivity: feedback_gain (alpha)") lines.append("=" * 85) lines.append( f"{'alpha':>8} {'Final F':>12} {'F (std)':>10} " f"{'Err Red':>10} {'Corrections':>12} " f"{'d_H':>8} {'Runtime':>10}" ) lines.append("-" * 85) for pt in self.points: marker = " <--" if abs(pt.alpha - self.optimal_alpha) < 1e-10 else "" lines.append( f"{pt.alpha:>8.4f} " f"{pt.final_fidelity:>12.6f} " f"{pt.final_fidelity_std:>10.6f} " f"{pt.error_reduction:>10.2f} " f"{pt.n_corrections:>12d} " f"{pt.hausdorff_dimension:>8.4f} " f"{pt.runtime_seconds:>10.4f}" f"{marker}" ) lines.append("-" * 85) lines.append(f"Optimal alpha: {self.optimal_alpha:.4f}") lines.append(f"Optimal fidelity: {self.optimal_fidelity:.8f}") return "\n".join(lines) @dataclass class ThresholdSensitivityPoint: """Metrics at a single correction_threshold value.""" threshold: float final_fidelity: float mean_fidelity: float error_reduction: float n_corrections: int hausdorff_dimension: float runtime_seconds: float final_fidelity_std: float = 0.0 @dataclass class ThresholdSensitivityResult: """Complete threshold sweep results.""" points: List[ThresholdSensitivityPoint] optimal_threshold: float optimal_fidelity: float def summary_table(self) -> str: """Format as publication-ready threshold sensitivity table.""" lines = [] lines.append("Parameter Sensitivity: correction_threshold") lines.append("=" * 85) lines.append( f"{'threshold':>10} {'Final F':>12} {'F (std)':>10} " f"{'Err Red':>10} {'Corrections':>12} " f"{'d_H':>8} {'Runtime':>10}" ) lines.append("-" * 85) for pt in self.points: marker = " <--" if abs(pt.threshold - self.optimal_threshold) < 1e-15 else "" lines.append( f"{pt.threshold:>10.5f} " f"{pt.final_fidelity:>12.6f} " f"{pt.final_fidelity_std:>10.6f} " f"{pt.error_reduction:>10.2f} " f"{pt.n_corrections:>12d} " f"{pt.hausdorff_dimension:>8.4f} " f"{pt.runtime_seconds:>10.4f}" f"{marker}" ) lines.append("-" * 85) lines.append(f"Optimal threshold: {self.optimal_threshold:.5f}") lines.append(f"Optimal fidelity: {self.optimal_fidelity:.8f}") return "\n".join(lines) @dataclass class Sweep2DPoint: """Metrics at a single (alpha, threshold) pair.""" alpha: float threshold: float final_fidelity: float n_corrections: int @dataclass class Sweep2DResult: """Complete 2D sweep (alpha x threshold) results.""" points: List[Sweep2DPoint] alpha_values: List[float] threshold_values: List[float] fidelity_grid: np.ndarray # shape (n_alpha, n_threshold) optimal_alpha: float optimal_threshold: float optimal_fidelity: float def summary_table(self) -> str: """Format as publication-ready 2D sweep table.""" lines = [] lines.append("2D Parameter Sweep: alpha x correction_threshold") lines.append("=" * 80) # Header row with threshold values header = f"{'alpha \\\\ thresh':>14}" for t in self.threshold_values: header += f" {t:>10.4f}" lines.append(header) lines.append("-" * 80) # Grid rows for i, alpha in enumerate(self.alpha_values): row = f"{alpha:>14.4f}" for j in range(len(self.threshold_values)): fid = self.fidelity_grid[i, j] row += f" {fid:>10.6f}" lines.append(row) lines.append("-" * 80) lines.append( f"Optimal: alpha={self.optimal_alpha:.4f}, " f"threshold={self.optimal_threshold:.5f}, " f"fidelity={self.optimal_fidelity:.8f}" ) return "\n".join(lines) def _sensitivity_worker(args: dict) -> dict: """ Top-level worker for parallel sensitivity sweep runs. Must be top-level for pickling by ProcessPoolExecutor. """ alpha = args["alpha"] seed = args["seed"] rep = args["rep"] H = args["H"] rho_init = args["rho_init"] system_params = args["system_params"] dt = args["dt"] n_steps = args["n_steps"] np.random.seed(seed + rep) config = EngineConfig( feedback_gain=alpha, correction_threshold=0.001, qec_enabled=(alpha > 0), enable_thermal_monitoring=False, ) engine = FractalCorrectionEngine( hamiltonian=H, system_params=system_params, config=config, ) t_start = time.perf_counter() result = engine.evolve(rho_init, dt, n_steps) t_end = time.perf_counter() fid = float(result.fidelities[-1]) err_corr = 1.0 - fid err_uncorr = 1.0 - float(result.fidelities_uncorrected[-1]) err_red = ( err_uncorr / err_corr if err_corr > 1e-15 else float('inf') ) if not np.isfinite(err_red): err_red = 0.0 return { "alpha": alpha, "rep": rep, "final_fidelity": fid, "error_reduction": err_red, "n_corrections": result.n_corrections, "hausdorff_dimension": result.hausdorff_dimension, "runtime_seconds": t_end - t_start, } def _threshold_worker(args: dict) -> dict: """ Top-level worker for parallel threshold sweep runs. Must be top-level for pickling by ProcessPoolExecutor. """ threshold = args["threshold"] alpha = args["alpha"] seed = args["seed"] rep = args["rep"] H = args["H"] rho_init = args["rho_init"] system_params = args["system_params"] dt = args["dt"] n_steps = args["n_steps"] np.random.seed(seed + rep) config = EngineConfig( feedback_gain=alpha, correction_threshold=threshold, qec_enabled=(alpha > 0), enable_thermal_monitoring=False, ) engine = FractalCorrectionEngine( hamiltonian=H, system_params=system_params, config=config, ) t_start = time.perf_counter() result = engine.evolve(rho_init, dt, n_steps) t_end = time.perf_counter() fid = float(result.fidelities[-1]) err_corr = 1.0 - fid err_uncorr = 1.0 - float(result.fidelities_uncorrected[-1]) err_red = ( err_uncorr / err_corr if err_corr > 1e-15 else float('inf') ) if not np.isfinite(err_red): err_red = 0.0 return { "threshold": threshold, "alpha": alpha, "rep": rep, "final_fidelity": fid, "error_reduction": err_red, "n_corrections": result.n_corrections, "hausdorff_dimension": result.hausdorff_dimension, "runtime_seconds": t_end - t_start, } class SensitivitySweep: """ Sweeps the FCE feedback gain (alpha) to map parameter sensitivity. Shows whether FCE genuinely improves accuracy as a function of correction strength, or whether results are insensitive to alpha. Supports parallel execution of independent alpha/repeat runs via ProcessPoolExecutor with dynamic worker reduction. """ 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, ): """ Args: hamiltonian: System Hamiltonian rho_init: Initial density matrix system_params: Decoherence parameters dt: Time step n_steps: Steps per evolution """ 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 def run( self, alpha_values: Optional[List[float]] = None, n_repeats: int = 1, seed: int = 42, parallel: bool = True, max_workers: Optional[int] = None, ) -> SensitivityResult: """ Run sensitivity sweep. Args: alpha_values: List of feedback_gain values to test. Default: [0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0] n_repeats: Number of repeats per alpha (for error bars) seed: RNG seed parallel: If True, run alpha values in parallel max_workers: Max parallel workers (auto-detected if None) Returns: SensitivityResult with metrics at each alpha """ if alpha_values is None: alpha_values = [0.0, 0.01, 0.02, 0.05, 0.10, 0.20, 0.50, 1.0] # Determine worker count monitor = get_default_monitor() if parallel: monitor.wait_until_cool() if max_workers is None: max_workers = monitor.get_recommended_workers() logger.info( "Sensitivity sweep: %d alphas x %d repeats, %d workers", len(alpha_values), n_repeats, max_workers, ) else: max_workers = 1 # Build all jobs jobs = [] for alpha in alpha_values: for rep in range(n_repeats): jobs.append({ "alpha": alpha, "seed": seed, "rep": rep, "H": self.H, "rho_init": self.rho_init, "system_params": self.system_params, "dt": self.dt, "n_steps": self.n_steps, }) # Execute jobs results_list = [] if parallel and max_workers > 1: with ProcessPoolExecutor(max_workers=max_workers) as executor: futures = { executor.submit(_sensitivity_worker, job): job for job in jobs } for future in as_completed(futures): try: results_list.append(future.result()) except Exception as exc: job = futures[future] logger.error( "Sensitivity run failed (alpha=%.4f, rep=%d): %s", job["alpha"], job["rep"], exc, ) monitor.cooldown_between_steps() else: for job in jobs: try: results_list.append(_sensitivity_worker(job)) except Exception as exc: logger.error( "Sensitivity run failed (alpha=%.4f): %s", job["alpha"], exc, ) # Group results by alpha alpha_results: Dict[float, List[dict]] = {a: [] for a in alpha_values} for res in results_list: alpha_results[res["alpha"]].append(res) points = [] for alpha in alpha_values: runs = alpha_results[alpha] if not runs: continue fids = [r["final_fidelity"] for r in runs] err_reds = [r["error_reduction"] for r in runs] corrections = [r["n_corrections"] for r in runs] d_Hs = [r["hausdorff_dimension"] for r in runs] runtimes = [r["runtime_seconds"] for r in runs] points.append(SensitivityPoint( alpha=alpha, final_fidelity=float(np.mean(fids)), mean_fidelity=float(np.mean(fids)), error_reduction=float(np.mean(err_reds)), n_corrections=int(np.mean(corrections)), hausdorff_dimension=float(np.mean(d_Hs)), runtime_seconds=float(np.mean(runtimes)), final_fidelity_std=float(np.std(fids)) if n_repeats > 1 else 0.0, )) # Find optimal alpha fidelities = [p.final_fidelity for p in points] best_idx = np.argmax(fidelities) return SensitivityResult( points=points, optimal_alpha=points[best_idx].alpha, optimal_fidelity=points[best_idx].final_fidelity, ) def run_threshold_sweep( self, threshold_values: Optional[List[float]] = None, alpha: float = 0.1, n_repeats: int = 1, seed: int = 42, parallel: bool = True, max_workers: Optional[int] = None, ) -> ThresholdSensitivityResult: """ Sweep correction_threshold at fixed alpha. Args: threshold_values: Thresholds to test. Default: [0.0, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1] alpha: Fixed feedback gain n_repeats: Repeats per threshold seed: RNG seed parallel: If True, run in parallel max_workers: Max parallel workers Returns: ThresholdSensitivityResult with metrics at each threshold """ if threshold_values is None: threshold_values = [0.0, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1] monitor = get_default_monitor() if parallel: monitor.wait_until_cool() if max_workers is None: max_workers = monitor.get_recommended_workers() else: max_workers = 1 jobs = [] for threshold in threshold_values: for rep in range(n_repeats): jobs.append({ "threshold": threshold, "alpha": alpha, "seed": seed, "rep": rep, "H": self.H, "rho_init": self.rho_init, "system_params": self.system_params, "dt": self.dt, "n_steps": self.n_steps, }) results_list = [] if parallel and max_workers > 1: with ProcessPoolExecutor(max_workers=max_workers) as executor: futures = { executor.submit(_threshold_worker, job): job for job in jobs } for future in as_completed(futures): try: results_list.append(future.result()) except Exception as exc: job = futures[future] logger.error( "Threshold sweep failed (t=%.5f, rep=%d): %s", job["threshold"], job["rep"], exc, ) monitor.cooldown_between_steps() else: for job in jobs: try: results_list.append(_threshold_worker(job)) except Exception as exc: logger.error( "Threshold sweep failed (t=%.5f): %s", job["threshold"], exc, ) # Group by threshold thresh_results: Dict[float, List[dict]] = {t: [] for t in threshold_values} for res in results_list: thresh_results[res["threshold"]].append(res) points = [] for threshold in threshold_values: runs = thresh_results[threshold] if not runs: continue fids = [r["final_fidelity"] for r in runs] err_reds = [r["error_reduction"] for r in runs] corrections = [r["n_corrections"] for r in runs] d_Hs = [r["hausdorff_dimension"] for r in runs] runtimes = [r["runtime_seconds"] for r in runs] points.append(ThresholdSensitivityPoint( threshold=threshold, final_fidelity=float(np.mean(fids)), mean_fidelity=float(np.mean(fids)), error_reduction=float(np.mean(err_reds)), n_corrections=int(np.mean(corrections)), hausdorff_dimension=float(np.mean(d_Hs)), runtime_seconds=float(np.mean(runtimes)), final_fidelity_std=float(np.std(fids)) if n_repeats > 1 else 0.0, )) fidelities = [p.final_fidelity for p in points] best_idx = np.argmax(fidelities) return ThresholdSensitivityResult( points=points, optimal_threshold=points[best_idx].threshold, optimal_fidelity=points[best_idx].final_fidelity, ) def run_2d_sweep( self, alpha_values: Optional[List[float]] = None, threshold_values: Optional[List[float]] = None, seed: int = 42, parallel: bool = True, max_workers: Optional[int] = None, ) -> Sweep2DResult: """ 2D parameter sweep over alpha x correction_threshold. Args: alpha_values: Feedback gain values. Default: [0.0, 0.05, 0.1, 0.2, 0.5] threshold_values: Threshold values. Default: [0.0, 0.001, 0.01, 0.1] seed: RNG seed parallel: If True, run in parallel max_workers: Max parallel workers Returns: Sweep2DResult with fidelity grid and optimal point """ if alpha_values is None: alpha_values = [0.0, 0.05, 0.1, 0.2, 0.5] if threshold_values is None: threshold_values = [0.0, 0.001, 0.01, 0.1] monitor = get_default_monitor() if parallel: monitor.wait_until_cool() if max_workers is None: max_workers = monitor.get_recommended_workers() else: max_workers = 1 jobs = [] for alpha in alpha_values: for threshold in threshold_values: jobs.append({ "threshold": threshold, "alpha": alpha, "seed": seed, "rep": 0, "H": self.H, "rho_init": self.rho_init, "system_params": self.system_params, "dt": self.dt, "n_steps": self.n_steps, }) results_list = [] if parallel and max_workers > 1: with ProcessPoolExecutor(max_workers=max_workers) as executor: futures = { executor.submit(_threshold_worker, job): job for job in jobs } for future in as_completed(futures): try: results_list.append(future.result()) except Exception as exc: job = futures[future] logger.error( "2D sweep failed (a=%.4f, t=%.5f): %s", job["alpha"], job["threshold"], exc, ) monitor.cooldown_between_steps() else: for job in jobs: try: results_list.append(_threshold_worker(job)) except Exception as exc: logger.error( "2D sweep failed (a=%.4f, t=%.5f): %s", job["alpha"], job["threshold"], exc, ) # Build grid fidelity_grid = np.zeros((len(alpha_values), len(threshold_values))) points = [] for res in results_list: a_idx = alpha_values.index(res["alpha"]) t_idx = threshold_values.index(res["threshold"]) fidelity_grid[a_idx, t_idx] = res["final_fidelity"] points.append(Sweep2DPoint( alpha=res["alpha"], threshold=res["threshold"], final_fidelity=res["final_fidelity"], n_corrections=res["n_corrections"], )) # Find optimal best_flat = np.argmax(fidelity_grid) best_a_idx, best_t_idx = np.unravel_index(best_flat, fidelity_grid.shape) return Sweep2DResult( points=points, alpha_values=alpha_values, threshold_values=threshold_values, fidelity_grid=fidelity_grid, optimal_alpha=alpha_values[best_a_idx], optimal_threshold=threshold_values[best_t_idx], optimal_fidelity=float(fidelity_grid[best_a_idx, best_t_idx]), )