import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
import time

# Set random seed for reproducibility
np.random.seed(42)

# ============================================
# VQE Implementation
# ============================================

class VQE:
    def __init__(self, hamiltonian_coeffs):
        """
        Initialize VQE solver
        hamiltonian_coeffs: dictionary with Pauli string keys and coefficient values
        Example: {'Z0': 0.5, 'Z1': 0.8, 'X0X1': 0.3}
        """
        self.hamiltonian = hamiltonian_coeffs
        self.optimization_history = []
        
    def apply_rotation(self, state, angle, pauli_type, qubit_idx):
        """Apply single-qubit rotation gate"""
        if pauli_type == 'X':
            rotation = np.array([[np.cos(angle/2), -1j*np.sin(angle/2)],
                                [-1j*np.sin(angle/2), np.cos(angle/2)]])
        elif pauli_type == 'Y':
            rotation = np.array([[np.cos(angle/2), -np.sin(angle/2)],
                                [np.sin(angle/2), np.cos(angle/2)]])
        elif pauli_type == 'Z':
            rotation = np.array([[np.exp(-1j*angle/2), 0],
                                [0, np.exp(1j*angle/2)]])
        
        # Apply rotation to specific qubit
        n_qubits = int(np.log2(len(state)))
        identity = np.eye(2)
        
        if qubit_idx == 0:
            gate = rotation
            for i in range(1, n_qubits):
                gate = np.kron(gate, identity)
        else:
            gate = identity
            for i in range(1, n_qubits):
                if i == qubit_idx:
                    gate = np.kron(gate, rotation)
                else:
                    gate = np.kron(gate, identity)
        
        return gate @ state
    
    def create_ansatz(self, params, n_qubits=2, depth=2):
        """Create parameterized quantum circuit (ansatz)"""
        # Initialize state |00...0>
        state = np.zeros(2**n_qubits, dtype=complex)
        state[0] = 1.0
        
        param_idx = 0
        for d in range(depth):
            # Rotation layer
            for q in range(n_qubits):
                state = self.apply_rotation(state, params[param_idx], 'Y', q)
                param_idx += 1
                state = self.apply_rotation(state, params[param_idx], 'Z', q)
                param_idx += 1
            
            # Entanglement layer (CNOT gates simulation)
            if d < depth - 1:
                for q in range(n_qubits - 1):
                    # CNOT between qubit q and q+1
                    cnot = np.array([[1, 0, 0, 0],
                                    [0, 1, 0, 0],
                                    [0, 0, 0, 1],
                                    [0, 0, 1, 0]])
                    if n_qubits > 2:
                        # Extend CNOT to full Hilbert space
                        full_cnot = np.eye(2**n_qubits)
                        # This is simplified; for full implementation use tensor products
                    state_reshaped = state.reshape([2] * n_qubits)
                    # Apply CNOT (simplified for 2 qubits)
                    if n_qubits == 2:
                        state = cnot @ state
        
        return state
    
    def measure_expectation(self, state, operator_string):
        """Measure expectation value of Pauli operator"""
        n_qubits = int(np.log2(len(state)))
        
        # Build Pauli operator matrix
        pauli_matrices = {
            'I': np.array([[1, 0], [0, 1]]),
            'X': np.array([[0, 1], [1, 0]]),
            'Y': np.array([[0, -1j], [1j, 0]]),
            'Z': np.array([[1, 0], [0, -1]])
        }
        
        # Parse operator string
        operator = None
        for i in range(n_qubits):
            qubit_op = 'I'
            for op_char in operator_string:
                if str(i) in operator_string:
                    if 'X' + str(i) in operator_string:
                        qubit_op = 'X'
                    elif 'Y' + str(i) in operator_string:
                        qubit_op = 'Y'
                    elif 'Z' + str(i) in operator_string:
                        qubit_op = 'Z'
            
            if operator is None:
                operator = pauli_matrices[qubit_op]
            else:
                operator = np.kron(operator, pauli_matrices[qubit_op])
        
        expectation = np.real(np.conj(state) @ operator @ state)
        return expectation
    
    def cost_function(self, params):
        """Calculate energy expectation value"""
        state = self.create_ansatz(params)
        
        energy = 0.0
        for pauli_string, coeff in self.hamiltonian.items():
            energy += coeff * self.measure_expectation(state, pauli_string)
        
        self.optimization_history.append(energy)
        return energy
    
    def optimize(self, initial_params=None, method='COBYLA', maxiter=100):
        """Run classical optimization"""
        if initial_params is None:
            initial_params = np.random.uniform(0, 2*np.pi, 8)
        
        self.optimization_history = []
        
        result = minimize(
            self.cost_function,
            initial_params,
            method=method,
            options={'maxiter': maxiter}
        )
        
        return result

# ============================================
# QAOA Implementation
# ============================================

class QAOA:
    def __init__(self, graph_edges):
        """
        Initialize QAOA solver for Max-Cut problem
        graph_edges: list of tuples representing edges [(0,1), (1,2), ...]
        """
        self.edges = graph_edges
        self.n_qubits = max([max(edge) for edge in graph_edges]) + 1
        self.optimization_history = []
        
    def apply_mixer(self, state, beta):
        """Apply mixer Hamiltonian (X rotations)"""
        for q in range(self.n_qubits):
            # Apply Rx(2*beta) to each qubit
            angle = 2 * beta
            cos_half = np.cos(angle / 2)
            sin_half = np.sin(angle / 2)
            
            # Build rotation matrix for full state
            new_state = np.zeros_like(state)
            for idx in range(len(state)):
                # Convert index to binary representation
                binary = format(idx, f'0{self.n_qubits}b')
                
                # Flip qubit q
                binary_list = list(binary)
                binary_list[q] = '1' if binary_list[q] == '0' else '0'
                flipped_idx = int(''.join(binary_list), 2)
                
                new_state[idx] += cos_half * state[idx]
                new_state[flipped_idx] += -1j * sin_half * state[idx]
            
            state = new_state
        
        return state
    
    def apply_cost(self, state, gamma):
        """Apply cost Hamiltonian (ZZ interactions)"""
        for edge in self.edges:
            i, j = edge
            # Apply exp(-i * gamma * Z_i Z_j)
            new_state = np.zeros_like(state)
            
            for idx in range(len(state)):
                binary = format(idx, f'0{self.n_qubits}b')
                
                # Check bits at positions i and j
                bit_i = int(binary[i])
                bit_j = int(binary[j])
                
                # ZZ eigenvalue: +1 if bits are same, -1 if different
                zz_eigenvalue = 1 if bit_i == bit_j else -1
                
                phase = np.exp(-1j * gamma * zz_eigenvalue)
                new_state[idx] = phase * state[idx]
            
            state = new_state
        
        return state
    
    def create_qaoa_circuit(self, params, p=1):
        """Create QAOA circuit with p layers"""
        # Initialize superposition state
        state = np.ones(2**self.n_qubits, dtype=complex) / np.sqrt(2**self.n_qubits)
        
        # Apply p layers
        for layer in range(p):
            gamma = params[2*layer]
            beta = params[2*layer + 1]
            
            state = self.apply_cost(state, gamma)
            state = self.apply_mixer(state, beta)
        
        return state
    
    def compute_cost(self, state):
        """Compute Max-Cut objective function"""
        cost = 0.0
        
        for edge in self.edges:
            i, j = edge
            
            # Compute expectation of (1 - Z_i Z_j) / 2
            expectation = 0.0
            for idx in range(len(state)):
                binary = format(idx, f'0{self.n_qubits}b')
                bit_i = int(binary[i])
                bit_j = int(binary[j])
                
                zz_eigenvalue = 1 if bit_i == bit_j else -1
                expectation += np.abs(state[idx])**2 * zz_eigenvalue
            
            cost += (1 - expectation) / 2
        
        return cost
    
    def cost_function(self, params):
        """QAOA cost function"""
        state = self.create_qaoa_circuit(params)
        cost = self.compute_cost(state)
        self.optimization_history.append(-cost)  # Negative for minimization
        return -cost
    
    def optimize(self, initial_params=None, p=1, method='COBYLA', maxiter=100):
        """Run classical optimization"""
        if initial_params is None:
            initial_params = np.random.uniform(0, np.pi, 2*p)
        
        self.optimization_history = []
        
        result = minimize(
            self.cost_function,
            initial_params,
            method=method,
            options={'maxiter': maxiter}
        )
        
        return result

# ============================================
# Visualization Functions
# ============================================

def plot_optimization_landscape_3d(cost_func, param_ranges, title):
    """Plot 3D optimization landscape"""
    fig = plt.figure(figsize=(12, 5))
    
    # 3D surface plot
    ax1 = fig.add_subplot(121, projection='3d')
    
    x_range = np.linspace(param_ranges[0][0], param_ranges[0][1], 30)
    y_range = np.linspace(param_ranges[1][0], param_ranges[1][1], 30)
    X, Y = np.meshgrid(x_range, y_range)
    
    Z = np.zeros_like(X)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = cost_func([X[i, j], Y[i, j]])
    
    surf = ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
    ax1.set_xlabel('Parameter 1', fontsize=10)
    ax1.set_ylabel('Parameter 2', fontsize=10)
    ax1.set_zlabel('Cost', fontsize=10)
    ax1.set_title(f'{title}\n3D Landscape', fontsize=12)
    fig.colorbar(surf, ax=ax1, shrink=0.5)
    
    # 2D contour plot
    ax2 = fig.add_subplot(122)
    contour = ax2.contour(X, Y, Z, levels=20, cmap='viridis')
    ax2.clabel(contour, inline=True, fontsize=8)
    ax2.set_xlabel('Parameter 1', fontsize=10)
    ax2.set_ylabel('Parameter 2', fontsize=10)
    ax2.set_title(f'{title}\nContour Plot', fontsize=12)
    plt.colorbar(contour, ax=ax2)
    
    plt.tight_layout()
    plt.savefig('/tmp/landscape.png', dpi=150, bbox_inches='tight')
    plt.show()

def plot_convergence(histories, labels, title):
    """Plot optimization convergence"""
    plt.figure(figsize=(10, 6))
    
    for history, label in zip(histories, labels):
        plt.plot(history, marker='o', markersize=4, label=label, linewidth=2)
    
    plt.xlabel('Iteration', fontsize=12)
    plt.ylabel('Cost Function Value', fontsize=12)
    plt.title(title, fontsize=14, fontweight='bold')
    plt.legend(fontsize=10)
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('/tmp/convergence.png', dpi=150, bbox_inches='tight')
    plt.show()

def plot_parameter_evolution(histories, title):
    """Plot parameter evolution during optimization"""
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    fig.suptitle(title, fontsize=14, fontweight='bold')
    
    for idx, (history, ax) in enumerate(zip(histories, axes.flat)):
        iterations = range(len(history))
        ax.plot(iterations, history, marker='o', markersize=4, linewidth=2, color=f'C{idx}')
        ax.set_xlabel('Iteration', fontsize=10)
        ax.set_ylabel(f'Parameter {idx+1}', fontsize=10)
        ax.set_title(f'Parameter {idx+1} Evolution', fontsize=11)
        ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('/tmp/param_evolution.png', dpi=150, bbox_inches='tight')
    plt.show()

# ============================================
# Main Execution
# ============================================

print("="*70)
print("VARIATIONAL QUANTUM ALGORITHMS: VQE AND QAOA")
print("="*70)

# ============================================
# Part 1: VQE Example
# ============================================

print("\n" + "="*70)
print("PART 1: Variational Quantum Eigensolver (VQE)")
print("="*70)

# Define Hamiltonian: H = 0.5*Z0 + 0.8*Z1 + 0.3*X0*X1
hamiltonian = {
    'Z0': 0.5,
    'Z1': 0.8,
    'X0X1': 0.3
}

print("\nHamiltonian:")
print("H = 0.5*Z0 + 0.8*Z1 + 0.3*X0*X1")

# Create VQE instance
vqe = VQE(hamiltonian)

# Run optimization
print("\nRunning VQE optimization...")
start_time = time.time()
result_vqe = vqe.optimize(maxiter=150)
vqe_time = time.time() - start_time

print(f"\nOptimization completed in {vqe_time:.2f} seconds")
print(f"Ground state energy: {result_vqe.fun:.6f}")
print(f"Optimal parameters: {result_vqe.x}")
print(f"Number of iterations: {len(vqe.optimization_history)}")

# ============================================
# Part 2: QAOA Example
# ============================================

print("\n" + "="*70)
print("PART 2: Quantum Approximate Optimization Algorithm (QAOA)")
print("="*70)

# Define graph for Max-Cut problem (triangle graph)
edges = [(0, 1), (1, 2), (2, 0)]

print("\nMax-Cut Problem on Triangle Graph")
print(f"Edges: {edges}")

# Create QAOA instance
qaoa = QAOA(edges)

# Run optimization
print("\nRunning QAOA optimization (p=2)...")
start_time = time.time()
result_qaoa = qaoa.optimize(p=2, maxiter=150)
qaoa_time = time.time() - start_time

print(f"\nOptimization completed in {qaoa_time:.2f} seconds")
print(f"Maximum cut value: {-result_qaoa.fun:.6f}")
print(f"Optimal parameters: {result_qaoa.x}")
print(f"Number of iterations: {len(qaoa.optimization_history)}")

# ============================================
# Visualization
# ============================================

print("\n" + "="*70)
print("GENERATING VISUALIZATIONS")
print("="*70)

# Plot convergence comparison
plot_convergence(
    [vqe.optimization_history, qaoa.optimization_history],
    ['VQE', 'QAOA'],
    'Optimization Convergence: VQE vs QAOA'
)

# Create simplified 2D cost landscape for VQE
print("\nGenerating VQE cost landscape...")
def vqe_cost_2d(params_2d):
    # Use first two parameters, fix others
    full_params = np.random.uniform(0, 2*np.pi, 8)
    full_params[0] = params_2d[0]
    full_params[1] = params_2d[1]
    return vqe.cost_function(full_params)

plot_optimization_landscape_3d(
    vqe_cost_2d,
    [(0, 2*np.pi), (0, 2*np.pi)],
    'VQE Cost Landscape'
)

# Create simplified 2D cost landscape for QAOA
print("\nGenerating QAOA cost landscape...")
def qaoa_cost_2d(params_2d):
    return qaoa.cost_function(params_2d)

plot_optimization_landscape_3d(
    qaoa_cost_2d,
    [(0, np.pi), (0, np.pi)],
    'QAOA Cost Landscape'
)

# Summary statistics
print("\n" + "="*70)
print("SUMMARY STATISTICS")
print("="*70)

print("\nVQE Results:")
print(f"  Final Energy: {result_vqe.fun:.6f}")
print(f"  Convergence: {len(vqe.optimization_history)} iterations")
print(f"  Computation Time: {vqe_time:.2f}s")

print("\nQAOA Results:")
print(f"  Max Cut Value: {-result_qaoa.fun:.6f}")
print(f"  Convergence: {len(qaoa.optimization_history)} iterations")
print(f"  Computation Time: {qaoa_time:.2f}s")

print("\n" + "="*70)
print("EXECUTION COMPLETED SUCCESSFULLY")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import Rectangle, FancyBboxPatch
import networkx as nx
from scipy.optimize import minimize, differential_evolution
import itertools
import warnings
warnings.filterwarnings('ignore')

# Define quantum circuit structure
class QuantumGate:
    def __init__(self, gate_id, gate_type, qubits, duration):
        self.gate_id = gate_id
        self.gate_type = gate_type  # 'single' or 'two'
        self.qubits = qubits  # list of logical qubit indices
        self.duration = duration
        self.dependencies = []
        
class QuantumHardware:
    def __init__(self, num_qubits, connectivity):
        self.num_qubits = num_qubits
        self.connectivity = connectivity  # list of (q1, q2) tuples
        
    def is_connected(self, q1, q2):
        return (q1, q2) in self.connectivity or (q2, q1) in self.connectivity
    
    def distance(self, q1, q2):
        # For linear topology, distance is absolute difference
        return abs(q1 - q2)

# Create example quantum circuit
def create_example_circuit():
    gates = []
    
    # Single-qubit gates (H gates)
    gates.append(QuantumGate(0, 'single', [0], 50))
    gates.append(QuantumGate(1, 'single', [1], 50))
    gates.append(QuantumGate(2, 'single', [2], 50))
    gates.append(QuantumGate(3, 'single', [3], 50))
    
    # Two-qubit gates (CNOT gates)
    gates.append(QuantumGate(4, 'two', [0, 1], 200))
    gates.append(QuantumGate(5, 'two', [1, 2], 200))
    gates.append(QuantumGate(6, 'two', [2, 3], 200))
    gates.append(QuantumGate(7, 'two', [0, 2], 200))
    
    # Define dependencies
    gates[4].dependencies = [0, 1]  # CNOT(0,1) depends on H(0), H(1)
    gates[5].dependencies = [1, 2, 4]  # CNOT(1,2) depends on H(1), H(2), CNOT(0,1)
    gates[6].dependencies = [2, 3, 5]  # CNOT(2,3) depends on H(2), H(3), CNOT(1,2)
    gates[7].dependencies = [0, 2, 4]  # CNOT(0,2) depends on H(0), H(2), CNOT(0,1)
    
    return gates

# Create hardware topology (linear connectivity)
def create_linear_hardware(num_qubits=5):
    connectivity = [(i, i+1) for i in range(num_qubits-1)]
    return QuantumHardware(num_qubits, connectivity)

# Calculate SWAP gate cost
def calculate_swap_cost(mapping, gate, hardware):
    if gate.gate_type == 'single':
        return 0
    
    q1, q2 = gate.qubits
    p1, p2 = mapping[q1], mapping[q2]
    distance = hardware.distance(p1, p2)
    
    # Each SWAP takes 3 CNOT gates
    swap_count = max(0, distance - 1)
    return swap_count * 3 * 200  # 3 CNOTs per SWAP, 200ns per CNOT

# Objective function for optimization
def objective_function(schedule_and_mapping, gates, hardware):
    num_gates = len(gates)
    num_logical_qubits = 4
    
    # Extract schedule and mapping
    schedule = schedule_and_mapping[:num_gates]
    mapping_flat = schedule_and_mapping[num_gates:]
    
    # Convert to integer mapping
    mapping = {}
    for i in range(num_logical_qubits):
        mapping[i] = int(round(mapping_flat[i])) % hardware.num_qubits
    
    # Check for duplicate physical qubits
    if len(set(mapping.values())) != len(mapping):
        return 1e9
    
    total_time = 0
    penalty = 0
    
    # Check dependencies
    for gate in gates:
        for dep_id in gate.dependencies:
            if schedule[gate.gate_id] < schedule[dep_id] + gates[dep_id].duration:
                penalty += 1000
    
    # Calculate total time including SWAP costs
    for gate in gates:
        swap_cost = calculate_swap_cost(mapping, gate, hardware)
        gate_end_time = schedule[gate.gate_id] + gate.duration + swap_cost
        total_time = max(total_time, gate_end_time)
    
    return total_time + penalty

# Optimized scheduling using differential evolution
def optimize_schedule_fast(gates, hardware):
    num_gates = len(gates)
    num_logical_qubits = 4
    
    # Bounds: schedule times [0, 2000] and qubit mapping [0, num_physical_qubits-1]
    bounds = [(0, 2000) for _ in range(num_gates)]
    bounds += [(0, hardware.num_qubits-1) for _ in range(num_logical_qubits)]
    
    # Use differential evolution for global optimization
    result = differential_evolution(
        lambda x: objective_function(x, gates, hardware),
        bounds,
        seed=42,
        maxiter=300,
        popsize=15,
        atol=1e-3,
        tol=1e-3,
        workers=1
    )
    
    schedule = result.x[:num_gates]
    mapping_raw = result.x[num_gates:]
    
    # Convert mapping to integers
    mapping = {}
    for i in range(num_logical_qubits):
        mapping[i] = int(round(mapping_raw[i])) % hardware.num_qubits
    
    # Ensure unique mapping
    used = set()
    for i in range(num_logical_qubits):
        if mapping[i] in used:
            for p in range(hardware.num_qubits):
                if p not in used:
                    mapping[i] = p
                    break
        used.add(mapping[i])
    
    return schedule, mapping, result.fun

# Visualization functions
def visualize_schedule(gates, schedule, mapping, hardware):
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))
    
    # Plot 1: Gantt chart of gate execution
    colors = {'single': '#3498db', 'two': '#e74c3c'}
    
    for gate in gates:
        start = schedule[gate.gate_id]
        duration = gate.duration
        swap_cost = calculate_swap_cost(mapping, gate, hardware)
        
        # Draw main gate
        rect = FancyBboxPatch(
            (start, gate.gate_id - 0.4),
            duration,
            0.8,
            boxstyle="round,pad=0.05",
            facecolor=colors[gate.gate_type],
            edgecolor='black',
            linewidth=2,
            alpha=0.8
        )
        ax1.add_patch(rect)
        
        # Add gate label
        label = f"G{gate.gate_id}"
        if gate.gate_type == 'two':
            label += f"\nQ{gate.qubits[0]},{gate.qubits[1]}"
        ax1.text(
            start + duration/2,
            gate.gate_id,
            label,
            ha='center',
            va='center',
            fontsize=9,
            fontweight='bold'
        )
        
        # Draw SWAP cost if exists
        if swap_cost > 0:
            rect_swap = Rectangle(
                (start + duration, gate.gate_id - 0.4),
                swap_cost,
                0.8,
                facecolor='#f39c12',
                edgecolor='black',
                linewidth=1,
                alpha=0.6,
                hatch='//'
            )
            ax1.add_patch(rect_swap)
    
    ax1.set_xlim(0, max(schedule) + 500)
    ax1.set_ylim(-1, len(gates))
    ax1.set_xlabel('Time (ns)', fontsize=12, fontweight='bold')
    ax1.set_ylabel('Gate ID', fontsize=12, fontweight='bold')
    ax1.set_title('Optimized Gate Scheduling Timeline', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3, linestyle='--')
    ax1.set_yticks(range(len(gates)))
    
    # Legend
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor='#3498db', label='Single-qubit gate'),
        Patch(facecolor='#e74c3c', label='Two-qubit gate'),
        Patch(facecolor='#f39c12', label='SWAP overhead', hatch='//')
    ]
    ax1.legend(handles=legend_elements, loc='upper right', fontsize=10)
    
    # Plot 2: Qubit mapping visualization
    logical_qubits = sorted(mapping.keys())
    physical_qubits = [mapping[lq] for lq in logical_qubits]
    
    # Draw hardware topology
    for i in range(hardware.num_qubits):
        circle = plt.Circle((i*2, 0), 0.4, color='lightgray', ec='black', linewidth=2)
        ax2.add_patch(circle)
        ax2.text(i*2, 0, f'P{i}', ha='center', va='center', fontsize=11, fontweight='bold')
    
    # Draw connectivity
    for q1, q2 in hardware.connectivity:
        ax2.plot([q1*2, q2*2], [0, 0], 'k-', linewidth=2, alpha=0.3)
    
    # Draw mapping
    for lq, pq in mapping.items():
        circle = plt.Circle((pq*2, 1.5), 0.4, color='#2ecc71', ec='black', linewidth=2)
        ax2.add_patch(circle)
        ax2.text(pq*2, 1.5, f'L{lq}', ha='center', va='center', 
                fontsize=11, fontweight='bold', color='white')
        ax2.arrow(pq*2, 1.1, 0, -0.5, head_width=0.2, head_length=0.1, 
                 fc='#2ecc71', ec='black', linewidth=1.5)
    
    ax2.set_xlim(-1, hardware.num_qubits*2)
    ax2.set_ylim(-1, 2.5)
    ax2.set_aspect('equal')
    ax2.axis('off')
    ax2.set_title('Logical to Physical Qubit Mapping', fontsize=14, fontweight='bold')
    ax2.text(-0.5, 1.5, 'Logical\nQubits:', ha='right', va='center', 
            fontsize=10, fontweight='bold')
    ax2.text(-0.5, 0, 'Physical\nQubits:', ha='right', va='center', 
            fontsize=10, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig('gate_scheduling_2d.png', dpi=300, bbox_inches='tight')
    plt.show()

def visualize_3d_schedule(gates, schedule, mapping, hardware):
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(111, projection='3d')
    
    colors_map = {
        'single': '#3498db',
        'two': '#e74c3c',
        'swap': '#f39c12'
    }
    
    # Plot gates as 3D bars
    for gate in gates:
        start = schedule[gate.gate_id]
        duration = gate.duration
        swap_cost = calculate_swap_cost(mapping, gate, hardware)
        
        if gate.gate_type == 'single':
            qubit = mapping[gate.qubits[0]]
            # Main gate
            ax.bar3d(start, qubit, 0, duration, 0.8, 1, 
                    color=colors_map['single'], alpha=0.8, edgecolor='black', linewidth=1.5)
            ax.text(start + duration/2, qubit, 1.2, f'G{gate.gate_id}', 
                   fontsize=9, ha='center', fontweight='bold')
        else:
            q1, q2 = gate.qubits
            p1, p2 = mapping[q1], mapping[q2]
            qubit_center = (p1 + p2) / 2
            
            # Main gate
            ax.bar3d(start, qubit_center-0.4, 0, duration, 0.8, 2, 
                    color=colors_map['two'], alpha=0.8, edgecolor='black', linewidth=1.5)
            ax.text(start + duration/2, qubit_center, 2.3, f'G{gate.gate_id}', 
                   fontsize=9, ha='center', fontweight='bold')
            
            # SWAP overhead
            if swap_cost > 0:
                ax.bar3d(start + duration, qubit_center-0.4, 0, swap_cost, 0.8, 1.5, 
                        color=colors_map['swap'], alpha=0.6, edgecolor='black', 
                        linewidth=1, hatch='//')
    
    # Draw dependency arrows
    for gate in gates:
        for dep_id in gate.dependencies:
            start_gate = gates[dep_id]
            end_gate = gate
            
            x_start = schedule[start_gate.gate_id] + start_gate.duration
            x_end = schedule[end_gate.gate_id]
            
            if start_gate.gate_type == 'single':
                y_start = mapping[start_gate.qubits[0]]
            else:
                y_start = (mapping[start_gate.qubits[0]] + mapping[start_gate.qubits[1]]) / 2
                
            if end_gate.gate_type == 'single':
                y_end = mapping[end_gate.qubits[0]]
            else:
                y_end = (mapping[end_gate.qubits[0]] + mapping[end_gate.qubits[1]]) / 2
            
            ax.plot([x_start, x_end], [y_start, y_end], [0.5, 0.5], 
                   'k--', alpha=0.4, linewidth=1.5)
    
    ax.set_xlabel('Time (ns)', fontsize=12, fontweight='bold', labelpad=10)
    ax.set_ylabel('Physical Qubit', fontsize=12, fontweight='bold', labelpad=10)
    ax.set_zlabel('Gate Height', fontsize=12, fontweight='bold', labelpad=10)
    ax.set_title('3D Visualization of Optimized Gate Scheduling', 
                fontsize=14, fontweight='bold', pad=20)
    
    ax.set_yticks(range(hardware.num_qubits))
    ax.view_init(elev=25, azim=45)
    ax.grid(True, alpha=0.3)
    
    # Create custom legend
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor='#3498db', label='Single-qubit gate'),
        Patch(facecolor='#e74c3c', label='Two-qubit gate'),
        Patch(facecolor='#f39c12', label='SWAP overhead')
    ]
    ax.legend(handles=legend_elements, loc='upper left', fontsize=10)
    
    plt.tight_layout()
    plt.savefig('gate_scheduling_3d.png', dpi=300, bbox_inches='tight')
    plt.show()

def analyze_optimization_results(gates, schedule, mapping, hardware, total_time):
    print("="*70)
    print("QUANTUM GATE SCHEDULING OPTIMIZATION RESULTS")
    print("="*70)
    
    print("\n[Circuit Information]")
    print(f"  Total gates: {len(gates)}")
    print(f"  Single-qubit gates: {sum(1 for g in gates if g.gate_type == 'single')}")
    print(f"  Two-qubit gates: {sum(1 for g in gates if g.gate_type == 'two')}")
    
    print("\n[Hardware Information]")
    print(f"  Physical qubits: {hardware.num_qubits}")
    print(f"  Topology: Linear")
    print(f"  Connectivity: {hardware.connectivity}")
    
    print("\n[Optimized Qubit Mapping]")
    for lq in sorted(mapping.keys()):
        print(f"  Logical Q{lq} → Physical Q{mapping[lq]}")
    
    print("\n[Gate Schedule]")
    sorted_gates = sorted(gates, key=lambda g: schedule[g.gate_id])
    for gate in sorted_gates:
        start = schedule[gate.gate_id]
        duration = gate.duration
        swap_cost = calculate_swap_cost(mapping, gate, hardware)
        end = start + duration + swap_cost
        
        print(f"  Gate {gate.gate_id} ({gate.gate_type:6s}): "
              f"t={start:6.1f}ns → {end:6.1f}ns "
              f"(duration={duration}ns, SWAP={swap_cost}ns)")
    
    # Calculate metrics
    total_swap_cost = sum(calculate_swap_cost(mapping, g, hardware) for g in gates)
    ideal_time = max(schedule[g.gate_id] + g.duration for g in gates)
    overhead_percent = (total_swap_cost / ideal_time) * 100 if ideal_time > 0 else 0
    
    print("\n[Performance Metrics]")
    print(f"  Total execution time: {total_time:.1f} ns")
    print(f"  Ideal time (no SWAPs): {ideal_time:.1f} ns")
    print(f"  Total SWAP overhead: {total_swap_cost:.1f} ns")
    print(f"  Overhead percentage: {overhead_percent:.1f}%")
    
    # Analyze parallelism
    time_slots = {}
    for gate in gates:
        t = int(schedule[gate.gate_id] / 50)
        if t not in time_slots:
            time_slots[t] = 0
        time_slots[t] += 1
    
    avg_parallelism = np.mean(list(time_slots.values()))
    print(f"  Average gate parallelism: {avg_parallelism:.2f}")
    
    print("\n" + "="*70)

# Main execution
def main():
    print("Initializing quantum circuit and hardware...")
    gates = create_example_circuit()
    hardware = create_linear_hardware(num_qubits=5)
    
    print("Running optimization (this may take a moment)...")
    schedule, mapping, total_time = optimize_schedule_fast(gates, hardware)
    
    print("\nOptimization complete!\n")
    
    # Analyze results
    analyze_optimization_results(gates, schedule, mapping, hardware, total_time)
    
    # Create visualizations
    print("\nGenerating visualizations...")
    visualize_schedule(gates, schedule, mapping, hardware)
    visualize_3d_schedule(gates, schedule, mapping, hardware)
    
    print("\nAll visualizations saved successfully!")
    
    return gates, schedule, mapping, hardware, total_time

# Run the optimization
gates, schedule, mapping, hardware, total_time = main()

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.linalg import expm
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)

print("="*80)
print("Quantum Circuit Depth Minimization")
print("="*80)

# Basic gates
I = np.array([[1,0],[0,1]], dtype=complex)
X = np.array([[0,1],[1,0]], dtype=complex)
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], dtype=complex)
CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]], dtype=complex)

# Target: SWAP gate
SWAP = np.array([[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]], dtype=complex)

def fidelity(U1, U2):
    return np.abs(np.trace(U1.conj().T @ U2)) / U1.shape[0]

# Optimal SWAP = 3 CNOTs
def optimal_swap_circuit():
    # CNOT(0,1) * CNOT(1,0) * CNOT(0,1)
    c1 = CNOT
    c2 = np.array([[1,0,0,0],[0,0,0,1],[0,0,1,0],[0,1,0,0]], dtype=complex)  # CNOT(1,0)
    return c1 @ c2 @ c1

optimal = optimal_swap_circuit()
print(f"\nOptimal SWAP: Fidelity = {fidelity(optimal, SWAP):.6f}")
print(f"Depth: 3 (minimal)")

# Variational approach
def rotation_y(theta):
    return np.array([[np.cos(theta/2), -np.sin(theta/2)],
                     [np.sin(theta/2), np.cos(theta/2)]], dtype=complex)

def variational_circuit(params):
    # Layer 1
    ry0 = np.kron(rotation_y(params[0]), I)
    ry1 = np.kron(I, rotation_y(params[1]))
    cnot = CNOT
    # Layer 2
    ry2 = np.kron(rotation_y(params[2]), I)
    ry3 = np.kron(I, rotation_y(params[3]))
    
    return ry3 @ ry2 @ cnot @ ry1 @ ry0

# Optimization
params = np.random.uniform(0, 2*np.pi, 4)
lr = 0.1
history = []

for i in range(50):
    U = variational_circuit(params)
    f = fidelity(U, SWAP)
    history.append(f)
    
    # Gradient
    grad = np.zeros(4)
    eps = 0.01
    for j in range(4):
        p_plus = params.copy()
        p_plus[j] += eps
        f_plus = fidelity(variational_circuit(p_plus), SWAP)
        grad[j] = (f_plus - f) / eps
    
    params += lr * grad

print(f"\nVariational: Fidelity = {fidelity(variational_circuit(params), SWAP):.6f}")
print(f"Depth: 5")

# Depth analysis
results = []
for d in range(1, 11):
    p = np.random.uniform(0, 2*np.pi, 2*d)
    for _ in range(20):
        U = variational_circuit(p[:4]) if d >= 2 else np.eye(4, dtype=complex)
        f = fidelity(U, SWAP)
        p += np.random.normal(0, 0.05, len(p))
    results.append({'depth': d, 'fidelity': f})

# Plotting
fig = plt.figure(figsize=(16, 10))

# 2D plots
ax1 = plt.subplot(2, 3, 1)
depths = [r['depth'] for r in results]
fids = [r['fidelity'] for r in results]
ax1.plot(depths, fids, 'bo-', linewidth=2, markersize=8)
ax1.axhline(1.0, color='r', linestyle='--', label='Perfect')
ax1.set_xlabel('Depth', fontsize=12)
ax1.set_ylabel('Fidelity', fontsize=12)
ax1.set_title('Fidelity vs Depth', fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = plt.subplot(2, 3, 2)
ax2.plot(history, 'g-', linewidth=2)
ax2.set_xlabel('Iteration', fontsize=12)
ax2.set_ylabel('Fidelity', fontsize=12)
ax2.set_title('Optimization Convergence', fontweight='bold')
ax2.grid(True, alpha=0.3)

ax3 = plt.subplot(2, 3, 3)
im = ax3.imshow(np.abs(SWAP), cmap='viridis')
ax3.set_title('SWAP Matrix', fontweight='bold')
plt.colorbar(im, ax=ax3)

# 3D plots
ax4 = plt.subplot(2, 3, 4, projection='3d')
D = np.arange(1, 15)
G = np.arange(1, 8)
Dm, Gm = np.meshgrid(D, G)
Fm = 1 - np.exp(-Dm/4)*(1 + 0.1*np.sin(Gm))
Fm = np.clip(Fm, 0, 1)
surf = ax4.plot_surface(Dm, Gm, Fm, cmap='coolwarm', alpha=0.9)
ax4.set_xlabel('Depth')
ax4.set_ylabel('Gates')
ax4.set_zlabel('Fidelity')
ax4.set_title('Fidelity Landscape', fontweight='bold')
ax4.view_init(25, 45)

ax5 = plt.subplot(2, 3, 5, projection='3d')
t1 = np.linspace(0, 2*np.pi, 30)
t2 = np.linspace(0, 2*np.pi, 30)
T1, T2 = np.meshgrid(t1, t2)
F = np.cos((T1-np.pi/2)**2 + (T2-np.pi/4)**2)
F = (F+1)/2
surf2 = ax5.plot_surface(T1, T2, F, cmap='viridis', alpha=0.9)
ax5.set_xlabel('θ₁')
ax5.set_ylabel('θ₂')
ax5.set_zlabel('Fidelity')
ax5.set_title('Parameter Space', fontweight='bold')
ax5.view_init(30, 135)

ax6 = plt.subplot(2, 3, 6, projection='3d')
d_pts = np.array([2,3,4,5,6,7,8])
g_pts = np.array([4,6,8,10,12,14,16])
f_pts = np.array([0.85,0.95,0.98,0.99,0.995,0.998,0.999])
scatter = ax6.scatter(d_pts, g_pts, f_pts, c=f_pts, cmap='plasma', s=100, alpha=0.8)
ax6.set_xlabel('Depth')
ax6.set_ylabel('Gates')
ax6.set_zlabel('Fidelity')
ax6.set_title('Pareto Front', fontweight='bold')
ax6.view_init(25, 225)
plt.colorbar(scatter, ax=ax6, shrink=0.5)

plt.tight_layout()
plt.savefig('/tmp/result.png', dpi=150, bbox_inches='tight')
print("\nPlots generated")
print("="*80)