Exploring the Geometric Phase of Robinson Congruences in Quantum Circuits on IBM’s 127-Qubit Quantum Computer

Code:

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

# Logging
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

# Parse calibration data
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_2025-01-02T01_49_03Z.csv'
calibration_data = load_calibration_data(calibration_file)

# Select best qubits
num_qubits = 7  
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)

# Quantum and classical registers
qr = QuantumRegister(num_qubits, 'light_directions')
cr = ClassicalRegister(num_qubits, 'meas')
qc = QuantumCircuit(qr, cr)

# Initialize the quantum states representing light-ray directions
for i in range(num_qubits):
    qc.h(qr[i])  

# Vary a curvature parameter to encode geometric phase
curvature_param = 1.5  
for i in range(1, num_qubits): 
    qc.cp(curvature_param * (np.pi / (i + 1)), qr[0], qr[i])

# Apply an interferometric protocol for Berry phase measurement
qc.h(qr[0])  
for i in range(1, num_qubits):
    qc. cx(qr[0], qr[i])  

# Measure
qc.measure(qr, cr)

# Transpile 
logger. info("Transpiling the quantum circuit for the backend")
qc_transpiled = transpile(qc, backend=backend, optimization_level=3)
logger. info("Circuit transpilation complete")

# Execute 
shots = 16384
with Session(service=service, backend=backend) as session:
    sampler = SamplerV2(session=session)
    logger. info("Executing the circuit on the backend")
    job = sampler. run([qc_transpiled], shots=shots)
    job_result = job.result()

    # Extract BitArray counts
    data_bin = job_result._pub_results[0]['__value__']['data']
    bit_array = data_bin['meas']  # Extract the BitArray
    counts = bit_array.get_counts()  # Use get_counts method

    # Json
    results_data = {"raw_counts": counts}
    file_path = '/Users/Documents/Geometric_Phase_RC_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
plot_histogram(counts)
plt.title("Geometric Phase of Robinson Congruences")
plt. show()

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

Code for All Visuals
# Imports
import json
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# File path
file_path = '/Users/Documents/Geometric_Phase_RC_0.json'

# Load data
def load_data(file_path):
    with open(file_path, 'r') as f:
        data = json.load(f)
    return data["raw_counts"]

# Histogram of Measured States
def plot_histogram_measured_states(counts):
    total_counts = sum(counts.values())
    sorted_probs = sorted(counts.items(), key=lambda x: x[1], reverse=True)
    bitstrings, _ = zip(*sorted_probs)
    plt.figure(figsize=(12, 6))
    plt. bar(bitstrings[:50], [counts[b] for b in bitstrings[:50]], color='blue', alpha=0.7)
    plt.xticks(rotation=90, fontsize=8)
    plt.xlabel("Bitstrings")
    plt.ylabel("Counts")
    plt.title("Histogram of Measured States")
    plt.tight_layout()
    plt. show()

# Heatmap of Probabilities
def plot_heatmap_probabilities(probabilities):
    sorted_probs = sorted(probabilities.items(), key=lambda x: x[1], reverse=True)
    _, probs = zip(*sorted_probs)
    plt.figure(figsize=(12, 6))
    plt.imshow(np.array([list(probs)]), cmap="viridis", aspect="auto")
    plt.colorbar(label="Probability")
    plt.xlabel("Bitstring Index (sorted by probability)")
    plt.ylabel("Probability Density")
    plt.title("Heatmap of Probabilities")
    plt.tight_layout()
    plt. show()

# Dominant States
def plot_dominant_states(counts, top_n=10):
    sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)
    top_bitstrings, top_counts = zip(*sorted_counts[:top_n])
    plt.figure(figsize=(10, 5))
    plt. bar(top_bitstrings, top_counts, color='green', alpha=0.8)
    plt.xticks(rotation=45, fontsize=10)
    plt.xlabel("Top Bitstrings")
    plt.ylabel("Counts")
    plt.title(f"Top {top_n} Dominant States")
    plt.tight_layout()
    plt. show()

# Cumulative Distribution Function
def plot_cumulative_distribution(counts):
    total_counts = sum(counts.values())
    sorted_counts = np.cumsum([counts[b] for b in sorted(counts.keys())])
    sorted_counts_prob = sorted_counts / total_counts
    plt.figure(figsize=(10, 5))
    plt.plot(range(len(sorted_counts_prob)), sorted_counts_prob, color='red')
    plt.xlabel("Bitstring Index (sorted)")
    plt.ylabel("Cumulative Probability")
    plt.title("Cumulative Distribution Function (CDF)")
    plt.grid(True)
    plt.tight_layout()
    plt. show()

# Probability Distribution Across All States (Line Plot)
def plot_probability_distribution(probabilities):
    sorted_probs = sorted(probabilities.items(), key=lambda x: x[1], reverse=True)
    _, probs = zip(*sorted_probs)
    plt.figure(figsize=(12, 6))
    plt.plot(range(len(probs)), probs, marker='o', linestyle='-', color='blue', alpha=0.7)
    plt.xlabel("Bitstring Index (sorted by probability)")
    plt.ylabel("Probability")
    plt.title("Probability Distribution Across All States")
    plt.grid(True)
    plt.tight_layout()
    plt. show()

# Logarithmic Histogram of Measured States
def plot_logarithmic_histogram(counts):
    sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)
    bitstrings, _ = zip(*sorted_counts)
    plt.figure(figsize=(12, 6))
    plt. bar(bitstrings[:50], [np.log10(counts[b]) for b in bitstrings[:50]], color='orange', alpha=0.8)
    plt.xticks(rotation=90, fontsize=8)
    plt.xlabel("Bitstrings")
    plt.ylabel("Log10(Counts)")
    plt.title("Logarithmic Histogram of Measured States")
    plt.tight_layout()
    plt. show()

# State Transition Analysis (Scatter Plot)
def plot_state_transition_analysis(probabilities):
    sorted_probs = sorted(probabilities.items(), key=lambda x: x[1], reverse=True)
    _, probs = zip(*sorted_probs)
    plt.figure(figsize=(12, 6))
    plt.scatter(range(len(probs)), probs, color='green', alpha=0.6)
    plt.xlabel("Bitstring Index (sorted by probability)")
    plt.ylabel("Probability")
    plt.title("State Transition Analysis")
    plt.grid(True)
    plt.tight_layout()
    plt. show()

# Probability Deviation Analysis
def plot_probability_deviation_analysis(probabilities):
    bitstrings = list(probabilities.keys())
    probs = list(probabilities.values())
    n = len(probs)
    theoretical_phases = [np.cos(i * np.pi / n)**2 for i in range(n)]
    scale_factor = max(probs) / max(theoretical_phases)
    theoretical_phases_scaled = [p * scale_factor for p in theoretical_phases]
    deviations = [abs(p - t) for p, t in zip(probs, theoretical_phases_scaled)]
    plt.figure(figsize=(12, 6))
    plt.scatter(range(len(deviations)), deviations, color='purple', alpha=0.6)
    plt.xlabel("Bitstring Index")
    plt.ylabel("Deviation (|Measured - Theoretical|)")
    plt.title("Probability Deviation Analysis")
    plt.grid(True)
    plt.tight_layout()
    plt. show()

# Full Probability Heatmap
def plot_full_probability_heatmap(probabilities):
    sorted_probs = sorted(probabilities.items(), key=lambda x: x[0])
    bitstrings, probs = zip(*sorted_probs)
    heatmap_data = np.array(probs).reshape(1, -1)
    plt.figure(figsize=(12, 6))
    plt.imshow(heatmap_data, cmap="viridis", aspect="auto")
    plt.colorbar(label="Probability")
    plt.xlabel("Bitstring Index")
    plt.ylabel("Probability Density")
    plt.title("Full Probability Heatmap (Macro Overview)")
    plt.tight_layout()
    plt. show()

# Phase Contribution vs. Curvature vs. Counts
def plot_phase_contribution_vs_curvature(counts):
    sorted_counts = sorted(counts.items(), key=lambda x: x[0])
    bitstrings, counts_sorted = zip(*sorted_counts)
    curvatures = [1.5 * (i + 1) / len(bitstrings) for i in range(len(bitstrings))]
    phase_contributions = [np.cos(curv)**2 for curv in curvatures]
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection="3d")
    for idx, (curv, phase, count) in enumerate(zip(curvatures, phase_contributions, counts_sorted)):
        ax. bar([curv], [count], zs=phase, zdir="y", alpha=0.8, color=plt. cm.viridis(phase), edgecolor="k")
    ax.set_xlabel("Curvature Parameter")
    ax.set_ylabel("Phase Contribution")
    ax.set_zlabel("Measured Counts")
    ax.set_title("Phase Contribution vs. Curvature vs. Counts")
    plt.tight_layout()
    plt. show()

# 3D Surface Plot of Probabilities
def plot_3d_surface_probabilities(probabilities):
    sorted_probs = sorted(probabilities.items(), key=lambda x: x[0])
    bitstrings, probs = zip(*sorted_probs)
    indices = np.arange(len(probs))
    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(111, projection="3d")
    X, Y = np.meshgrid(indices, indices)
    Z = np.outer(probs, probs)
    ax.plot_surface(X, Y, Z, cmap="plasma", alpha=0.8)
    ax.set_xlabel("Bitstring Index (X)")
    ax.set_ylabel("Bitstring Index (Y)")
    ax.set_zlabel("Probability")
    ax.set_title("3D Surface Plot of Probabilities")
    plt.tight_layout()
    plt. show()

# Bitstring Pair Correlation (Scatter Plot)
def plot_bitstring_pair_correlation(counts):
    sorted_counts = sorted(counts.items(), key=lambda x: x[0])
    _, counts_sorted = zip(*sorted_counts)
    correlations = [counts_sorted[i] * counts_sorted[i + 1] for i in range(len(counts_sorted) - 1)]
    plt.figure(figsize=(12, 6))
    plt.scatter(range(len(correlations)), correlations, color="orange", alpha=0.6)
    plt.xlabel("Bitstring Pair Index")
    plt.ylabel("Correlation (Count[i] × Count[i+1])")
    plt.title("Bitstring Pair Correlation (Scatter Plot)")
    plt.grid(True)
    plt.tight_layout()
    plt. show()

# Main 
def main():
    counts = load_data(file_path)
    
    total_counts = sum(counts.values())
    probabilities = {k: v / total_counts for k, v in counts.items()}

    plot_histogram_measured_states(counts)
    plot_heatmap_probabilities(probabilities)
    plot_dominant_states(counts)
    plot_cumulative_distribution(counts)
    plot_probability_distribution(probabilities)
    plot_logarithmic_histogram(counts)
    plot_state_transition_analysis(probabilities)
    plot_probability_deviation_analysis(probabilities)
    plot_full_probability_heatmap(probabilities)
    plot_phase_contribution_vs_curvature(counts)
    plot_3d_surface_probabilities(probabilities)
    plot_bitstring_pair_correlation(counts)

if __name__ == "__main__":
    main()

# End