Approximating the Speed of Light from 4-bit E=mc² Interference on a 133-Qubit Quantum Computer

Code:

# Main circuit
# Imports
import logging, json, math
import numpy as np
import pandas as pd
from math import gcd, isfinite
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit.circuit.library import UnitaryGate, QFT
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2

# IBMQ 
TOKEN = "YOUR_IBMQ_API_KEY "
INSTANCE = "YOUR_IBMQ_CRN"
BACKEND = "ibm_torino"
CAL_CSV = "/Users/steventippeconnic/Downloads/ibm_torino_calibrations_2025-10-04T19_20_42Z.csv"
SHOTS   = 8192

# 4-bit modulus and layout
ORDER = 16
N_Q   = 4
N_Q_TOTAL = N_Q * 3  # a, b, p

# Physical inputs (independent E and m)
object_mass_kg      = 1.0
object_energy_joule = 8.98755179e16  # independent E 

# Ticks / scaling 
DM_PER_TICK = 1.0      # kg per tick  (Δm)
DE_PER_TICK = 1.0e16   # J  per tick  (ΔE) ~1e16 J so E_idx ≈ 9 mod 16

# Speed search window
MIN_SPEED = 1.0e8
MAX_SPEED = 4.0e8
MAX_CANDIDATES_PER_K = 3

# Logging
logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
log = logging.getLogger("E_mc2_4bit")

# Choose best physical qubits from calibration CSV
def best_qubits(csv_path: str, n: int) -> list[int]:
    df = pd.read_csv(csv_path)
    df.columns = df.columns.str.strip()
    winners = (
        df.sort_values(["√x (sx) error", "T1 (us)", "T2 (us)"],
                       ascending=[True, False, False])
        ["Qubit"].head(n).tolist()
    )
    log.info("Best physical qubits: %s", winners)
    return winners

PHYSICAL = best_qubits(CAL_CSV, N_Q_TOTAL)

# Encode m and E into Z_N
def encode_indices(m_kg: float, E_joule: float, dm: float, dE: float, mod: int):
    """
    Returns (M_IDX, E_IDX, C2_UNIT) with:
      M_IDX = round(m/dm) mod N
      E_IDX = round(E/dE) mod N
      C2_UNIT = dE/dm  (so c² ≈ (k + tN)*C2_UNIT)
    """
    M_IDX = int(round(m_kg / dm)) % mod
    E_IDX = int(round(E_joule / dE)) % mod
    C2_UNIT = dE / dm
    if M_IDX == 0 or E_IDX == 0:
        log.warning("Encoded M_IDX=%d, E_IDX=%d; adjust Δm/ΔE to avoid zeros mod %d.", M_IDX, E_IDX, mod)
    return M_IDX, E_IDX, C2_UNIT

M_IDX, E_IDX, C2_UNIT = encode_indices(object_mass_kg, object_energy_joule, DM_PER_TICK, DE_PER_TICK, ORDER)
log.info("Encoded → M_IDX=%d, E_IDX=%d (mod %d).  c² unit scale ΔE/Δm = %.6g J/kg",
         M_IDX, E_IDX, ORDER, C2_UNIT)

# Modular constant adders
def add_const_modN_gate(c: int, mod: int) -> UnitaryGate:
    dim = mod
    mat = np.zeros((dim, dim))
    for x in range(dim):
        mat[(x + c) % mod, x] = 1.0
    return UnitaryGate(mat, label=f"+{c}_mod{mod}")

ADDERS = {c: add_const_modN_gate(c, ORDER) for c in range(1, ORDER)}  # c=0 is identity

def controlled_add(qc: QuantumCircuit, ctrl_qubit, point_reg, constant):
    qc.append(ADDERS[constant].control(), [ctrl_qubit, *point_reg])

# Oracle: f(a,b) = a*E_IDX + b*M_IDX (mod 16)
def emc_oracle(qc: QuantumCircuit, a_reg, b_reg, p_reg):
    # add a * E
    for i in range(N_Q):
        constant = (E_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, a_reg[i], p_reg, constant)
    # add b * m
    for i in range(N_Q):
        constant = (M_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, b_reg[i], p_reg, constant)

# Full circuit
def e_mc2_circuit() -> QuantumCircuit:
    a = QuantumRegister(N_Q, "a")
    b = QuantumRegister(N_Q, "b")
    p = QuantumRegister(N_Q, "p")
    c = ClassicalRegister(N_Q * 2, "c")
    qc = QuantumCircuit(a, b, p, c, name="E_equals_m_c2_mod16")

    qc.h(a)
    qc.h(b)
    emc_oracle(qc, a, b, p)
    qc.barrier()

    qc.append(QFT(N_Q, do_swaps=False), a)
    qc.append(QFT(N_Q, do_swaps=False), b)

    qc.measure(a, c[:N_Q])
    qc.measure(b, c[N_Q:])
    return qc

# Run on IBM Runtime v2
service = QiskitRuntimeService(channel="ibm_cloud", token=TOKEN, instance=INSTANCE)
backend = service.backend(BACKEND)
log.info("Backend → %s", backend.name)

qc_raw = e_mc2_circuit()
trans = transpile(qc_raw, backend=backend, initial_layout=PHYSICAL, optimization_level=3)
log.info("Circuit depth %d, gate counts %s", trans.depth(), trans.count_ops())

sampler = SamplerV2(mode=backend)
job = sampler.run([trans], shots=SHOTS)
result = job.result()

# Post-processing
creg_name = trans.cregs[0].name
counts_raw = result[0].data.__getattribute__(creg_name).get_counts()

def bits_to_int(bs: str) -> int:
    return int(bs[::-1], 2)  # reverse endian to match convention

# Map raw "b a" bitstrings -> (a_val, b_val)
counts_pairs = {(bits_to_int(k[N_Q:]), bits_to_int(k[:N_Q])): v for k, v in counts_raw.items()}
top_pairs = sorted(counts_pairs.items(), key=lambda kv: kv[1], reverse=True)

# Recover k from top invertible b rows: k ≡ -u*v^{-1} (mod N)
top_invertibles = []
for (a_val, b_val), freq in top_pairs:
    if gcd(b_val, ORDER) != 1:
        continue
    inv_b = pow(b_val, -1, ORDER)
    k_mod = (-a_val * inv_b) % ORDER
    top_invertibles.append(((a_val, b_val), k_mod, freq))
    if len(top_invertibles) == 100:
        break

if not top_invertibles:
    print("\nWARNING - No invertible b rows found in top measurements.\n")

print("\nTop invertible (a,b) → recovered k ≡ E·m^{-1} (mod 16):")
for (a_val, b_val), k_mod, freq in top_invertibles[:32]:
    print(f"  (a={a_val:2d}, b={b_val:2d})  →  k = {k_mod:2d}   (count = {freq})")

# Lift residues to c (m/s)
def c_candidates_from_k(k_mod: int,
                        order: int,
                        c2_unit: float,
                        vmin: float,
                        vmax: float,
                        max_k: int = 3):
    """
    Return a few (t, c) with c² = (k_mod + t*order) * c2_unit inside [vmin², vmax²].
    """
    c2_min = vmin * vmin
    c2_max = vmax * vmax
    t_low  = math.ceil((c2_min / c2_unit - k_mod) / order)
    t_high = math.floor((c2_max / c2_unit - k_mod) / order)
    if t_low > t_high:
        return []
    picks = [t_low, (t_low + t_high) // 2, t_high]
    seen, out = set(), []
    mid = 0.5*(vmin + vmax)
    for t in picks:
        if t in seen:
            continue
        seen.add(t)
        c2 = (k_mod + t*order) * c2_unit
        if c2 <= 0 or not isfinite(c2):
            continue
        c_val = math.sqrt(c2)
        if vmin <= c_val <= vmax:
            out.append((t, c_val))
    out.sort(key=lambda x: abs(x[1] - mid))
    return out[:max_k]

# Aggregate by k value (heaviest first)
k_by_weight = {}
for (_, _), k_mod, freq in top_invertibles:
    k_by_weight[k_mod] = k_by_weight.get(k_mod, 0) + freq

unique_k_sorted = sorted(k_by_weight.items(), key=lambda kv: kv[1], reverse=True)

print("\nCandidate speed-of-light estimates (lifted from residues):")
print(f"  scale ΔE/Δm = {C2_UNIT:.6g} J/kg, window [{MIN_SPEED:.2e}, {MAX_SPEED:.2e}] m/s")
printed = 0
for k_mod, total_w in unique_k_sorted[:8]:
    cands = c_candidates_from_k(k_mod, ORDER, C2_UNIT, MIN_SPEED, MAX_SPEED, MAX_CANDIDATES_PER_K)
    if not cands:
        continue
    print(f"  k = {k_mod:2d} (weight={total_w}) → ", end="")
    pretty = [f"t={t}: c≈{c:0.6e} m/s" for (t, c) in cands]
    print(";  ".join(pretty))
    printed += 1

if printed == 0:
    print("  (No candidates in the current window; widen MIN_SPEED..MAX_SPEED or refine ΔE/Δm.)")

# Save raw data
OUT_JSON = "/Users/steventippeconnic/Documents/QC/E_mc2_4_bit_Run_0.json"
payload = {
    "experiment": "E_equals_m_c2_mod16",
    "backend": backend.name,
    "physical_qubits": PHYSICAL,
    "shots": SHOTS,
    "order": ORDER,
    "N_Q": N_Q,
    "scales": {"DM_PER_TICK": DM_PER_TICK, "DE_PER_TICK": DE_PER_TICK, "C2_UNIT": C2_UNIT},
    "inputs": {"mass_kg": object_mass_kg, "energy_joule": object_energy_joule},
    "encoded": {"M_IDX": M_IDX, "E_IDX": E_IDX},
    "counts": counts_raw
}
with open(OUT_JSON, "w") as fp:
    json.dump(payload, fp, indent=2)
log.info("Results saved → %s", OUT_JSON)

# End


//////////////////////////////////////
# Code for all analysis from backend data


# Renders Circuit 0 (Full Backend Result Download)
# Imports
import json, math, io, base64, zlib
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# Load & Decode all shots
PATH = '/Users/steventippeconnic/Documents/QC/job-d3hviu9b641c738mlkmg-result.json'
with open(PATH, 'r') as fp:
    raw = json.load(fp)

def _find_bitarray_nodes(obj, out):
    """Recursively collect dicts that look like a BitArray holder with 'num_bits' and 'array'."""
    if isinstance(obj, dict):
        if 'num_bits' in obj and 'array' in obj:
            out.append(obj)
        for v in obj.values():
            _find_bitarray_nodes(v, out)
    elif isinstance(obj, list):
        for v in obj:
            _find_bitarray_nodes(v, out)
    return out

nodes = _find_bitarray_nodes(raw, [])
if not nodes:
    raise RuntimeError("Could not find BitArray data in the result JSON.")

# Concatenate all shots across any BitArray nodes.
all_uint8 = []
num_bits = None
for node in nodes:
    nb = int(node.get('num_bits', 0))
    if num_bits is None:
        num_bits = nb
    elif nb != num_bits:
        raise RuntimeError(f"Inconsistent num_bits across nodes: {num_bits} vs {nb}")

    arr_node = node['array']
    # Unwrap potential nesting like {'__type__':'ndarray','__value__':'...'}
    if isinstance(arr_node, dict) and '__value__' in arr_node:
        arr_payload = arr_node['__value__']
    else:
        arr_payload = arr_node

    if isinstance(arr_payload, str):
        by = base64.b64decode(arr_payload)
        try:
            buf = zlib.decompress(by)
        except Exception:
            buf = by
        # Try to load as .npy
        try:
            arr = np.load(io.BytesIO(buf), allow_pickle=False)
        except Exception:
            # Fallback: treat as raw uint8 bytes
            arr = np.frombuffer(buf, dtype=np.uint8)
    elif isinstance(arr_payload, list):
        arr = np.array(arr_payload)
    else:
        arr = np.array(arr_payload)

    # Accept common shapes
    if arr.ndim == 1 and np.issubdtype(arr.dtype, np.integer):
        all_uint8.append(arr.astype(np.uint8))
    elif arr.ndim == 2 and arr.shape[1] == num_bits:
        # Convert bits -> uint8 integers (MSB-first weight)
        weights = (1 << np.arange(num_bits-1, -1, -1, dtype=np.uint8))
        vals = (arr.astype(np.uint8) * weights[None, :]).sum(axis=1).astype(np.uint8)
        all_uint8.append(vals)
    else:
        # Try unpackbits if looks like packed bits
        arr_u8 = arr.view(np.uint8).ravel()
        bits = np.unpackbits(arr_u8, bitorder='big')
        if bits.size % num_bits != 0:
            raise RuntimeError("Cannot interpret array payload as shots.")
        shots = bits.reshape(-1, num_bits)
        weights = (1 << np.arange(num_bits-1, -1, -1, dtype=np.uint8))
        vals = (shots * weights[None, :]).sum(axis=1).astype(np.uint8)
        all_uint8.append(vals)

if not all_uint8:
    raise RuntimeError("No shots decoded.")
shots_uint8 = np.concatenate(all_uint8, axis=0)
N_SHOTS = shots_uint8.size

# Build (v,u) weights 
N_Q = num_bits // 2
ORDER = 1 << N_Q  # e.g., 16

def bit_reverse(x, nbits):
    y = 0
    for i in range(nbits):
        y = (y << 1) | ((x >> i) & 1)
    return y

# Split into halves (MSB..LSB), then reverse each half's bit order
u_vals = np.empty(N_SHOTS, dtype=np.int16)
v_vals = np.empty(N_SHOTS, dtype=np.int16)
mask_u = (1 << N_Q) - 1
for i, val in enumerate(shots_uint8):
    v_raw = (val >> N_Q) & mask_u   # upper 4 bits
    u_raw = val & mask_u            # lower 4 bits
    v_vals[i] = bit_reverse(v_raw, N_Q)
    u_vals[i] = bit_reverse(u_raw, N_Q)

# Aggregate to grid: grid_vu[v,u] = weight
grid_vu = np.zeros((ORDER, ORDER), dtype=float)
for u, v in zip(u_vals, v_vals):
    grid_vu[v, u] += 1.0

# 3D topology-only surface 
U, V = np.meshgrid(np.arange(ORDER), np.arange(ORDER), indexing='xy')  # x=v, y=u

fig = plt.figure(figsize=(11, 7.5))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(V, U, grid_vu.T, rstride=1, cstride=1, linewidth=0, antialiased=True,
                cmap='viridis', alpha=0.98)  # transpose: z = grid[u,v] under x=v,y=u

ax.set_xlabel('v')
ax.set_ylabel('u')
ax.set_zlabel('counts')
ax.set_title('All-shots interference landscape — topology only')
plt.tight_layout()
plt.show()

# End.

# Renders Circuit 1 (Full Backend Result Download)
# Imports
import os, json, base64, zlib, io
from collections import Counter
import numpy as np
import matplotlib.pyplot as plt

# Config 
# Default to user's path; allow override via EMC2_PATH for portability.
USER_PATH_DEFAULT = "/Users/steventippeconnic/Documents/QC/job-d3hviu9b641c738mlkmg-result.json"
PATH = os.environ.get("EMC2_PATH", USER_PATH_DEFAULT)

ORDER = 16
N_Q = 4
N_BITS = 8

# Helpers 
def decode_ndarray(ndarray_obj):
    """Decode a Qiskit JSON-encoded ndarray where '__value__' is base64 of zlib-compressed .npy bytes."""
    b64 = ndarray_obj["__value__"]
    raw = base64.b64decode(b64)
    decomp = zlib.decompress(raw)
    arr = np.load(io.BytesIO(decomp), allow_pickle=False)
    return arr

def ensure_bit_rows(arr, n_bits=N_BITS):
    a = np.array(arr)
    if a.ndim == 2 and a.shape[1] == n_bits:
        return (a != 0).astype(np.uint8)
    if a.ndim == 2 and a.dtype == np.uint8 and a.shape[1]*8 >= n_bits:
        unpack = np.unpackbits(a, axis=1)
        return unpack[:, :n_bits]
    if a.ndim == 1 and a.dtype == np.uint8:
        nbytes = (n_bits + 7)//8
        shots = a.size // nbytes
        a2 = a[:shots*nbytes].reshape(shots, nbytes)
        unpack = np.unpackbits(a2, axis=1)
        return unpack[:, :n_bits]
    if a.ndim == 2 and a.shape[1] % 8 == 0:
        unpack = np.unpackbits(a.astype(np.uint8), axis=1)
        return unpack[:, :n_bits]
    raise ValueError(f"Unrecognized bit array shape: {a.shape}, dtype={a.dtype}")

def bits_to_int_le(bs: str) -> int:
    """Little-endian bitstring to int (reverse string then base-2)."""
    return int(bs[::-1], 2)

def smooth2d(mat, sigma=0.9, ksize=5):
    """Separable Gaussian smoothing (no SciPy dependency)."""
    xs = np.linspace(-2, 2, ksize)
    kern1d = np.exp(-(xs**2)/(2*(sigma**2)))
    kern1d /= kern1d.sum()
    # Row convolution
    tmp = np.pad(mat, ((ksize//2, ksize//2), (0, 0)), mode='reflect')
    row_conv = np.zeros_like(mat, dtype=np.float64)
    for i in range(mat.shape[0]):
        sl = tmp[i:i+ksize, :]
        row_conv[i, :] = (sl * kern1d[:, None]).sum(axis=0)
    # Column convolution
    tmp2 = np.pad(row_conv, ((0, 0), (ksize//2, ksize//2)), mode='reflect')
    out = np.zeros_like(mat, dtype=np.float64)
    for j in range(mat.shape[1]):
        sl = tmp2[:, j:j+ksize]
        out[:, j] = (sl * kern1d[None, :]).sum(axis=1)
    return out

# Load & reconstruct 
with open(PATH, 'r') as f:
    raw = json.load(f)

pub0   = raw["__value__"]["pub_results"][0]["__value__"]
data   = pub0["data"]["__value__"]
fields = data["fields"]
bitarr = fields["c"]["__value__"]
assert bitarr["num_bits"] == N_BITS, f"Expected {N_BITS} bits, got {bitarr['num_bits']}"

arr      = decode_ndarray(bitarr["array"])
bitrows  = ensure_bit_rows(arr, N_BITS)
bits     = (bitrows > 0).astype(np.uint8)
bitstrs  = [''.join(str(b) for b in bits[i]) for i in range(bits.shape[0])]
counts   = Counter(bitstrs)

# Map raw "b a" bitstrings -> (a_val, b_val)
pairs = Counter()
for k, v in counts.items():
    a_bits = k[:N_Q]
    b_bits = k[N_Q:]
    a_val = bits_to_int_le(a_bits)
    b_val = bits_to_int_le(b_bits)
    pairs[(a_val, b_val)] += v

# Build lattice grid
grid = np.zeros((ORDER, ORDER), dtype=np.float64)
for (a, b), w in pairs.items():
    grid[a, b] += w

# 2D Phase Interference Lattice 
grid_sm = smooth2d(grid, sigma=0.9, ksize=5)

plt.figure(figsize=(7, 6), dpi=180)
im = plt.imshow(grid_sm.T, origin='lower', extent=[0, ORDER, 0, ORDER], aspect='equal')
plt.xlabel('a')
plt.ylabel('b')
plt.title('Phase Interference Lattice (smoothed)')
plt.contour(grid_sm.T, levels=8, linewidths=0.7)
plt.colorbar(im, fraction=0.046, pad=0.04)
plt.show()

# End.

# Renders Circuit 2 (Saved Json Backend Result)
# Imports
import json, math, random
from collections import Counter, defaultdict
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# Load results and helpers
PATH = '/Users/steventippeconnic/Documents/QC/E_mc2_4_bit_Run_8k_0.json'
with open(PATH, 'r') as fp:
    data = json.load(fp)

ORDER = int(data.get("order", 16))
N_Q   = int(data.get("N_Q", 4))
COUNTS_RAW = data["counts"]                  # bitstring -> frequency
M_IDX = int(data["encoded"]["M_IDX"])
E_IDX = int(data["encoded"]["E_IDX"])
C2_UNIT = float(data["scales"]["C2_UNIT"])   # ΔE/Δm (J/kg)
counts = data["counts"]

def bits_to_int(bs: str) -> int:
    # Reverse endian 
    return int(bs[::-1], 2)

# Recover (u,v) = (a,b) from raw bitstrings, same mapping: k ≡ -u * v^{-1} (mod N)
pairs = defaultdict(int)
grid_vu = np.zeros((ORDER, ORDER), dtype=float)  # rows=v, cols=u
for bitstr, cnt in COUNTS_RAW.items():
    a_val = bits_to_int(bitstr[N_Q:])   # right half = u (a)
    b_val = bits_to_int(bitstr[:N_Q])   # left  half = v (b)
    u = a_val % ORDER
    v = b_val % ORDER
    w = int(cnt)
    pairs[(u, v)] += w
    grid_vu[v, u] += w

# Reconstruct k’s per outcome where v=b is invertible
def inv_mod(x, N):
    # N is 2^n=16 invert only odd x (coprime to N)
    if math.gcd(x, N) != 1: 
        return None
    return pow(x, -1, N)

k_weights = Counter()
for (u, v), w in pairs.items():
    iv = inv_mod(v, ORDER)
    if iv is None: 
        continue
    k_res = (-u * iv) % ORDER
    k_weights[k_res] += w

# Theoretical target slope k_th = E_idx * M_idx^{-1} mod N
M_inv = inv_mod(M_IDX, ORDER)
k_th = (E_IDX * M_inv) % ORDER if M_inv is not None else None

# Recovered k per sample (skipping non-invertible v)
def recovered_k(u, v, N):
    iv = inv_mod(v, N)
    if iv is None:
        return None
    return (-u * iv) % N

# Expand weighted samples (8192-size typical, OK for bootstrap)
samples = []
for (u, v), w in pairs.items():
    samples.extend([(u, v)] * w)
samples = np.array(samples, dtype=np.int16)

# Histogram of recovered k residues 
ks = np.arange(ORDER)
ys = [k_weights.get(int(k), 0) for k in ks]

plt.figure(figsize=(8, 4))
plt.bar(ks, ys, width=0.8, align='center')
if k_th is not None:
    plt.axvline(k_th, color='r', linestyle='--', linewidth=2, label=f'theory k={k_th}')
    plt.legend()
plt.xticks(ks)
plt.xlabel('k residue  (E_idx * m_idx^{-1} mod N)')
plt.ylabel('weighted frequency')
plt.title('Recovered slope residues k from invertible rows')
plt.tight_layout()
plt.show()

# Polar residue wheel
theta = 2 * np.pi * ks / ORDER
r = np.array(ys, dtype=float)

plt.figure(figsize=(6, 6))
ax = plt.subplot(111, projection='polar')
ax.bar(theta, r, width=2*np.pi/ORDER, bottom=0.0, align='edge', edgecolor='k', linewidth=0.5)
if k_th is not None:
    ax.plot([2*np.pi*k_th/ORDER, 2*np.pi*k_th/ORDER], 
	[0, r.max()*1.05], 'r--', linewidth=2, label=f'theory k={k_th}')
    ax.legend(loc='upper left', bbox_to_anchor=(1.05, 1.05))
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.set_title('Residue wheel: weight by k on Z_N', va='bottom')
plt.tight_layout()
plt.show()

# Lifted c estimates from dominant residues 
V_MIN, V_MAX = 1.0e8, 4.0e8
def lifted_c_candidates_for_k(k, C2, N, vmin, vmax):
    c2_min, c2_max = vmin*vmin, vmax*vmax
    # smallest non-negative t that places c into window, else nearest plausible t
    t_low  = math.ceil((c2_min / C2 - k) / N)
    t_high = math.floor((c2_max / C2 - k) / N)
    if t_low <= t_high:
        # pick central t in allowable band and also show edges (up to 3 candidates)
        picks = [t_low, (t_low + t_high)//2, t_high]
    else:
        # nothing lands in window pick the nearest t to c ≈ sqrt(k*C2)
        t_guess = max(0, round(((vmin+vmax)/2)**2 / C2 / N))
        picks = [t_guess]
    out = []
    seen = set()
    for t in picks:
        if t in seen:
            continue
        seen.add(t)
        c2 = (k + t*N) * C2
        if c2 <= 0:
            continue
        c = math.sqrt(c2)
        out.append((t, c))
    return out[:3]

# Take top residues by weight (up to 8) and compute one best candidate each (closest to mid-window)
mid = 0.5*(V_MIN + V_MAX)
top_k = [kv[0] for kv in k_weights.most_common(8)]   # ensure we use the same weights source
cand_list = []
for k in top_k:
    cands = lifted_c_candidates_for_k(k, C2_UNIT, ORDER, V_MIN, V_MAX)
    if not cands:
        continue
    # prefer candidate closest to mid-window
    t_best, c_best = sorted(cands, key=lambda tc: abs(tc[1]-mid))[0]
    cand_list.append((k, int(k_weights[k]), t_best, c_best))  # show the exact weight for this k

# Sort by weight then by closeness to mid-window
cand_list.sort(key=lambda row: (-row[1], abs(row[3]-mid)))

# Horizontal bar c (m/s) vs residue k
if cand_list:
    labels = [f'k={k} (t={t})' for (k, w, t, c) in cand_list]
    cs = [c for (_, _, _, c) in cand_list]
    ws = [w for (_, w, _, _) in cand_list]

    plt.figure(figsize=(8, 5))
    y = np.arange(len(cs))
    plt.barh(y, cs)
    for i, (lbl, cval, w) in enumerate(zip(labels, cs, ws)):
        plt.text(
            cval + 0.05e8,  # small offset to the right
            i,
            f'{lbl}\nc≈{cval:0.3e} m/s | weight={w}',
            va='center',
            ha='left',
            fontsize=9,
            bbox=dict(facecolor='white', edgecolor='black', boxstyle='round,pad=0.5')
        )
    if k_th is not None:
        # Overlay theoretical c for k_th 
        t_th = 0
        c_th = math.sqrt((k_th + t_th*ORDER)*C2_UNIT) if (k_th + t_th*ORDER) > 0 else None
        if c_th is None or not (V_MIN <= c_th <= V_MAX):
            t_th = 1
            c_th = math.sqrt((k_th + t_th*ORDER)*C2_UNIT)
        plt.axvline(c_th, color='r', linestyle='--', linewidth=2, label=f'theory c for k={k_th}, t={t_th}')
        plt.legend(loc='lower right')
    plt.xlabel('c (m/s)')
    plt.yticks(y, [f'#{i+1}' for i in range(len(cs))])
    plt.title('Lifted c estimates from dominant residues (ΔE/Δm = {:0.3e} J/kg)'.format(C2_UNIT))
    plt.xlim(V_MIN*0.9, V_MAX*1.05)
    plt.tight_layout()
    plt.show()
else:
    plt.figure()
    plt.text(0.5, 0.5, 'No lifted c candidates in window.\nAdjust V_MIN/V_MAX or scaling.',
             ha='center', va='center')
    plt.axis('off')
    plt.show()

# Residuals
def centered_residual(u, v, k, N):
    r = (u + (k * v)) % N
    return ((r + N//2) % N) - N//2  # in [-N/2, N/2)

invertible_vs = [v for v in range(ORDER) if math.gcd(v, ORDER) == 1]

# Ridge Phase Ribbon per-v mean phase & coherence
angles = []
magnitudes = []
for v in invertible_vs:
    num = 0+0j
    den = 0.0
    for (u, vv), w in pairs.items():
        if vv != v:
            continue
        theta = 2*np.pi * ((u + (k_th * v)) % ORDER) / ORDER
        num += w * np.exp(1j*theta)
        den += w
    z = num / max(den, 1.0)
    angles.append(np.angle(z))       # in [-pi, pi]
    magnitudes.append(np.abs(z))     # resultant length in [0,1]

angles = np.unwrap(np.array(angles))  # smooth the phase along v
magnitudes = np.array(magnitudes)

plt.figure(figsize=(10, 4))
plt.plot(invertible_vs, angles, lw=2, label='mean phase after k_th unwrapping')
upper = angles + (1 - magnitudes) * 0.25
lower = angles - (1 - magnitudes) * 0.25
plt.fill_between(invertible_vs, lower, upper, alpha=0.25, step='mid', label='coherence band (thinner = sharper)')
plt.axhline(0.0, ls='--', lw=1, color='k', alpha=0.6)
plt.yticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi], [r'$-\pi$', r'$-\pi/2$', '0', r'$\pi/2$', r'$\pi$'])
plt.xlabel('v (invertible only)')
plt.ylabel('phase (rad)')
plt.title('Ridge Phase Ribbon (per-v phase & coherence after k_th unwrapping)')
plt.legend()
plt.tight_layout()
plt.show()

# Oblique Ridge Tomography: stacked profiles along r0
r0_list = [ -2, -1, 0, 1, 2 ]  # small band around the main ridge
profiles = {r0: np.zeros(ORDER, dtype=float) for r0 in r0_list}  # profile over v
mass_v = Counter()
for (u, v), w in pairs.items():
    mass_v[v] += w
    for r0 in r0_list:
        if centered_residual(u, v, k_th, ORDER) == r0:
            profiles[r0][v] += w

# normalize each profile by column mass to remove sampling bias
for r0 in r0_list:
    for v in range(ORDER):
        m = mass_v[v] if mass_v[v] > 0 else 1.0
        profiles[r0][v] = profiles[r0][v] / m
plt.figure(figsize=(10, 5))
offset = 0.0
delta  = 0.18
for r0 in r0_list:
    # small vertical offset so bands don't overlap
    y = profiles[r0].copy()
    plt.plot(range(ORDER), y + offset, label=f'r0={r0}')
    offset += delta

plt.xticks(range(ORDER))
plt.xlabel('v (all columns)')
plt.ylabel('normalized line-integral (stacked)')
plt.title('Oblique Ridge Tomography (parallel lines near the ridge)')
plt.legend(ncol=len(r0_list), fontsize=9)
plt.tight_layout()
plt.show()

# Bootstrap confidence for the modal k (with 95% CI)
B = 400                   
rng = np.random.default_rng(7)
modes = []

def mode_of_k_from_sample(arr, N):
    # compute k for invertible v only; return mode (ties broken by weight)
    ks = []
    for (u, v) in arr:
        k = recovered_k(int(u), int(v), N)
        if k is not None:
            ks.append(k)
    if not ks:
        return None
    counts = Counter(ks)
    
    max_w = max(counts.values())
    cands = [k for k, w in counts.items() if w == max_w]
    return int(min(cands))

for _ in range(B):
    idx = rng.integers(0, len(samples), size=len(samples))
    boot = samples[idx]
    m = mode_of_k_from_sample(boot, ORDER)
    modes.append(m)

# summarize bootstrap modal probabilities
mode_counts = Counter([m for m in modes if m is not None])
ks_sorted = sorted(mode_counts.keys())
probs = np.array([mode_counts[k]/B for k in ks_sorted])

errs = np.sqrt(probs*(1-probs)/B)

plt.figure(figsize=(8, 4))
plt.bar(ks_sorted, probs, yerr=errs, capsize=3)
if k_th is not None:
    plt.axvline(k_th, linestyle='--', linewidth=2, label=f'theory k={k_th}')
    plt.legend()
plt.xticks(range(ORDER))
plt.ylim(0, 1)
plt.xlabel('k residue (modal outcome)')
plt.ylabel('bootstrap probability  P(mode = k)')
plt.title('Bootstrap confidence for the modal slope k (B = {})'.format(B))
plt.tight_layout()
plt.show()

# End

# Renders Circuit 3 (Full Backend Result Download)
import os, json, base64, zlib, io, math
from collections import Counter, defaultdict
import numpy as np
import matplotlib.pyplot as plt

# Config 
USER_PATH_DEFAULT = "/Users/steventippeconnic/Documents/QC/job-d3hviu9b641c738mlkmg-result.json"
PATH = os.environ.get("EMC2_PATH", USER_PATH_DEFAULT)

ORDER = 16
N_Q = 4
N_BITS = 8
DM_PER_TICK = 1.0
DE_PER_TICK = 1.0e16
C2_UNIT = DE_PER_TICK / DM_PER_TICK
MIN_SPEED = 1.0e8
MAX_SPEED = 4.0e8

# Helpers 
def decode_ndarray(ndarray_obj):
    b64 = ndarray_obj["__value__"]
    raw = base64.b64decode(b64)
    decomp = zlib.decompress(raw)
    arr = np.load(io.BytesIO(decomp), allow_pickle=False)
    return arr

def ensure_bit_rows(arr, n_bits=N_BITS):
    a = np.array(arr)
    if a.ndim == 2 and a.shape[1] == n_bits:
        return (a != 0).astype(np.uint8)
    if a.ndim == 2 and a.dtype == np.uint8 and a.shape[1]*8 >= n_bits:
        unpack = np.unpackbits(a, axis=1)
        return unpack[:, :n_bits]
    if a.ndim == 1 and a.dtype == np.uint8:
        nbytes = (n_bits + 7)//8
        shots = a.size // nbytes
        a2 = a[:shots*nbytes].reshape(shots, nbytes)
        unpack = np.unpackbits(a2, axis=1)
        return unpack[:, :n_bits]
    if a.ndim == 2 and a.shape[1] % 8 == 0:
        unpack = np.unpackbits(a.astype(np.uint8), axis=1)
        return unpack[:, :n_bits]
    raise ValueError(f"Unrecognized bit array shape: {a.shape}, dtype={a.dtype}")

def bits_to_int_le(bs):
    return int(bs[::-1], 2)

def k_from_pair(a, b, mod=ORDER):
    if math.gcd(b, mod) != 1:
        return None
    inv_b = pow(b, -1, mod)
    return (-a * inv_b) % mod

def c_candidates_from_k(k, c2_unit=C2_UNIT, vmin=MIN_SPEED, vmax=MAX_SPEED, order=ORDER):
    c2_min = vmin*vmin
    c2_max = vmax*vmax
    t_low  = math.ceil((c2_min / c2_unit - k) / order)
    t_high = math.floor((c2_max / c2_unit - k) / order)
    out = []
    for t in range(t_low, t_high+1):
        val = (k + t*order) * c2_unit
        if val <= 0:
            continue
        c = math.sqrt(val)
        if vmin <= c <= vmax:
            out.append((t, c))
    return out

def gaussian(x, mu, sigma):
    return np.exp(-0.5*((x-mu)/sigma)**2)

# Load backend data 
with open(PATH, 'r') as f:
    raw = json.load(f)

pub0   = raw["__value__"]["pub_results"][0]["__value__"]
data   = pub0["data"]["__value__"]
fields = data["fields"]
bitarr = fields["c"]["__value__"]
assert bitarr["num_bits"] == N_BITS, "Unexpected classical width"

arr      = decode_ndarray(bitarr["array"])
bitrows  = ensure_bit_rows(arr, N_BITS)
bits     = (bitrows > 0).astype(np.uint8)

bitstrings = [''.join(str(b) for b in bits[i]) for i in range(bits.shape[0])]
counts     = Counter(bitstrings)

# Compute k-weights 
pairs = Counter()
for s, v in counts.items():
    a_bits = s[:N_Q]
    b_bits = s[N_Q:]
    a_val  = bits_to_int_le(a_bits)
    b_val  = bits_to_int_le(b_bits)
    pairs[(a_val, b_val)] += v

k_weight_backend = Counter()
for (a, b), w in pairs.items():
    k = k_from_pair(a, b, ORDER)
    if k is not None:
        k_weight_backend[k] += w

# Build continuation density 
k_all = list(range(ORDER))
c_grid = np.linspace(MIN_SPEED, MAX_SPEED, 1600)
sigma_c = 0.012*(MAX_SPEED - MIN_SPEED)
density_per_k = defaultdict(lambda: np.zeros_like(c_grid))

for k in k_all:
    w = k_weight_backend.get(k, 0)
    if w <= 0:
        continue
    for (_, c) in c_candidates_from_k(k):
        density_per_k[k] += gaussian(c_grid, c, sigma_c) * w
    m = density_per_k[k].max()
    if m > 0:
        density_per_k[k] /= m  # per-k normalize for contrast

# Visual: k–c continuation map 
kc_img = np.zeros((ORDER, c_grid.size), dtype=float)
for k in range(ORDER):
    kc_img[k] = density_per_k.get(k, np.zeros_like(c_grid))

plt.figure(figsize=(9, 6))
plt.imshow(kc_img, origin='lower', aspect='auto',
           extent=[c_grid[0], c_grid[-1], 0, ORDER])
plt.xlabel('c (m/s)')
plt.ylabel('k (mod 16)')
plt.title('k–c Continuation Map (Per-k Spectral Density)')
plt.colorbar(label='normalized intensity')
plt.show()

# Renders Circuit 4 (Backend Json)
import json, math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  
from collections import defaultdict

# Load results + helpers
PATH = '/Users/steventippeconnic/Documents/QC/E_mc2_4_bit_Run_8k_0.json'
with open(PATH, 'r') as fp:
    data = json.load(fp)

ORDER  = int(data.get("order", 16))
N_Q    = int(data.get("N_Q", 4))
counts = data["counts"]
M_IDX  = int(data["encoded"]["M_IDX"])
E_IDX  = int(data["encoded"]["E_IDX"])

def bits_to_int(bs: str) -> int:
    # reverse-endian to match your earlier convention
    return int(bs[::-1], 2)

def inv_mod(x, N):
    return pow(x, -1, N) if math.gcd(x, N) == 1 else None

# Build weighted (u,v) grid
grid = np.zeros((ORDER, ORDER), dtype=float)  # rows=u (a), cols=v (b)
pairs = defaultdict(int)
for bitstr, cnt in counts.items():
    u = bits_to_int(bitstr[N_Q:]) % ORDER  # measured a
    v = bits_to_int(bitstr[:N_Q])  % ORDER  # measured b
    grid[u, v] += cnt
    pairs[(u, v)] += int(cnt)

# Theory slope k_th = E_idx * M_idx^{-1} mod N
M_inv = inv_mod(M_IDX, ORDER)
k_th = (E_IDX * M_inv) % ORDER if M_inv is not None else None

# Small helpers
def ridge_u_of_v(k, v, N):
    return (-k * v) % N

def take_along_ridge(grid, k, N):
    # For each v, pick u=-k v (mod N) and return heights
    vs = np.arange(N)
    us = ridge_u_of_v(k, vs, N)
    h  = grid[us, vs]
    return vs, us, h

# 3D surface of counts over (v,u) with the predicted ridge as an elevated ribbon
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

V, U = np.meshgrid(np.arange(ORDER), np.arange(ORDER))
ax.plot_surface(V, U, grid, rstride=1, cstride=1, linewidth=0, antialiased=True)

# Overlay predicted ridge as a lifted ribbon (height = counts along line)
if k_th is not None:
    vs, us, hs = take_along_ridge(grid, k_th, ORDER)
    ax.plot(vs, us, hs, linewidth=3)

ax.set_xlabel('v  (measured b)')
ax.set_ylabel('u  (measured a)')
ax.set_zlabel('counts')
ax.set_title('3D surface: interference landscape with ridge ribbon')
plt.tight_layout()
plt.show()

# Ridge tube vs anti-ridge: 3D line profiles
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

# Plot the entire surface faintly to give context (wire)
ax.plot_wireframe(V, U, grid, rstride=2, cstride=2, linewidth=0.5)

if k_th is not None:
    vs, us, hs = take_along_ridge(grid, k_th, ORDER)
    # anti-ridge: half-cycle away
    us_anti = (us + ORDER//2) % ORDER
    hs_anti = grid[us_anti, vs]

    ax.plot(vs, us, hs, linewidth=3)         # ridge tube
    ax.plot(vs, us_anti, hs_anti, linewidth=3)  # anti-ridge tube

ax.set_xlabel('v')
ax.set_ylabel('u')
ax.set_zlabel('counts')
ax.set_title('3D ridge tube vs anti-ridge tube')
plt.tight_layout()
plt.show()


# Torus embedding of Z_N × Z_N: points on a torus with ridge geodesic overlay
R_major, r_minor = 3.0, 1.0
phi = 2*np.pi*np.arange(ORDER)/ORDER    # for u
psi = 2*np.pi*np.arange(ORDER)/ORDER    # for v

# Point cloud on torus
Xs, Ys, Zs, Ss = [], [], [], []
for u in range(ORDER):
    for v in range(ORDER):
        w = grid[u, v]
        if w <= 0:
            continue
        Φ = phi[u]
        Ψ = psi[v]
        X = (R_major + r_minor*np.cos(Φ)) * np.cos(Ψ)
        Y = (R_major + r_minor*np.cos(Φ)) * np.sin(Ψ)
        Z =  r_minor*np.sin(Φ)
        Xs.append(X); Ys.append(Y); Zs.append(Z); Ss.append(w)

Ss = np.array(Ss, dtype=float)
Ss = 10.0 * (Ss / Ss.max())**0.5 * 20  # gentle normalization for marker size

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(Xs, Ys, Zs, s=Ss)

# Overlay predicted ridge geodesic: φ = -k_th ψ (mod 2π)
if k_th is not None:
    t = np.linspace(0, 2*np.pi, 1000)
    Ψ = t
    Φ = (-k_th * t) % (2*np.pi)
    X = (R_major + r_minor*np.cos(Φ)) * np.cos(Ψ)
    Y = (R_major + r_minor*np.cos(Φ)) * np.sin(Ψ)
    Z =  r_minor*np.sin(Φ)
    ax.plot(X, Y, Z, linewidth=3)

ax.set_title('Torus embedding: interference ridge as a geodesic')
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
plt.tight_layout()
plt.show()

# End

# Renders Circuit 5 (Full Backend Result Download)
# Imports
import os, json, base64, zlib, io
from collections import Counter
import numpy as np
import matplotlib.pyplot as plt

# Config 
USER_PATH_DEFAULT = "/Users/steventippeconnic/Documents/QC/job-d3hviu9b641c738mlkmg-result.json"
PATH = os.environ.get("EMC2_PATH", USER_PATH_DEFAULT)

ORDER = 16
N_Q = 4
N_BITS = 8

# Helpers 
def decode_ndarray(ndarray_obj):
    """Decode a Qiskit JSON-encoded ndarray where '__value__' is base64 of zlib-compressed .npy bytes."""
    b64 = ndarray_obj["__value__"]
    raw = base64.b64decode(b64)
    decomp = zlib.decompress(raw)
    arr = np.load(io.BytesIO(decomp), allow_pickle=False)
    return arr

def ensure_bit_rows(arr, n_bits=N_BITS):
    a = np.array(arr)
    if a.ndim == 2 and a.shape[1] == n_bits:
        return (a != 0).astype(np.uint8)
    if a.ndim == 2 and a.dtype == np.uint8 and a.shape[1]*8 >= n_bits:
        unpack = np.unpackbits(a, axis=1)
        return unpack[:, :n_bits]
    if a.ndim == 1 and a.dtype == np.uint8:
        nbytes = (n_bits + 7)//8
        shots = a.size // nbytes
        a2 = a[:shots*nbytes].reshape(shots, nbytes)
        unpack = np.unpackbits(a2, axis=1)
        return unpack[:, :n_bits]
    if a.ndim == 2 and a.shape[1] % 8 == 0:
        unpack = np.unpackbits(a.astype(np.uint8), axis=1)
        return unpack[:, :n_bits]
    raise ValueError(f"Unrecognized bit array shape: {a.shape}, dtype={a.dtype}")

def bits_to_int_le(bs: str) -> int:
    """Little-endian bitstring to int (reverse string then base-2)."""
    return int(bs[::-1], 2)

def smooth2d(mat, sigma=0.9, ksize=5):
    """Separable Gaussian smoothing (no SciPy dependency)."""
    xs = np.linspace(-2, 2, ksize)
    kern1d = np.exp(-(xs**2)/(2*(sigma**2)))
    kern1d /= kern1d.sum()
    # Row convolution
    tmp = np.pad(mat, ((ksize//2, ksize//2), (0, 0)), mode='reflect')
    row_conv = np.zeros_like(mat, dtype=np.float64)
    for i in range(mat.shape[0]):
        sl = tmp[i:i+ksize, :]
        row_conv[i, :] = (sl * kern1d[:, None]).sum(axis=0)
    # Column convolution
    tmp2 = np.pad(row_conv, ((0, 0), (ksize//2, ksize//2)), mode='reflect')
    out = np.zeros_like(mat, dtype=np.float64)
    for j in range(mat.shape[1]):
        sl = tmp2[:, j:j+ksize]
        out[:, j] = (sl * kern1d[None, :]).sum(axis=1)
    return out

# Load & reconstruct lattice 
with open(PATH, 'r') as f:
    raw = json.load(f)

pub0   = raw["__value__"]["pub_results"][0]["__value__"]
data   = pub0["data"]["__value__"]
fields = data["fields"]
bitarr = fields["c"]["__value__"]
assert bitarr["num_bits"] == N_BITS, f"Expected {N_BITS} bits, got {bitarr['num_bits']}"

arr      = decode_ndarray(bitarr["array"])
bitrows  = ensure_bit_rows(arr, N_BITS)
bits     = (bitrows > 0).astype(np.uint8)
bitstrs  = [''.join(str(b) for b in bits[i]) for i in range(bits.shape[0])]
counts   = Counter(bitstrs)

# Map raw "b a" bitstrings -> (a_val, b_val)
pairs = Counter()
for s, v in counts.items():
    a_bits = s[:N_Q]
    b_bits = s[N_Q:]
    a_val = bits_to_int_le(a_bits)
    b_val = bits_to_int_le(b_bits)
    pairs[(a_val, b_val)] += v

grid = np.zeros((ORDER, ORDER), dtype=np.float64)
for (a, b), w in pairs.items():
    grid[a, b] += w

# Visual: Phase-Slope Vector Field 
grid_sm = smooth2d(grid, sigma=0.9, ksize=5)
Gy, Gx = np.gradient(grid_sm)  # rows=a (y), cols=b (x)

# Downsample vectors to avoid clutter
step = 2
A_coords, B_coords = np.mgrid[0:ORDER:step, 0:ORDER:step]
Ux = Gx[::step, ::step]
Uy = Gy[::step, ::step]

plt.figure(figsize=(7, 6), dpi=180)
plt.quiver(B_coords, A_coords, Ux, Uy, angles='xy', scale_units='xy', scale=1)
plt.xlim(0, ORDER-1); plt.ylim(0, ORDER-1)
plt.xlabel('b')
plt.ylabel('a')
plt.title('Phase-Slope Vector Field (gradient of smoothed lattice)')
plt.show()

# End