Skip to content

Added cultivation simulation code #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
stevenhofmeyr merged 1 commit into from
Jul 10, 2026
Merged

Added cultivation simulation code #3

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
237 changes: 237 additions & 0 deletions src/mcmc_cultivation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,237 @@
#!/usr/bin/env python3

"""
MCMC simulation of magic state cultivation time distribution.

Survival fractions are read from Figure 15 of Gidney et al.
(arXiv:2409.17595). Each array entry is the proportion of shots surviving
to that round (i.e. the y-value of the survival curve at that x position).
Values must be non-increasing and start at 1.0.

The x-axis of Figure 15 has:
- d1=3 case: 18 labelled cycles
- d1=5 case: 25 labelled cycles

The final ~10 cycles are the escape stage gap-wait, during which survival
holds flat (no postselection). The final survival rates from Figure 2 are:
- d1=3: ~20% (discard rate 80%)
- d1=5: ~1% (discard rate 99%)
"""

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

# ---------------------------------------------------------------------------
# Cumulative survival fractions from Figure 15 of Gidney et al.
# Each entry is the proportion of shots surviving TO that round
# (the y-value of the survival curve at that x position).
# Values must be non-increasing; final value = end-to-end survival rate.
#
# Structure for d1=5 (25 cycles total):
# - ~3 injection cycles (survival drops sharply)
# - ~12 cultivation cycles (check-grow-stabilize, survival drops each)
# - ~10 escape/gap-wait cycles (survival flat - no postselection)
#
# Final values from Figure 2:
# d1=3: ~0.20 (discard rate 80%)
# d1=5: ~0.01 (discard rate 99%)
# ---------------------------------------------------------------------------

# d1=3 case: 18 cycles, final survival ~20%
# Rough structure: 3 injection + 5 cultivation + 10 gap-wait
# We distribute the 80% total loss across the first 8 cycles.
# Adjust these values to match the actual curve shape in Figure 15.
per_cycle_survival_d3 = [
1.0, # encode T
1.0, # stabilize
0.83, # check T
0.83, # check T
0.74, # stabilize
0.64, # stabilize
0.6, # stabilize
0.57, # escaped
0.55, # wait
0.55, # wait
0.55, # wait
0.55, # wait
0.55, # wait
0.55, # wait
0.55, # wait
0.55, # wait
0.55, # wait
0.25 # ready
]
assert len(per_cycle_survival_d3) == 18, f"Expected 18 cycles, got {len(per_cycle_survival_d3)}"

# d1=5 case: 25 cycles, final survival ~1%
# Rough structure: 3 injection + 12 cultivation + 10 gap-wait
# We distribute the 99% total loss across the first 15 cycles.
# Adjust these values to match the actual curve shape in Figure 15.
per_cycle_survival_d5 = [
1.0, # encode T
1.0, # stabilize
0.83, # check T
0.83, # check T
0.72, # stabilze
0.53, # stabilize
0.4, # stabilize
0.3, # check T
0.3, # check T
0.19, # stabilize
0.12, # stabilize
0.1, # stabilize
0.08, # stabilize
0.075, # stabilize
0.07, # escaped
0.06, # wait
0.06, # wait
0.06, # wait
0.06, # wait
0.06, # wait
0.06, # wait
0.06, # wait
0.06, # wait
0.06, # wait
0.018, # ready
]
assert len(per_cycle_survival_d5) == 25, f"Expected 25 cycles, got {len(per_cycle_survival_d5)}"

# Convert cumulative survival fractions to conditional per-cycle probabilities
# cond[i] = P(survive cycle i | survived all previous cycles)
# = cumulative[i] / cumulative[i-1]
def cumulative_to_conditional(cumulative):
arr = np.array(cumulative, dtype=float)
cond = np.empty_like(arr)
cond[0] = arr[0]
for i in range(1, len(arr)):
cond[i] = arr[i] / arr[i - 1] if arr[i - 1] > 0 else 0.0
return cond

# Report end-to-end survival (final value of each array)
cum_d3 = per_cycle_survival_d3[-1]
cum_d5 = per_cycle_survival_d5[-1]
print(f"d1=3 end-to-end survival: {cum_d3:.4f} (target ~0.20)")
print(f"d1=5 end-to-end survival: {cum_d5:.4f} (target ~0.01)")

cond_d3 = cumulative_to_conditional(per_cycle_survival_d3)
cond_d5 = cumulative_to_conditional(per_cycle_survival_d5)


def run_mcmc(conditional_probs, n_samples=1_000_000):
"""
Run MCMC simulation of cultivation time.

conditional_probs[i] = probability of surviving cycle i given all
previous cycles survived (derived from cumulative survival fractions).

Each attempt proceeds cycle by cycle. If it fails at any cycle, the
attempt is discarded and a new one starts from cycle 0. If it survives
all cycles, the magic state is successfully produced.

Returns an array of total cycle counts until success (including all
failed attempts).
"""
probs = np.array(conditional_probs)
success_times = []
report_every = max(1, n_samples // 20)

for i in range(n_samples):
if i % report_every == 0:
pct = 100.0 * i / n_samples
print(f" progress: {i:>{len(str(n_samples))}}/{n_samples} ({pct:5.1f}%)", end="\r", flush=True)
total_cycles = 0
while True:
survived = True
for p in probs:
total_cycles += 1
if rng.random() > p:
survived = False
break
if survived:
success_times.append(total_cycles)
break

print(f" progress: {n_samples}/{n_samples} (100.0%)", flush=True)
return np.array(success_times)


def fit_exponential(times):
"""Fit an exponential distribution to the cultivation times."""
# MLE for exponential: lambda = 1/mean
mean_time = np.mean(times)
lam = 1.0 / mean_time
return lam, mean_time


def plot_results(times_d3, times_d5, lam_d3, lam_d5):
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for ax, times, lam, label, color in zip(
axes,
[times_d3, times_d5],
[lam_d3, lam_d5],
["d1=3", "d1=5"],
["steelblue", "darkorange"],
):
mean_t = 1.0 / lam
# CDF
sorted_t = np.sort(times)
cdf = np.arange(1, len(sorted_t) + 1) / len(sorted_t)
ax.plot(sorted_t, cdf, color=color, lw=1.5, label="Empirical CDF")

# Fitted exponential CDF
t_range = np.linspace(0, sorted_t[-1], 1000)
ax.plot(t_range, 1 - np.exp(-lam * t_range), "k--", lw=1.5,
label=f"Exp fit: λ={lam:.5f}\nmean={mean_t:.1f} cycles")

ax.set_xlabel("Circuit cycles until success")
ax.set_ylabel("CDF")
ax.set_title(f"Cultivation time distribution ({label})")
ax.legend()
ax.grid(True, alpha=0.4)

plt.tight_layout()
plt.savefig("cultivation_time_distribution.png", dpi=150)
plt.show()
print("Saved cultivation_time_distribution.png")


if __name__ == "__main__":
N_SAMPLES = 100_000

n_cycles_d3 = len(per_cycle_survival_d3)
n_cycles_d5 = len(per_cycle_survival_d5)

print(f"\nRunning MCMC for d1=3 ({N_SAMPLES:,} samples)...")
times_d3 = run_mcmc(cond_d3, n_samples=N_SAMPLES)
lam_d3, mean_d3 = fit_exponential(times_d3)
print(f" Mean cultivation time: {mean_d3:.1f} circuit cycles")
print(f" Fitted λ: {lam_d3:.6f} per circuit cycle")
print(f" Median: {np.median(times_d3):.1f} cycles")
print(f" 90th percentile: {np.percentile(times_d3, 90):.1f} cycles")

print(f"\nRunning MCMC for d1=5 ({N_SAMPLES:,} samples)...")
times_d5 = run_mcmc(cond_d5, n_samples=N_SAMPLES)
lam_d5, mean_d5 = fit_exponential(times_d5)
print(f" Mean cultivation time: {mean_d5:.1f} circuit cycles")
print(f" Fitted λ: {lam_d5:.6f} per circuit cycle")
print(f" Median: {np.median(times_d5):.1f} cycles")
print(f" 90th percentile: {np.percentile(times_d5, 90):.1f} cycles")

# Compare to paper's claimed value
paper_lam = 0.00227
paper_mean_cycles = 1.0 / paper_lam
paper_d = 17
paper_mean_logical = paper_mean_cycles / paper_d
print(f"\nPaper's claimed λ={paper_lam}, mean={paper_mean_cycles:.1f} circuit cycles")
print(f" = {paper_mean_logical:.1f} logical cycles (at d={paper_d})")

print(f"\nAnalytic estimate (N_cycles / end-to-end survival):")
print(f" d1=3: {n_cycles_d3} / {cum_d3:.4f} = {n_cycles_d3 / cum_d3:.1f} circuit cycles"
f" = {n_cycles_d3 / cum_d3 / paper_d:.1f} logical cycles")
print(f" d1=5: {n_cycles_d5} / {cum_d5:.4f} = {n_cycles_d5 / cum_d5:.1f} circuit cycles"
f" = {n_cycles_d5 / cum_d5 / paper_d:.1f} logical cycles")

plot_results(times_d3, times_d5, lam_d3, lam_d5)
Loading