""" GBS → GFFK ODE-CCT Pipeline (FIXED v4.1) ======================================== Fixed: Variable name conflict (time variable vs time module) """ import numpy as np from math import factorial import time # ← This is the time MODULE import sys def hafnian(M): n = M.shape[0] if n % 2 != 0: return 0.0 if n == 0: return 1.0 if n == 2: return M[0, 1] return _hafnian_direct(M) def _hafnian_direct(M): n = M.shape[0] if n == 0: return 1.0 if n == 2: return M[0, 1] haf = 0.0 for j in range(1, n): indices = list(range(1, n))[:j] + list(range(j+1, n)) sub_M = M[np.ix_(indices, indices)] haf += M[0, j] * _hafnian_direct(sub_M) return haf def progress_bar(current, total, prefix="", suffix="", length=50): percent = (current / total) * 100 filled = int(length * current // total) bar = '█' * filled + '-' * (length - filled) sys.stdout.write(f'\r{prefix} |{bar}| {percent:.1f}% {suffix}') sys.stdout.flush() if current == total: print() class GBSGFFK_Fixed: def __init__(self, n_modes, squeezing_params=None, unitary=None): self.n = n_modes if squeezing_params is None: self.squeezing = np.ones(n_modes) * 0.8 else: self.squeezing = np.array(squeezing_params) if unitary is None: self.U = self._random_haar_unitary(self.n) else: self.U = np.array(unitary, dtype=complex) D_squeeze = np.diag(2 * self.squeezing + 1) self.Sigma = self.U @ D_squeeze @ self.U.conj().T D_H = np.diag(self.squeezing) self.H = self.U @ D_H @ self.U.conj().T self.F = -1j * self.H self.log_det_sigma = self._safe_log_det() def _safe_log_det(self): try: eigenvalues = np.linalg.eigvalsh(self.Sigma) eigenvalues = np.maximum(eigenvalues, 1e-10) return np.sum(np.log(eigenvalues)) except: return 0.0 def _random_haar_unitary(self, n): Z = np.random.randn(n, n) + 1j * np.random.randn(n, n) Q, R = np.linalg.qr(Z) d = np.diagonal(R) return Q @ np.diag(d / np.abs(d)) def evolve_state(self, initial_state, t_end, n_steps=1000): # ← RENAMED: 'time' → 't_end' to avoid conflict with time module print("\n[1/3] Evolving state (ODE integration)...", end=" ", flush=True) t_start = time.time() # ← Now works correctly trajectory = np.zeros((n_steps + 1, self.n), dtype=complex) trajectory[0] = initial_state.copy() dt = t_end / n_steps # ← Use t_end for step in range(n_steps): z = trajectory[step] k1 = self.F @ z k2 = self.F @ (z + 0.5 * dt * k1) k3 = self.F @ (z + 0.5 * dt * k2) k4 = self.F @ (z + dt * k3) trajectory[step + 1] = z + (dt / 6) * (k1 + 2*k2 + 2*k3 + k4) if step % 100 == 0: progress_bar(step, n_steps, " Step", f" / {n_steps}") progress_bar(n_steps, n_steps, " Step", f" / {n_steps}", length=30) elapsed = time.time() - t_start # ← Now works correctly print(f" Done! ({elapsed:.2f}s)") profiles = {} for mode in range(self.n): m = np.abs(trajectory[:, mode]) a = np.mean(m) b = np.std(m) / (np.mean(m) + 1e-10) profiles[mode] = {'a': a, 'b': b} return trajectory, profiles def compute_output_distribution(self, max_photons=3): print("\n[2/3] Computing output distribution...") A = self.Sigma - np.eye(self.n) distribution = {} patterns = self._enumerate_patterns(max_photons) total_patterns = len(patterns) for i, pattern in enumerate(patterns): progress_bar(i+1, total_patterns, f" Pattern", f" / {total_patterns}", length=30) prob = self._compute_pattern_probability(pattern, A) if prob is not None and prob > 1e-15: distribution[pattern] = prob progress_bar(total_patterns, total_patterns, f" Pattern", f" / {total_patterns}", length=30) total = sum(distribution.values()) if total > 0: for k in distribution: distribution[k] /= total print(f"\n Found {len(distribution)} valid patterns out of {total_patterns}") return distribution def _compute_pattern_probability(self, pattern, A): total_photons = sum(pattern) if total_photons == 0: try: prob = np.exp(-0.5 * self.log_det_sigma) return prob except: return 0.0 A_sub = self._build_sub_matrix(pattern, A) haf = hafnian(A_sub) try: det_sigma = np.linalg.det(self.Sigma) if abs(det_sigma) < 1e-20: norm = np.exp(0.5 * self.log_det_sigma) else: norm = np.sqrt(abs(det_sigma)) except: norm = 1.0 if norm == 0: return None prob = np.abs(haf) ** 2 / ( np.prod([factorial(s) for s in pattern]) * norm ) return prob def _build_sub_matrix(self, pattern, A): n = len(pattern) total = sum(pattern) if total == 0: return np.array([[]]) A_sub = np.zeros((total, total), dtype=complex) row_offset = 0 for j in range(n): n_j = pattern[j] if n_j > 0: for k in range(n_j): row = row_offset + k A_sub[row, row] = A[j, j] for m in range(row + 1, total): A_sub[row, m] = A[j, j] A_sub[m, row] = A[j, j] row_offset += n_j A_sub = (A_sub + A_sub.T) / 2 return A_sub def _enumerate_patterns(self, max_photons): patterns = [] def recurse(mode, current): if mode == self.n: patterns.append(tuple(current)) return for k in range(max_photons + 1): current.append(k) recurse(mode + 1, current) current.pop() recurse(0, []) return patterns def get_quantum_advantage_indicator(self, distribution): probs = np.array([p for p in distribution.values() if p > 1e-10]) if len(probs) == 0: return {'advantage_score': 0, 'n_patterns': 0, 'entropy': 0} entropy = 0.0 for p in probs: if p > 0: entropy -= p * np.log2(p) n_patterns = len(probs) advantage = entropy * np.log2(n_patterns + 1) return { 'advantage_score': advantage, 'n_patterns': n_patterns, 'entropy': entropy, 'total_patterns': len(distribution), } if __name__ == "__main__": n = 4 squeezing = np.ones(n) * 0.8 print("=" * 60) print("GBS → GFFK ODE-CCT PIPELINE (v4.1 — Fixed)") print("=" * 60) gbs = GBSGFFK_Fixed(n_modes=n, squeezing_params=squeezing) print(f"\n[Setup] Modes: {n}, Squeezing: {squeezing}") print(f" Covariance det: {np.linalg.det(gbs.Sigma):.6f}") print(f" Log-det: {gbs.log_det_sigma:.6f}") print("=" * 60) initial_state = np.ones(n, dtype=complex) * 0.5 t_end = 2.0 # ← Use t_end instead of time # PASS t_end, not time (which is now the module) trajectory, profiles = gbs.evolve_state(initial_state, t_end, n_steps=500) print("\n--- Trajectory Profiles ---") for mode in range(n): p = profiles[mode] print(f" Mode {mode}: a={p['a']:.4f}, b={p['b']:.4f}") distribution = gbs.compute_output_distribution(max_photons=3) print(f"\n[3/3] Results:") print(f" Patterns: {len(distribution)}") sorted_dist = sorted(distribution.items(), key=lambda x: x[1], reverse=True) print(f" Top 10 patterns:") for pattern, prob in sorted_dist[:10]: print(f" {pattern}: P = {prob:.8f}") total_prob = sum(distribution.values()) print(f"\n Total probability: {total_prob:.8f}") advantage = gbs.get_quantum_advantage_indicator(distribution) print(f"\n Quantum advantage:") print(f" Score: {advantage['advantage_score']:.4f}") print(f" Entropy: {advantage['entropy']:.4f} bits") print("=" * 60)