Null-Light Ray Tomography with a Twistor Inspired Circuit on IBM’s 127-Qubit Quantum Computer

Code:

# Imports
import numpy as np
import pandas as pd
import json
import logging
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit_ibm_runtime import QiskitRuntimeService, Session, SamplerV2
from qiskit.visualization import plot_histogram, plot_bloch_multivector
from qiskit.quantum_info import Statevector
import matplotlib.pyplot as plt
import seaborn as sns
from qiskit.quantum_info import DensityMatrix

# Logging setup
logging.basicConfig(level=logging. INFO, format='%(asctime)s - %(levelname)s - %(message)s')
logger = logging.getLogger(__name__)

# Load calibration data
def load_calibration_data(file_path):
    logger. info("Loading calibration data from %s", file_path)
    calibration_data = pd. read_csv(file_path)
    calibration_data.columns = calibration_data.columns.str.strip()
    logger. info("Calibration data loaded successfully")
    return calibration_data

# Select the best qubits based on T1, T2, and error rates
def select_best_qubits(calibration_data, n_qubits):
    logger. info("Selecting the best qubits based on T1, T2, and error rates")
    qubits_sorted = calibration_data.sort_values(by=['\u221ax (sx) error', 'T1 (us)', 'T2 (us)'], ascending=[True, False, False])
    best_qubits = qubits_sorted['Qubit'].head(n_qubits).tolist()
    logger. info("Selected qubits: %s", best_qubits)
    return best_qubits

# Load backend calibration data
calibration_file = '/Users/Downloads/ibm_brisbane_calibrations_2024-12-30T17_39_45Z.csv'
calibration_data = load_calibration_data(calibration_file)

# Select best qubits
num_qubits = 4  
best_qubits = select_best_qubits(calibration_data, num_qubits)

# IBMQ 
logger. info("Setting up IBM Q service")
service = QiskitRuntimeService(
    channel='ibm_quantum',
    instance='ibm-q/open/main',
    token='YOUR_IBMQ_API_KEY_O-`'  
)

backend_name = 'ibm_brisbane'
backend = service.backend(backend_name)
logger. info("Backend selected: %s", backend_name)

# Parametric null-ray preparation circuit 
def prepare_null_ray_circuit(theta, phi):
    qr = QuantumRegister(num_qubits, 'q')
    cr = ClassicalRegister(num_qubits, 'c')
    qc = QuantumCircuit(qr, cr)

    # Prepare initial state using sx and rz (decomposition of u gate)
    for i in range(0, num_qubits, 2):
        qc.rz(phi, qr[i])
        qc. sx(qr[i])
        qc.rz(theta, qr[i])

        qc.rz(-phi, qr[i + 1])
        qc. sx(qr[i + 1])
        qc.rz(-theta, qr[i + 1])

    # Decomposition with rz, sx, and x
    for i in range(0, num_qubits, 2):
        qc.rz(np.pi / 2, qr[i + 1])
        qc. sx(qr[i + 1])
        qc.rz(-np.pi / 2, qr[i + 1])
        qc. sx(qr[i])
        qc.rz(np.pi / 2, qr[i])
        qc. sx(qr[i])
        qc.rz(-np.pi / 2, qr[i])
        qc. sx(qr[i + 1])
        qc.rz(np.pi / 2, qr[i + 1])
        qc. sx(qr[i + 1])
        qc.rz(-np.pi / 2, qr[i + 1])

    return qc

# Generate tomography circuits
basis_gates = ['rz', 'sx', 'x', 'id', 'measure']
angles = np.linspace(0, 2 * np.pi, 6)  
circuits = []
for theta in angles:
    for phi in angles:
        qc = prepare_null_ray_circuit(theta, phi)
        qc.measure_all()
        transpiled_qc = transpile(qc, backend=backend, basis_gates=basis_gates, optimization_level=3)
        circuits.append(transpiled_qc)

# Execute 
shots = 16384
logger. info("Starting quantum tomography on %d circuits", len(circuits))
with Session(service=service, backend=backend) as session:
    sampler = SamplerV2(session=session)
    job = sampler. run(circuits, shots=shots)
    job_result = job.result()

    # Results
    results_data = []
    for idx, circuit_result in enumerate(job_result._pub_results):
        counts = circuit_result["__value__"]["data"].meas.get_counts()
        results_data.append({"circuit_idx": idx, "counts": counts})

    # Json
    file_path = '/Users/Documents/NullRay_Twister_Result_0.json'
    with open(file_path, 'w') as f:
        json.dump(results_data, f, indent=4)
    logger. info("Results saved to %s", file_path)

# Visual
for result in results_data:
    counts = result["counts"]
    plot_histogram(counts)
    plt.title(f"Circuit {result['circuit_idx']} Measurement Results")
    plt. show()

/////////////////

Code for All Visuals
# Imports
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import json
from qiskit.quantum_info import DensityMatrix, Statevector
from qiskit.visualization import plot_bloch_multivector
from scipy.stats import entropy

# Load run data
file_path = '/Users/Documents/NullRay_Twister_Result_0.json'
try:
    with open(file_path, 'r') as file:
        data = json.load(file)
except FileNotFoundError:
    raise FileNotFoundError(f"The file at {file_path} was not found. Please check the path.")

# Labels
labels = sorted(data[0]['counts'].keys())

# Global variables 
dominant_states = ["0111", "1111"]

# Measurement Distribution Heatmap
def plot_heatmap():
    heatmap_data = []
    for circuit in data:
        row = [circuit['counts'].get(label, 0) for label in labels]
        heatmap_data.append(row)

    plt.figure(figsize=(15, 10))
    sns.heatmap(np.array(heatmap_data), xticklabels=labels, yticklabels=[f"Circuit {c['circuit_idx']}" for c in data], cmap='viridis', cbar=True)
    plt.title("Measurement Distribution Heatmap", fontsize=16)
    plt.xlabel("States", fontsize=14)
    plt.ylabel("Circuits", fontsize=14)
    plt. show()

# Bloch Sphere Projections for Selected Circuits
def plot_bloch_spheres():
    selected_indices = [0, 5, 10, 15]
    for idx in selected_indices:
        counts = data[idx]['counts']
        total = sum(counts.values())
        probs = [np.sqrt(counts.get(state, 0) / total) for state in labels]  

        # Create Statevector approximation
        sv = Statevector(probs)
        dm = DensityMatrix(sv)

        plot_bloch_multivector(dm)
        plt.title(f"Bloch Sphere Projection: Circuit {idx}", fontsize=16)
        plt. show()

# Interference Strength Across Circuits
def plot_interference_strength():
    interference_strength = []

    for circuit in data:
        total_counts = sum(circuit['counts'].values())
        dominant_counts = sum(circuit['counts'].get(state, 0) for state in dominant_states)
        interference_strength.append(dominant_counts / total_counts)

    plt.figure(figsize=(10, 6))
    plt.plot(range(len(data)), interference_strength, marker='o', color='b')
    plt.title("Interference Strength Across Circuits", fontsize=16)
    plt.xlabel("Circuit Index", fontsize=14)
    plt.ylabel("Interference Strength (Dominant State Proportion)", fontsize=14)
    plt.grid(True)
    plt. show()

# Error Rate Analysis
def plot_error_rate():
    expected_distribution = [1/16] * 16  
    def calculate_error(observed, expected):
        return sum(abs(o - e) for o, e in zip(observed, expected))

    def normalize_distribution(counts, labels):
        total = sum(counts.values())
        return [counts.get(label, 0) / total for label in labels]

    errors = []
    for circuit in data:
        observed_distribution = normalize_distribution(circuit['counts'], labels)
        errors.append(calculate_error(observed_distribution, expected_distribution))

    plt.figure(figsize=(10, 6))
    plt.plot(range(len(data)), errors, marker='x', color='r')
    plt.title("Error Rate Analysis Across Circuits", fontsize=16)
    plt.xlabel("Circuit Index", fontsize=14)
    plt.ylabel("Total Error (Observed vs Expected)", fontsize=14)
    plt.grid(True)
    plt. show()

# State Population Histogram
def plot_state_population():
    plt.figure(figsize=(15, 8))
    state_populations = {label: 0 for label in labels}

    for circuit in data:
        for state, count in circuit['counts'].items():
            state_populations[state] += count

    state_counts = [state_populations[state] for state in labels]
    plt. bar(labels, state_counts, color='blue')
    plt.title("Overall State Population Histogram", fontsize=16)
    plt.xlabel("States", fontsize=14)
    plt.ylabel("Total Counts", fontsize=14)
    plt.xticks(rotation=90)
    plt. show()

# Cumulative Interference Strength
def plot_cumulative_interference():
    cumulative_interference = []
    cumulative_sum = 0

    for circuit in data:
        total_counts = sum(circuit['counts'].values())
        dominant_counts = sum(circuit['counts'].get(state, 0) for state in dominant_states)
        cumulative_sum += dominant_counts / total_counts
        cumulative_interference.append(cumulative_sum)

    plt.figure(figsize=(10, 6))
    plt.plot(range(len(data)), cumulative_interference, marker='o', color='green')
    plt.title("Cumulative Interference Strength", fontsize=16)
    plt.xlabel("Circuit Index", fontsize=14)
    plt.ylabel("Cumulative Strength", fontsize=14)
    plt.grid(True)
    plt. show()

# State Entropy Across Circuits
def plot_state_entropy():
    entropies = []

    for circuit in data:
        total_counts = sum(circuit['counts'].values())
        probabilities = [circuit['counts'].get(state, 0) / total_counts for state in labels]
        entropies.append(entropy(probabilities))

    plt.figure(figsize=(10, 6))
    plt.plot(range(len(data)), entropies, marker='x', color='orange')
    plt.title("State Entropy Across Circuits", fontsize=16)
    plt.xlabel("Circuit Index", fontsize=14)
    plt.ylabel("Entropy (bits)", fontsize=14)
    plt.grid(True)
    plt. show()

# Dominant State Evolution
def plot_dominant_state_evolution():
    dominant_state_counts = {state: [] for state in dominant_states}

    for circuit in data:
        total_counts = sum(circuit['counts'].values())
        for state in dominant_states:
            dominant_state_counts[state].append(circuit['counts'].get(state, 0) / total_counts)

    plt.figure(figsize=(10, 6))
    for state, counts in dominant_state_counts.items():
        plt.plot(range(len(data)), counts, label=f"State {state}")

    plt.title("Dominant State Evolution", fontsize=16)
    plt.xlabel("Circuit Index", fontsize=14)
    plt.ylabel("Proportion of Total Counts", fontsize=14)
    plt.legend()
    plt.grid(True)
    plt. show()

# Call
plot_heatmap()
plot_bloch_spheres()
plot_interference_strength()
plot_error_rate()
plot_state_population()
plot_cumulative_interference()
plot_state_entropy()
plot_dominant_state_evolution()

# End.