Process Tomography and Characterization

def qpt(U: Qobj, op_basis: list = None) -> Qobj:
    """
    Calculate quantum process tomography chi-matrix from unitary.
    
    Parameters:
    - U: Unitary operator or superoperator to characterize
    - op_basis: Operator basis for process representation (default: Pauli basis)
    
    Returns:
    - Qobj: Chi-matrix representation of quantum process
    """

def qpt_plot(chi: Qobj, lbls_list: list = None, title: str = None, 
             fig: object = None, threshold: float = None) -> None:
    """
    Plot quantum process tomography chi-matrix.
    
    Parameters:
    - chi: Chi-matrix from QPT
    - lbls_list: Labels for process basis elements
    - title: Plot title
    - fig: Matplotlib figure object
    - threshold: Minimum value to display
    """

def qpt_plot_combined(chi: Qobj, lbls_list: list = None, title: str = None,
                      fig: object = None, threshold: float = 0.01) -> None:
    """
    Plot combined real and imaginary parts of chi-matrix.
    
    Parameters:
    - chi: Chi-matrix from QPT
    - lbls_list: Labels for basis elements
    - title: Plot title
    - fig: Matplotlib figure
    - threshold: Display threshold for matrix elements
    """

def process_fidelity(chi1: Qobj, chi2: Qobj) -> float:
    """
    Calculate fidelity between two quantum processes.
    
    Parameters:
    - chi1: First process chi-matrix
    - chi2: Second process chi-matrix
    
    Returns:
    - float: Process fidelity (0 to 1)
    """

def average_gate_fidelity(chi: Qobj, target_U: Qobj = None) -> float:
    """
    Calculate average gate fidelity of quantum process.
    
    Parameters:
    - chi: Process chi-matrix
    - target_U: Target unitary (identity if None)
    
    Returns:
    - float: Average gate fidelity
    """

def chi_to_pauli_transfer(chi: Qobj) -> np.ndarray:
    """
    Convert chi-matrix to Pauli transfer matrix representation.
    
    Parameters:
    - chi: Chi-matrix in Pauli basis
    
    Returns:
    - ndarray: Pauli transfer matrix
    """

def pauli_transfer_to_chi(R: np.ndarray) -> Qobj:
    """
    Convert Pauli transfer matrix to chi-matrix.
    
    Parameters:
    - R: Pauli transfer matrix
    
    Returns:
    - Qobj: Chi-matrix representation
    """

import qutip as qt
import numpy as np
import matplotlib.pyplot as plt

# Single-qubit process tomography examples

# 1. Perfect gates
print("=== Perfect Gate Process Tomography ===")

# Identity gate
I_gate = qt.qeye(2)
chi_I = qt.qpt(I_gate)
print(f"Identity gate chi-matrix:")
print(f"  Shape: {chi_I.shape}")
print(f"  Real part diagonal: {np.diag(chi_I.full().real)}")

# Pauli-X gate
X_gate = qt.sigmax()
chi_X = qt.qpt(X_gate)

# Pauli-Y gate  
Y_gate = qt.sigmay()
chi_Y = qt.qpt(Y_gate)

# Pauli-Z gate
Z_gate = qt.sigmaz()
chi_Z = qt.qpt(Z_gate)

# Hadamard gate
H_gate = (qt.sigmax() + qt.sigmaz()).unit()
chi_H = qt.qpt(H_gate)

# Process fidelities with perfect gates
perfect_gates = {
    'I': (I_gate, chi_I),
    'X': (X_gate, chi_X), 
    'Y': (Y_gate, chi_Y),
    'Z': (Z_gate, chi_Z),
    'H': (H_gate, chi_H)
}

print("Perfect gate process fidelities (should be 1.0):")
for name, (gate, chi) in perfect_gates.items():
    # Compare with itself
    fid = qt.process_fidelity(chi, chi)
    print(f"  {name} gate: {fid:.6f}")

# Cross-fidelities between different gates
print("Cross-fidelities between different perfect gates:")
gate_names = list(perfect_gates.keys())
for i, name1 in enumerate(gate_names):
    for j, name2 in enumerate(gate_names):
        if i < j:  # Only upper triangle
            chi1 = perfect_gates[name1][1]
            chi2 = perfect_gates[name2][1]
            fid = qt.process_fidelity(chi1, chi2)
            print(f"  {name1}-{name2}: {fid:.3f}")

# 2. Noisy gates (depolarizing noise)
print("\n=== Noisy Gate Process Tomography ===")

def depolarizing_noise(U, p):
    """Add depolarizing noise to unitary gate."""
    # Kraus operators for depolarizing channel after unitary
    I = qt.qeye(2)
    X = qt.sigmax()
    Y = qt.sigmay()
    Z = qt.sigmaz()
    
    K0 = np.sqrt(1 - 3*p/4) * U
    K1 = np.sqrt(p/4) * X * U
    K2 = np.sqrt(p/4) * Y * U
    K3 = np.sqrt(p/4) * Z * U
    
    return [K0, K1, K2, K3]

# Noisy X gate
p_noise = 0.1  # 10% depolarizing noise
noisy_X_kraus = depolarizing_noise(X_gate, p_noise)

# Convert Kraus to superoperator for QPT
def kraus_to_super(kraus_ops):
    """Convert Kraus operators to superoperator."""
    N = kraus_ops[0].shape[0]
    S = qt.Qobj(np.zeros((N**2, N**2)), dims=[[[N],[N]], [[N],[N]]])
    
    for K in kraus_ops:
        S += qt.sprepost(K, K.dag())
    
    return S

S_noisy_X = kraus_to_super(noisy_X_kraus)
chi_noisy_X = qt.qpt(S_noisy_X)

# Compare noisy vs perfect X gate
fid_noisy = qt.process_fidelity(chi_X, chi_noisy_X)
print(f"Noisy X gate (p={p_noise}) vs perfect X gate fidelity: {fid_noisy:.3f}")

# Average gate fidelity
avg_fid = qt.average_gate_fidelity(chi_noisy_X, X_gate)
print(f"Average gate fidelity of noisy X gate: {avg_fid:.3f}")

# 3. Two-qubit process tomography
print("\n=== Two-Qubit Process Tomography ===")

# CNOT gate
CNOT = qt.cnot()
chi_CNOT = qt.qpt(CNOT)
print(f"CNOT chi-matrix shape: {chi_CNOT.shape}")

# Controlled-Z gate
CZ = qt.cz_gate()
chi_CZ = qt.qpt(CZ)

# Process fidelity between CNOT and CZ
fid_CNOT_CZ = qt.process_fidelity(chi_CNOT, chi_CZ)
print(f"CNOT vs CZ process fidelity: {fid_CNOT_CZ:.3f}")

# 4. Amplitude damping channel
print("\n=== Amplitude Damping Channel ===")

def amplitude_damping_kraus(gamma):
    """Amplitude damping Kraus operators."""
    K0 = qt.Qobj([[1, 0], [0, np.sqrt(1-gamma)]])
    K1 = qt.Qobj([[0, np.sqrt(gamma)], [0, 0]])
    return [K0, K1]

gamma = 0.2
AD_kraus = amplitude_damping_kraus(gamma)
S_AD = kraus_to_super(AD_kraus)
chi_AD = qt.qpt(S_AD)

print(f"Amplitude damping (γ={gamma}) chi-matrix computed")

# Test on different input states
test_states = [
    qt.basis(2, 0),      # |0⟩
    qt.basis(2, 1),      # |1⟩  
    (qt.basis(2,0) + qt.basis(2,1)).unit()  # |+⟩
]

print("Amplitude damping effect on different states:")
for i, psi in enumerate(test_states):
    rho_in = psi * psi.dag()
    
    # Apply channel using Kraus operators
    rho_out = sum(K * rho_in * K.dag() for K in AD_kraus)
    
    # Calculate excited state population
    pop_in = qt.expect(qt.num(2), rho_in)
    pop_out = qt.expect(qt.num(2), rho_out)
    
    state_names = ['|0⟩', '|1⟩', '|+⟩']
    print(f"  {state_names[i]}: population {pop_in:.1f} → {pop_out:.3f}")

# 5. Process tomography visualization
print("\n=== Process Visualization ===")

# Generate Pauli basis labels for single qubit
pauli_labels = ['I', 'X', 'Y', 'Z']

print("Chi-matrix elements for different processes:")

# Perfect X gate
print("Perfect X gate (should have chi_XX = 1):")
chi_X_full = chi_X.full()
for i, label_i in enumerate(pauli_labels):
    for j, label_j in enumerate(pauli_labels):
        element = chi_X_full[i, j]
        if abs(element) > 0.01:  # Only show significant elements
            print(f"  χ_{label_i}{label_j} = {element:.3f}")

# Noisy X gate
print(f"\nNoisy X gate (p={p_noise}):")
chi_noisy_X_full = chi_noisy_X.full()
for i, label_i in enumerate(pauli_labels):
    for j, label_j in enumerate(pauli_labels):
        element = chi_noisy_X_full[i, j]
        if abs(element) > 0.01:
            print(f"  χ_{label_i}{label_j} = {element:.3f}")

# 6. Pauli transfer matrix representation
print("\n=== Pauli Transfer Matrix ===")

# Convert chi-matrix to Pauli transfer matrix
R_X = qt.chi_to_pauli_transfer(chi_X)
print("Perfect X gate Pauli transfer matrix:")
print("(Maps Pauli expectation values: ⟨σᵢ⟩_out = Σⱼ Rᵢⱼ ⟨σⱼ⟩_in)")

pauli_ops = [qt.qeye(2), qt.sigmax(), qt.sigmay(), qt.sigmaz()]
for i, op_i in enumerate(pauli_ops):
    for j, op_j in enumerate(pauli_ops):
        if abs(R_X[i, j]) > 0.01:
            print(f"  R_{pauli_labels[i]}{pauli_labels[j]} = {R_X[i, j]:6.3f}")

# Verify: X gate should map σz → -σz, σy → -σy, σx → σx, I → I
print("Expected: X maps σz→-σz, σy→-σy, σx→σx, I→I")

# Convert back to chi-matrix
chi_X_reconstructed = qt.pauli_transfer_to_chi(R_X)
reconstruction_fidelity = qt.process_fidelity(chi_X, chi_X_reconstructed)
print(f"PTM reconstruction fidelity: {reconstruction_fidelity:.6f}")

# 7. Process tomography from experimental data (simulation)
print("\n=== Experimental Process Tomography Simulation ===")

def simulate_qpt_experiment(true_process_kraus, n_measurements=1000):
    """
    Simulate QPT experiment with finite measurement statistics.
    """
    # Input states for process tomography (Pauli eigenstates)
    input_states = [
        qt.basis(2, 0),  # |0⟩ (σz eigenstate)
        qt.basis(2, 1),  # |1⟩ (σz eigenstate) 
        (qt.basis(2,0) + qt.basis(2,1)).unit(),  # |+⟩ (σx eigenstate)
        (qt.basis(2,0) - qt.basis(2,1)).unit(),  # |-⟩ (σx eigenstate)
        (qt.basis(2,0) + 1j*qt.basis(2,1)).unit(),  # |+i⟩ (σy eigenstate)
        (qt.basis(2,0) - 1j*qt.basis(2,1)).unit(),  # |-i⟩ (σy eigenstate)
    ]
    
    # Measurement operators (Pauli measurements)
    measurement_ops = [qt.qeye(2), qt.sigmax(), qt.sigmay(), qt.sigmaz()]
    
    print(f"Simulating QPT with {n_measurements} measurements per configuration...")
    
    # For each input state, apply process and measure
    results = []
    for psi_in in input_states:
        rho_in = psi_in * psi_in.dag()
        
        # Apply process
        rho_out = sum(K * rho_in * K.dag() for K in true_process_kraus)
        
        # Measure expectation values
        expectations = []
        for M in measurement_ops:
            # True expectation value
            exp_true = qt.expect(M, rho_out)
            
            # Add statistical noise (simulate finite measurements)
            measurement_noise = np.random.normal(0, 1/np.sqrt(n_measurements))
            exp_measured = exp_true + measurement_noise
            
            expectations.append(exp_measured)
        
        results.append(expectations)
    
    return results

# Simulate experiment for noisy X gate
experimental_data = simulate_qpt_experiment(noisy_X_kraus, n_measurements=5000)

print("Experimental QPT simulation completed")
print("(In real experiment, this data would be used to reconstruct chi-matrix)")
print("Sample expectation values for first input state:")
for i, exp_val in enumerate(experimental_data[0]):
    print(f"  ⟨{pauli_labels[i]}⟩ = {exp_val:.3f}")

# 8. Gate set tomography concept (simplified)
print("\n=== Gate Set Analysis ===")

# Analyze a set of gates for process tomography
gate_set = {
    'I': qt.qeye(2),
    'X': qt.sigmax(),
    'Y': qt.sigmay(), 
    'Z': qt.sigmaz(),
    'H': (qt.sigmax() + qt.sigmaz()).unit(),
    'S': qt.phasegate(np.pi/2)
}

print("Gate set process fidelity matrix:")
gate_names = list(gate_set.keys())
chi_matrices = {name: qt.qpt(gate) for name, gate in gate_set.items()}

# Create fidelity matrix
fidelity_matrix = np.zeros((len(gate_names), len(gate_names)))
for i, name1 in enumerate(gate_names):
    for j, name2 in enumerate(gate_names):
        fid = qt.process_fidelity(chi_matrices[name1], chi_matrices[name2])
        fidelity_matrix[i, j] = fid

print("Process fidelity matrix (1.0 = identical, 0.0 = orthogonal):")
print("    ", " ".join(f"{name:^6}" for name in gate_names))
for i, name in enumerate(gate_names):
    row_str = f"{name:>3} "
    row_str += " ".join(f"{fidelity_matrix[i,j]:6.3f}" for j in range(len(gate_names)))
    print(row_str)

# Process tomography functions return:
# - qpt(): Qobj representing chi-matrix
# - process_fidelity(): float (0 to 1)  
# - average_gate_fidelity(): float (0 to 1)
# - chi_to_pauli_transfer(): numpy ndarray
# - Plotting functions return None (display plots)