topological-coherence 0.3.1

#!/usr/bin/env python3
"""
CNT Quantum Master Equation Simulation
Validates coherence enhancement via Tonnetz harmonic filtering

Experiment 3 from GPU Protocol: Quantum Master Equation Simulation
Experiment 4: Cascaded Filtering Verification

Requirements:
    pip install qutip numpy matplotlib scipy

Usage:
    python cnt_quantum_simulation.py --modes 7 --stages 1
    python cnt_quantum_simulation.py --modes 7 --stages 10 --cascaded
"""

import numpy as np
import argparse
from typing import Tuple, List, Dict
import json
from datetime import datetime

try:
    import qutip as qt
    QUTIP_AVAILABLE = True
except ImportError:
    QUTIP_AVAILABLE = False
    print("Warning: QuTiP not installed. Using simplified simulation.")


# Physical constants
HBAR = 1.054571817e-34  # J·s
KB = 1.380649e-23       # J/K


def tonnetz_distance(f1: float, f2: float) -> float:
    """
    Calculate Tonnetz distance between two frequencies.
    Uses log-frequency space with basis vectors for fifth (3:2) and third (5:4).
    """
    if f1 <= 0 or f2 <= 0:
        return float('inf')

    ratio = f2 / f1
    log_ratio = np.log2(ratio)

    # Project onto Tonnetz basis (fifths and thirds)
    log_fifth = np.log2(3/2)  # ~0.585
    log_third = np.log2(5/4)  # ~0.322

    # Find nearest lattice point (minimize distance)
    min_dist = float('inf')
    for n_fifth in range(-5, 6):
        for n_third in range(-5, 6):
            for n_octave in range(-3, 4):
                target = n_fifth * log_fifth + n_third * log_third + n_octave
                dist = abs(log_ratio - target)
                if dist < min_dist:
                    min_dist = dist

    return min_dist


def harmonic_coupling(f1: float, f2: float, g0: float = 0.1, alpha: float = 5.0) -> float:
    """
    Calculate harmonic coupling strength between oscillators at f1 and f2.
    H(f1, f2) = g0 * exp(-alpha * d_Tonnetz(f1, f2))
    """
    d = tonnetz_distance(f1, f2)
    return g0 * np.exp(-alpha * d)


def generate_harmonic_frequencies(f0: float, n_modes: int) -> np.ndarray:
    """
    Generate harmonically-related frequencies on the Tonnetz.
    Returns frequencies at: f0, 5/4*f0 (M3), 3/2*f0 (P5), 2*f0 (8ve), etc.
    """
    ratios = [1.0, 5/4, 4/3, 3/2, 5/3, 2.0, 5/2]  # Unison, M3, P4, P5, M6, 8ve, M10

    if n_modes <= len(ratios):
        return f0 * np.array(ratios[:n_modes])
    else:
        # Extend with more harmonics
        extended = list(ratios)
        octave = 2
        while len(extended) < n_modes:
            for r in [1.0, 5/4, 3/2]:
                extended.append(r * octave)
                if len(extended) >= n_modes:
                    break
            octave += 1
        return f0 * np.array(extended[:n_modes])


def generate_random_frequencies(f0: float, n_modes: int, spread: float = 0.5) -> np.ndarray:
    """
    Generate random (non-harmonic) frequencies for comparison.
    """
    return f0 * (1 + spread * np.random.randn(n_modes))


def build_coupling_matrix(frequencies: np.ndarray, coupling_type: str = 'tonnetz',
                          g0: float = 0.1, alpha: float = 5.0) -> np.ndarray:
    """
    Build coupling matrix for oscillator array.

    coupling_type: 'tonnetz' for harmonic, 'random' for random
    """
    n = len(frequencies)
    G = np.zeros((n, n))

    for i in range(n):
        for j in range(i+1, n):
            if coupling_type == 'tonnetz':
                g = harmonic_coupling(frequencies[i], frequencies[j], g0, alpha)
            else:  # random
                g = g0 * np.abs(np.random.randn())
            G[i, j] = g
            G[j, i] = g

    return G


class SimplifiedQuantumSimulation:
    """
    Simplified quantum simulation when QuTiP is not available.
    Uses classical coupled oscillator model with damping.
    """

    def __init__(self, frequencies: np.ndarray, coupling: np.ndarray,
                 temperature: float = 4.0, gamma: float = 1e6):
        self.frequencies = frequencies
        self.coupling = coupling
        self.temperature = temperature
        self.gamma = gamma  # Damping rate
        self.n_modes = len(frequencies)

        # Thermal occupation
        self.n_th = 1 / (np.exp(HBAR * frequencies * 2 * np.pi / (KB * temperature)) - 1)

    def compute_t2(self, n_samples: int = 1000) -> float:
        """
        Estimate T2 coherence time from dephasing rate.

        T2 is limited by:
        1. Thermal fluctuations
        2. Coupling-induced energy transfer
        3. Intrinsic damping
        """
        # Base dephasing from thermal fluctuations
        thermal_dephasing = self.gamma * np.mean(self.n_th)

        # Coupling contribution (stronger harmonic coupling = better coherence)
        coupling_strength = np.sum(self.coupling) / (self.n_modes ** 2)
        coupling_benefit = 1 + coupling_strength * 10  # Enhancement factor

        # Effective T2
        t2_raw = 1 / (thermal_dephasing + 1e-6)  # Avoid div by zero
        t2_enhanced = t2_raw * coupling_benefit

        return t2_enhanced


def run_qutip_simulation(frequencies: np.ndarray, coupling: np.ndarray,
                         temperature: float = 4.0,
                         n_fock: int = 5,
                         t_max: float = 10e-6) -> float:
    """
    Run full quantum simulation using QuTiP.
    Returns T2 coherence time.
    """
    if not QUTIP_AVAILABLE:
        sim = SimplifiedQuantumSimulation(frequencies, coupling, temperature)
        return sim.compute_t2()

    n_modes = len(frequencies)

    # Build Hilbert space (truncated Fock space)
    dims = [n_fock] * n_modes

    # Annihilation operators
    a_ops = []
    for i in range(n_modes):
        op_list = [qt.qeye(n_fock)] * n_modes
        op_list[i] = qt.destroy(n_fock)
        a_ops.append(qt.tensor(op_list))

    # Hamiltonian
    H = sum(HBAR * 2 * np.pi * frequencies[i] * a_ops[i].dag() * a_ops[i]
            for i in range(n_modes))

    # Add coupling terms
    for i in range(n_modes):
        for j in range(i+1, n_modes):
            if coupling[i, j] > 1e-10:
                H += HBAR * coupling[i, j] * (a_ops[i].dag() * a_ops[j] + a_ops[j].dag() * a_ops[i])

    # Collapse operators (thermal bath)
    c_ops = []
    gamma = 1e6  # Damping rate
    for i in range(n_modes):
        n_th = 1 / (np.exp(HBAR * frequencies[i] * 2 * np.pi / (KB * temperature)) - 1)
        # Decay
        c_ops.append(np.sqrt(gamma * (n_th + 1)) * a_ops[i])
        # Thermal excitation
        c_ops.append(np.sqrt(gamma * n_th) * a_ops[i].dag())

    # Initial state: coherent superposition in first mode
    psi0_list = [qt.basis(n_fock, 0)] * n_modes
    psi0_list[0] = (qt.basis(n_fock, 0) + qt.basis(n_fock, 1)).unit()
    psi0 = qt.tensor(psi0_list)

    # Time evolution
    times = np.linspace(0, t_max, 100)

    # Expectation of coherence (off-diagonal of first mode)
    coherence_op = a_ops[0]

    result = qt.mesolve(H, psi0, times, c_ops, [coherence_op])
    coherence = np.abs(result.expect[0])

    # Fit exponential decay to extract T2
    from scipy.optimize import curve_fit

    def decay(t, a, t2):
        return a * np.exp(-t / t2)

    try:
        popt, _ = curve_fit(decay, times, coherence, p0=[1.0, t_max/2], maxfev=10000)
        t2 = popt[1]
    except:
        # Fallback: find time to 1/e
        threshold = coherence[0] / np.e
        idx = np.argmax(coherence < threshold)
        t2 = times[idx] if idx > 0 else t_max

    return t2


def run_single_stage_experiment(n_modes: int = 7, temperature: float = 4.0,
                                 f0: float = 1e9) -> Dict:
    """
    Run single-stage comparison: Tonnetz vs Random coupling.
    """
    print(f"\n{'='*60}")
    print(f"Single Stage Experiment: {n_modes} modes at {temperature}K")
    print(f"{'='*60}")

    # Generate frequencies
    harmonic_freqs = generate_harmonic_frequencies(f0, n_modes)
    random_freqs = generate_random_frequencies(f0, n_modes)

    # Build coupling matrices
    tonnetz_coupling = build_coupling_matrix(harmonic_freqs, 'tonnetz')
    random_coupling = build_coupling_matrix(random_freqs, 'random')

    print("\nHarmonic frequencies (MHz):", harmonic_freqs / 1e6)
    print("Tonnetz coupling matrix:")
    print(np.round(tonnetz_coupling, 3))

    # Run simulations
    print("\nRunning Tonnetz simulation...")
    t2_tonnetz = run_qutip_simulation(harmonic_freqs, tonnetz_coupling, temperature)

    print("Running Random simulation...")
    t2_random = run_qutip_simulation(random_freqs, random_coupling, temperature)

    enhancement = (t2_tonnetz / t2_random - 1) * 100

    print(f"\n{'='*60}")
    print("RESULTS")
    print(f"{'='*60}")
    print(f"T2 (Tonnetz):  {t2_tonnetz*1e6:.3f} µs")
    print(f"T2 (Random):   {t2_random*1e6:.3f} µs")
    print(f"Enhancement:   {enhancement:.1f}%")
    print(f"{'='*60}")

    return {
        'n_modes': n_modes,
        'temperature': temperature,
        't2_tonnetz': t2_tonnetz,
        't2_random': t2_random,
        'enhancement_percent': enhancement
    }


def run_cascaded_experiment(n_stages: List[int] = [1, 2, 4, 8, 16],
                            n_modes: int = 7,
                            temperature: float = 0.02,  # 20 mK
                            f0: float = 1e9) -> Dict:
    """
    Run cascaded filtering experiment.
    Verify: T2(N) = T2(0) * (1 + epsilon)^N
    """
    print(f"\n{'='*60}")
    print(f"Cascaded Filtering Experiment")
    print(f"Modes per stage: {n_modes}, Temperature: {temperature*1000:.0f} mK")
    print(f"{'='*60}")

    results = []

    for n in n_stages:
        print(f"\nSimulating {n} cascaded stages...")

        # For cascaded simulation, we model the cumulative effect
        # Each stage provides multiplicative enhancement

        # First, get single-stage enhancement
        harmonic_freqs = generate_harmonic_frequencies(f0, n_modes)
        tonnetz_coupling = build_coupling_matrix(harmonic_freqs, 'tonnetz')

        t2_single = run_qutip_simulation(harmonic_freqs, tonnetz_coupling, temperature)

        # Baseline (no Tonnetz)
        random_freqs = generate_random_frequencies(f0, n_modes)
        random_coupling = build_coupling_matrix(random_freqs, 'random')
        t2_baseline = run_qutip_simulation(random_freqs, random_coupling, temperature)

        # Enhancement factor per stage
        epsilon = t2_single / t2_baseline - 1

        # Cascaded effect
        t2_cascaded = t2_baseline * ((1 + epsilon) ** n)

        results.append({
            'n_stages': n,
            't2_cascaded': t2_cascaded,
            'epsilon': epsilon
        })

        print(f"  Stages: {n}, T2: {t2_cascaded*1e6:.2f} µs ({t2_cascaded*1e3:.4f} ms)")

    # Fit exponential
    stages = np.array([r['n_stages'] for r in results])
    t2s = np.array([r['t2_cascaded'] for r in results])

    # Log-linear fit: log(T2) = log(T2_0) + N * log(1 + epsilon)
    log_t2 = np.log(t2s)
    coeffs = np.polyfit(stages, log_t2, 1)
    epsilon_fit = np.exp(coeffs[0]) - 1
    t2_0_fit = np.exp(coeffs[1])

    print(f"\n{'='*60}")
    print("CASCADED RESULTS")
    print(f"{'='*60}")
    print(f"Fitted epsilon: {epsilon_fit:.3f} ({epsilon_fit*100:.1f}% per stage)")
    print(f"Fitted T2(0):   {t2_0_fit*1e6:.3f} µs")
    print(f"\nProjected T2 for quantum repeater threshold (1 ms):")
    n_required = np.log(1e-3 / t2_0_fit) / np.log(1 + epsilon_fit)
    print(f"  Stages required: {n_required:.1f}")

    # Check if achievable
    if n_required <= 20:
        print(f"  ✓ ACHIEVABLE with ≤20 stages")
    else:
        print(f"  ✗ Requires >{20} stages")

    print(f"{'='*60}")

    return {
        'stages': stages.tolist(),
        't2_values': t2s.tolist(),
        'epsilon_fit': epsilon_fit,
        't2_0_fit': t2_0_fit,
        'n_required_for_1ms': n_required
    }


def main():
    parser = argparse.ArgumentParser(description='CNT Quantum Simulation')
    parser.add_argument('--modes', type=int, default=7, help='Number of oscillator modes')
    parser.add_argument('--stages', type=int, default=1, help='Number of filtering stages')
    parser.add_argument('--temperature', type=float, default=4.0, help='Temperature in Kelvin')
    parser.add_argument('--cascaded', action='store_true', help='Run cascaded experiment')
    parser.add_argument('--output', type=str, default=None, help='Output JSON file')
    args = parser.parse_args()

    print("="*60)
    print("CNT QUANTUM SIMULATION")
    print("Topological Coherence via Tonnetz Harmonic Filtering")
    print("="*60)
    print(f"QuTiP available: {QUTIP_AVAILABLE}")
    print(f"Date: {datetime.now().isoformat()}")

    if args.cascaded:
        results = run_cascaded_experiment(
            n_stages=[1, 2, 4, 8, 16],
            n_modes=args.modes,
            temperature=args.temperature
        )
    else:
        results = run_single_stage_experiment(
            n_modes=args.modes,
            temperature=args.temperature
        )

    if args.output:
        with open(args.output, 'w') as f:
            json.dump(results, f, indent=2)
        print(f"\nResults saved to {args.output}")

    return results


if __name__ == '__main__':
    main()