"""
Run the GBS → GFFK ODE-CCT Pipeline
"""

# --- Setup: 4-mode GBS circuit ---
n = 4
squeezing = np.ones(n) * 0.8  # 0.8 squeezing per mode

# Random interferometer (unitary matrix)
U = np.random.randn(n, n) + 1j * np.random.randn(n, n)
Q, R = np.linalg.qr(U)
d = np.diagonal(R)
U = Q @ np.diag(d / np.abs(d))

# Create the GBS-GFFK system
gbs = GBSGFFK(n_modes=n, squeezing_params=squeezing, unitary=U)

print("=" * 60)
print("GBS → GFFK ODE-CCT PIPELINE")
print("=" * 60)
print(f"Modes: {n}")
print(f"Squeezing: {squeezing}")
print(f"Interferometer: Random Haar unitary")
print(f"Hamiltonian eigenvalues: {np.linalg.eigvalsh(gbs.H)}")
print("=" * 60)

# --- Evolve the state ---
initial_state = np.ones(n, dtype=complex) * 0.5  # Coherent state input
time = 2.0
n_steps = 1000

trajectory, profiles = gbs.evolve_state(initial_state, time, n_steps)

print("\n--- Trajectory Profiles ---")
for mode in range(n):
    p = profiles[mode]
    classification = gbs.classify_attractor(trajectory, mode)
    print(f"Mode {mode}: a={p['a']:.3f}, b={p['b']:.3f} → {classification}")

# --- Compute output distribution ---
distribution = gbs.compute_output_distribution(max_photons=2)

print("\n--- Output Distribution (Top 10 patterns) ---")
sorted_dist = sorted(distribution.items(), key=lambda x: x[1], reverse=True)
for pattern, prob in sorted_dist[:10]:
    print(f"  Pattern {pattern}: P = {prob:.6f}")

# --- Quantum advantage indicator ---
advantage = gbs.get_quantum_advantage_indicator()
print("\n--- Quantum Advantage Indicator ---")
print(f"  Score: {advantage['advantage_score']:.3f}")
print(f"  Distinguishable patterns: {advantage['n_patterns']}")
print(f"  Distribution entropy: {advantage['entropy']:.3f} bits")
print(f"  Total possible patterns: {advantage['total_patterns']}")
print("=" * 60)