diff --git a/benchmarking/gari/README.md b/benchmarking/gari/README.md new file mode 100644 index 0000000..c1e0ede --- /dev/null +++ b/benchmarking/gari/README.md @@ -0,0 +1,62 @@ +# GARI Benchmarks + + +## Contents + +- `submit.sh`: Slurm submission script for the benchmark sweep. +- `submit_locally.sh`: Local bash script execution equivalent for the benchmark sweep. +- `aggregated_results.jsonl`: Aggregated benchmark results (not checked into source control). +- `interactive_plot.py`: Interactive html dashboard for plotting and visualizing performance tradeoffs dynamically in the browser. + +## steps before running jobs + +```bash +bazel build src:tesseract +``` + +Generate all the gari dems with all prior modes and detector ordering for the targeted circuits. +From the repository root: +# Run DEM Generation under Bazel: +``` +bazel run //src/py/_tesseract_py_util:gari_dem_utils -- "testdata/bivariatebicyclecodes/" +bazel run //src/py/_tesseract_py_util:gari_dem_utils -- "testdata/colorcodes/" +bazel run //src/py/_tesseract_py_util:gari_dem_utils -- "testdata/surfacecodes/" +``` + +# Run test simulation under Bazel: +``` +bazel run //src/py:gari_simulation_test +``` + +The gari_dem_utils and simulation scripts require `numpy`, `scipy`, `matplotlib`, and `stim`, which are managed by Bazel when run via `bazel run`. + +## Re-running jobs + +From the repository root: + +```bash +benchmarking/gari/submit.sh +``` + +The script assumes the Tesseract binary is available at +`./bazel-bin/src/tesseract`, reads circuits from `testdata/`, gari dems for the corresponding circuits, and submits jobs +with `sbatch`. The Slurm partition, memory, CPU count, and walltime are tuned +for the cluster used for the original PR benchmark and may need adjustment +before reuse. + +Per-job stats are written under `out/`. Those raw per-job JSON files are not +included here; `aggregated_results.jsonl` is the aggregated dataset used for the +published plots. + +## Plotting and Visualization + +From this directory: + +```bash +python3 interactive_plot.py +``` + +The `interactive_plot.py` script reads the local `aggregated_results.jsonl` data and outputs an interactive HTML dashboard (`interactive_plot.html`) that can be opened in any web browser. +It requires `bokeh`, `pandas`, `numpy`, and `jupyter_bokeh`. + +Within the interactive dashboard, you can filter against Code Families, Distances, and other hyperparmeters, dynamically generating visual speedup baselines and Logical Error Rate plots. diff --git a/benchmarking/gari/interactive_plot.py b/benchmarking/gari/interactive_plot.py new file mode 100644 index 0000000..5864258 --- /dev/null +++ b/benchmarking/gari/interactive_plot.py @@ -0,0 +1,877 @@ +import os +import math +import json +from plot import process_data, compute_metrics, format_gari_label +from bokeh.plotting import figure, save, output_file +from bokeh.models import ColumnDataSource, CustomJS, CustomJSFilter, CDSView, HoverTool, TapTool, CheckboxGroup, Select, Button, Div, LabelSet +from bokeh.layouts import column, row + +def generate_interactive_plot(): + # load metrics + data_groups = process_data('aggregated_results.jsonl') + metrics = compute_metrics(data_groups) + + # compute tradeoff data for all points + baselines = {} + for m in metrics: + if m['is_baseline'] and m.get('num_det_orders', 1) == 1 and m.get('decoder', 'tesseract') == 'tesseract': + key = (m['p'], m['type'], m['d'], m['q'], m['r']) + baselines[key] = m + + source_data = { + 'x': [], 'y': [], 'color': [], 'marker': [], 'size': [], 'alpha': [], 'line_color': [], 'line_width': [], + 'p': [], 'type': [], 'd': [], 'q': [], 'r': [], 'op_type': [], 'mode': [], 'order': [], 'beam': [], 'is_baseline': [], + 'sparsify_errors': [], 'sparsify_reactivate_limit': [], 'beam_climbing': [], 'decoder': [], + 'label': [], 'hover_label': [], 'ler_str': [], 'speedup_str': [], 'mid_x': [], 'mid_y': [], 'base_x': [], 'base_y': [], + 'ler_err_low': [], 'ler_err_high': [], 'base_err_low': [], 'base_err_high': [], 'low_conf_str': [], 'shots': [] + } + + color_map = {'surfacecodes': '#5D95E8', 'colorcodes': '#F6C644', 'bivariatebicyclecodes': 'fuchsia'} + markers = ['circle', 'square', 'triangle', 'diamond', 'hex', 'star', 'inverted_triangle'] + + unique_qd = sorted(list(set((m['d'], m['q']) for m in metrics))) + marker_map = {unique_qd[i]: markers[i % len(markers)] for i in range(len(unique_qd))} + + short_names = {'surfacecodes': 'sc', 'colorcodes': 'cc', 'bivariatebicyclecodes': 'bb'} + type_display_names = {'surfacecodes': 'Surface Codes (sc)', 'colorcodes': 'Color Codes (cc)', 'bivariatebicyclecodes': 'Bivariate Bicycle Codes (bb)'} + + for m in metrics: + # Omit Mode B and Mode C completely from plotting + if m['mode'] in ['modeB', 'modeC']: + continue + + if m['is_baseline'] and (m.get('num_det_orders', 1) > 1 or m.get('decoder', 'tesseract') != 'tesseract'): + continue + + x1, y1 = m['time_per_round'], m['ler'] + if y1 <= 0 or x1 <= 0: + continue + + c = color_map.get(m['type'], 'black') + mark = marker_map.get((m['d'], m['q']), 'circle') + + order_val = "Ensemble of all" if m['order'] == 'all' else m['order'] + + source_data['x'].append(x1) + source_data['y'].append(y1) + source_data['color'].append(c) + source_data['marker'].append(mark) + + if m['is_baseline']: + source_data['line_color'].append('black') + source_data['line_width'].append(2) + source_data['alpha'].append(1.0) + if m['num_det_orders'] == 1: + source_data['size'].append(15) + else: + source_data['size'].append(20) + else: + source_data['line_color'].append(c) + source_data['line_width'].append(0) + source_data['size'].append(10) + source_data['alpha'].append(0.5) + + source_data['p'].append(str(m['p'])) + source_data['type'].append(m['type']) + source_data['d'].append(str(m['d'])) + source_data['q'].append(str(m['q'])) + source_data['r'].append(str(m['r'])) + source_data['op_type'].append('Auxiliary Columns' if m['op_type'] == 'normal' else 'Original Columns') + source_data['mode'].append(m['mode']) + source_data['order'].append(order_val) + source_data['decoder'].append(m.get('decoder', 'tesseract')) + if m['beam'] == 'simplex': + source_data['beam'].append('simplex') + elif m['beam'] == 20: + source_data['beam'].append('longbeam') + else: + source_data['beam'].append(str(m['beam'])) + source_data['is_baseline'].append(str(m['is_baseline'])) + source_data['sparsify_errors'].append(str(m['sparsify_errors'])) + source_data['sparsify_reactivate_limit'].append(str(m['sparsify_reactivate_limit'])) + source_data['beam_climbing'].append(str(m['beam_climbing'])) + source_data['ler_err_low'].append(m['ler_err_low']) + source_data['ler_err_high'].append(m['ler_err_high']) + total_failures = int(round(m['ler'] * m['shots'] * m['r'])) + frac = m['num_low_confidence'] / total_failures if total_failures > 0 else 0 + num_high_conf_errors = max(0, total_failures - m['num_low_confidence']) + ler_no_low_conf = num_high_conf_errors / (m['shots'] * m['r']) if (m['shots'] * m['r']) > 0 else 0 + source_data['low_conf_str'].append(f"{frac*100:.1f}% (LER w/o low-conf: {ler_no_low_conf:.2e})") + source_data['shots'].append(str(m['shots'])) + + label = "Baseline" if m['is_baseline'] else format_gari_label(m['mode'], m['order'], m['beam'], m['beam_climbing'], m['sparsify_errors'], m['sparsify_reactivate_limit'], m.get('decoder', 'tesseract')) + source_data['label'].append(label) + + short = short_names.get(m['type'], m['type']) + code_dist_str = f"{short}, d={m['d']}" + if m['is_baseline']: + hover_label = f"Baseline [{code_dist_str}]" + else: + hover_label = f"{label} [{code_dist_str}]" + source_data['hover_label'].append(hover_label) + + base_key = (m['p'], m['type'], m['d'], m['q'], m['r']) + if not m['is_baseline'] and base_key in baselines: + base_21 = baselines[base_key] + x0, y0 = base_21['time_per_round'], base_21['ler'] + + speedup = x0 / x1 if x1 > 0 else 1 + + k0 = y0 * base_21['r'] * base_21['shots'] + n0 = base_21['shots'] + k1 = y1 * m['r'] * m['shots'] + n1 = m['shots'] + + k0_a, n0_a = k0 + 0.5, n0 + 0.5 + k1_a, n1_a = k1 + 0.5, n1 + 0.5 + + r_adj = (k1_a / n1_a) / (k0_a / n0_a) + se_log_r = math.sqrt(max(0, 1/k1_a - 1/n1_a + 1/k0_a - 1/n0_a)) + + r_low = r_adj * math.exp(-1.96 * se_log_r) + r_high = r_adj * math.exp(1.96 * se_log_r) + + if round(r_low, 2) == round(r_high, 2): + ler_str = f"{r_low:.2f}x err" + else: + ler_str = f"{r_low:.2f}-{r_high:.2f}x err" + + mid_x = math.exp((math.log(x0) + math.log(x1)) / 2) if x0>0 and x1>0 else (x0+x1)/2 + mid_y = math.exp((math.log(y0) + math.log(y1)) / 2) if y0>0 and y1>0 else (y0+y1)/2 + + source_data['speedup_str'].append(f"{speedup:.1f}x spd") + source_data['ler_str'].append(ler_str) + source_data['mid_x'].append(mid_x) + source_data['mid_y'].append(mid_y) + source_data['base_x'].append(x0) + source_data['base_y'].append(y0) + source_data['base_err_low'].append(base_21['ler_err_low']) + source_data['base_err_high'].append(base_21['ler_err_high']) + else: + source_data['speedup_str'].append("") + source_data['ler_str'].append("") + source_data['mid_x'].append(x1) + source_data['mid_y'].append(y1) + source_data['base_x'].append(x1) + source_data['base_y'].append(y1) + source_data['base_err_low'].append(m['ler_err_low']) + source_data['base_err_high'].append(m['ler_err_high']) + + source = ColumnDataSource(data=source_data) + + line_source = ColumnDataSource(data={'xs': [], 'ys': []}) + text_source = ColumnDataSource(data={'x': [], 'y': [], 'text': []}) + tap_details_source = ColumnDataSource(data={'x': [], 'y': [], 'text': [], 'x_offset': [], 'y_offset': [], 'text_align': []}) + + # Error bar sources + tap_error_source = ColumnDataSource(data={'x': [], 'y0': [], 'y1': [], 'color': []}) + hover_error_source = ColumnDataSource(data={'x': [], 'y0': [], 'y1': [], 'color': []}) + selected_baseline_source = ColumnDataSource(data={'x': [], 'y': [], 'color': [], 'marker': []}) + + unique_decoders = sorted(list(set(source_data['decoder']))) + unique_types = sorted(list(set(source_data['type']))) + unique_modes = sorted(list(set(m for m in source_data['mode'] if m != 'unknown'))) + unique_orders = sorted(list(set(o for o in source_data['order'] if o != 'unknown'))) + unique_op_types = sorted(list(set(source_data['op_type']))) + unique_beams = sorted(list(set(b for b in source_data['beam'] if b != '0'))) + unique_p = sorted(list(set(source_data['p']))) + unique_beam_climbing = sorted(list(set(source_data['beam_climbing']))) + unique_sparsify = sorted(list(set(source_data['sparsify_errors']))) + + # Native Bokeh UI Controls for Distance + code_sections = [] + d_chk_widgets = {} + d_chk_keys = [] + + for ct in unique_types: + d_list = sorted(list(set(source_data['d'][i] for i in range(len(source_data['d'])) if source_data['type'][i] == ct))) + if not d_list: + continue + + name = type_display_names.get(ct, ct) + btn = Button(label=f"{name} ▼", button_type="default", width=330, height=30) + + d_rows = [] + family_widgets = {} + for d in d_list: + td_val = f"{ct}_{d}" + d_chk_keys.append(td_val) + + c = color_map.get(ct, 'black') + + # find marker + mark = 'circle' + for i in range(len(source_data['type'])): + if source_data['type'][i] == ct and source_data['d'][i] == d: + mark = source_data['marker'][i] + break + + sym = '●' + if mark == 'square': sym = '■' + elif mark == 'triangle': sym = '▲' + elif mark == 'diamond': sym = '◆' + elif mark == 'hex': sym = '⬢' + elif mark == 'star': sym = '★' + elif mark == 'inverted_triangle': sym = '▼' + + chk = CheckboxGroup(labels=[f"Distance {d}"], active=[0], width=100) + mark_div = Div(text=f"
{sym}
", width=30) + + d_chk_widgets[td_val] = chk + family_widgets[f"chk_{d}"] = chk + d_rows.append(row(chk, mark_div, sizing_mode="fixed", width=150, height=30)) + + select_all_btn = Button(label="Select All", button_type="success", width=100, height=30, sizing_mode="fixed") + clear_all_btn = Button(label="Clear All", button_type="primary", width=100, height=30, sizing_mode="fixed") + + js_select_all = "; ".join([f"{k}.active = [0]" for k in family_widgets.keys()]) + js_clear_all = "; ".join([f"{k}.active = []" for k in family_widgets.keys()]) + + select_all_btn.js_on_click(CustomJS(args=family_widgets, code=js_select_all)) + clear_all_btn.js_on_click(CustomJS(args=family_widgets, code=js_clear_all)) + + d_col = column(row(select_all_btn, clear_all_btn), *d_rows, visible=False) + btn.js_on_click(CustomJS(args=dict(col=d_col), code="col.visible = !col.visible;")) + code_sections.append(column(btn, d_col)) + + accordion_container = column(*code_sections) + + # General UI Controls + decoder_group = CheckboxGroup(labels=[d.capitalize() for d in unique_decoders], active=list(range(len(unique_decoders)))) + mode_group = CheckboxGroup(labels=unique_modes, active=list(range(len(unique_modes)))) + order_group = CheckboxGroup(labels=unique_orders, active=list(range(len(unique_orders)))) + op_type_group = CheckboxGroup(labels=unique_op_types, active=list(range(len(unique_op_types)))) + beam_group = CheckboxGroup(labels=unique_beams, active=list(range(len(unique_beams)))) + beam_climbing_group = CheckboxGroup(labels=[f"Beam Climbing: {v}" for v in unique_beam_climbing], active=list(range(len(unique_beam_climbing)))) + sparsify_group = CheckboxGroup(labels=[f"Sparsify: {v}" for v in unique_sparsify], active=list(range(len(unique_sparsify)))) + p_select = Select(title="Physical Error Rate (p)", value=unique_p[0] if unique_p else "", options=unique_p) + + mode_select_all = Button(label="Select All", button_type="success", width=70, height=30, sizing_mode="fixed") + mode_clear_all = Button(label="Clear All", button_type="primary", width=70, height=30, sizing_mode="fixed") + order_select_all = Button(label="Select All", button_type="success", width=70, height=30, sizing_mode="fixed") + order_clear_all = Button(label="Clear All", button_type="primary", width=70, height=30, sizing_mode="fixed") + + filter_args = dict( + decoder_group=decoder_group, unique_decoders=unique_decoders, + mode_group=mode_group, order_group=order_group, + op_type_group=op_type_group, beam_group=beam_group, + beam_climbing_group=beam_climbing_group, sparsify_group=sparsify_group, + p_select=p_select, + unique_modes=unique_modes, unique_orders=unique_orders, + unique_op_types=unique_op_types, unique_beams=unique_beams, + unique_beam_climbing=unique_beam_climbing, unique_sparsify=unique_sparsify, + d_chk_keys=d_chk_keys, d_chk_widgets=d_chk_widgets + ) + + js_code = """ + const indices = []; + const active_decoders = decoder_group.active.map(i => unique_decoders[i]); + const active_modes = mode_group.active.map(i => unique_modes[i]); + const active_orders = order_group.active.map(i => unique_orders[i]); + const active_op_types = op_type_group.active.map(i => unique_op_types[i]); + const active_beams = beam_group.active.map(i => unique_beams[i]); + const active_beam_climbing = beam_climbing_group.active.map(i => unique_beam_climbing[i]); + const active_sparsify = sparsify_group.active.map(i => unique_sparsify[i]); + const active_p = p_select.value; + + let active_td = []; + for (let i = 0; i < d_chk_keys.length; i++) { + let key = d_chk_keys[i]; + let chk = d_chk_widgets[key]; + if (chk.active.includes(0)) { + active_td.push(key); + } + } + + for (let i = 0; i < source.get_length(); i++) { + if (source.data['p'][i] !== active_p) continue; + if (!active_decoders.includes(source.data['decoder'][i])) continue; + + const t = source.data['type'][i]; + const d = source.data['d'][i]; + const td_key = t + "_" + d; + + if (!active_td.includes(td_key)) continue; + + const is_base = source.data['is_baseline'][i] === 'True'; + if (is_base) { + indices.push(i); + } else { + const dec = source.data['decoder'][i]; + if (dec === 'simplex') { + if (active_modes.includes(source.data['mode'][i]) && + active_op_types.includes(source.data['op_type'][i]) && + active_sparsify.includes(source.data['sparsify_errors'][i])) { + indices.push(i); + } + } else { + if (active_modes.includes(source.data['mode'][i]) && + active_orders.includes(source.data['order'][i]) && + active_op_types.includes(source.data['op_type'][i]) && + active_beams.includes(source.data['beam'][i]) && + active_beam_climbing.includes(source.data['beam_climbing'][i]) && + active_sparsify.includes(source.data['sparsify_errors'][i])) { + indices.push(i); + } + } + } + } + return indices; + """ + + js_filter = CustomJSFilter(args=filter_args, code=js_code) + view = CDSView(filter=js_filter) + + fastest_data = {'x': [], 'y': [], 'marker': []} + accurate_data = {'x': [], 'y': [], 'marker': []} + best_x = {} + best_y = {} + active_p = unique_p[0] if unique_p else "" + for i in range(len(source_data['x'])): + if source_data['p'][i] != active_p: + continue + if source_data['is_baseline'][i] == 'True': + continue + + t = source_data['type'][i] + d = source_data['d'][i] + key = f"{t}_{d}" + + if key not in best_x: + best_x[key] = {'min_x': float('inf'), 'x': -1, 'y': -1, 'marker': 'circle'} + best_y[key] = {'min_y': float('inf'), 'x': -1, 'y': -1, 'marker': 'circle'} + + x, y = source_data['x'][i], source_data['y'][i] + if x < best_x[key]['min_x']: + best_x[key]['min_x'] = x + best_x[key]['x'] = x + best_x[key]['y'] = y + best_x[key]['marker'] = source_data['marker'][i] + + if y < best_y[key]['min_y']: + best_y[key]['min_y'] = y + best_y[key]['x'] = x + best_y[key]['y'] = y + best_y[key]['marker'] = source_data['marker'][i] + + for k in best_x: + if best_x[k]['x'] != -1: + fastest_data['x'].append(best_x[k]['x']) + fastest_data['y'].append(best_x[k]['y']) + fastest_data['marker'].append(best_x[k]['marker']) + for k in best_y: + if best_y[k]['x'] != -1: + accurate_data['x'].append(best_y[k]['x']) + accurate_data['y'].append(best_y[k]['y']) + accurate_data['marker'].append(best_y[k]['marker']) + + fastest_source = ColumnDataSource(data=fastest_data) + accurate_source = ColumnDataSource(data=accurate_data) + + update_args = dict( + source=source, fastest_source=fastest_source, accurate_source=accurate_source, + decoder_group=decoder_group, unique_decoders=unique_decoders, + mode_group=mode_group, order_group=order_group, + op_type_group=op_type_group, beam_group=beam_group, + beam_climbing_group=beam_climbing_group, sparsify_group=sparsify_group, + p_select=p_select, + unique_modes=unique_modes, unique_orders=unique_orders, + unique_op_types=unique_op_types, unique_beams=unique_beams, + unique_beam_climbing=unique_beam_climbing, unique_sparsify=unique_sparsify, + d_chk_keys=d_chk_keys, d_chk_widgets=d_chk_widgets + ) + + update_js_code = """ + const active_decoders = decoder_group.active.map(i => unique_decoders[i]); + const active_modes = mode_group.active.map(i => unique_modes[i]); + const active_orders = order_group.active.map(i => unique_orders[i]); + const active_op_types = op_type_group.active.map(i => unique_op_types[i]); + const active_beams = beam_group.active.map(i => unique_beams[i]); + const active_beam_climbing = beam_climbing_group.active.map(i => unique_beam_climbing[i]); + const active_sparsify = sparsify_group.active.map(i => unique_sparsify[i]); + const active_p = p_select.value; + + let active_td = []; + for (let i = 0; i < d_chk_keys.length; i++) { + let key = d_chk_keys[i]; + let chk = d_chk_widgets[key]; + if (chk.active.includes(0)) { + active_td.push(key); + } + } + + let best_x = {}; + let best_y = {}; + + for (let i = 0; i < source.get_length(); i++) { + if (source.data['p'][i] !== active_p) continue; + if (!active_decoders.includes(source.data['decoder'][i])) continue; + + const t = source.data['type'][i]; + const d = source.data['d'][i]; + const key = t + "_" + d; + + if (!active_td.includes(key)) continue; + + const is_base = source.data['is_baseline'][i] === 'True'; + if (!is_base) { + const dec = source.data['decoder'][i]; + let is_active = false; + if (dec === 'simplex') { + is_active = (active_modes.includes(source.data['mode'][i]) && + active_op_types.includes(source.data['op_type'][i]) && + active_sparsify.includes(source.data['sparsify_errors'][i])); + } else { + is_active = (active_modes.includes(source.data['mode'][i]) && + active_orders.includes(source.data['order'][i]) && + active_op_types.includes(source.data['op_type'][i]) && + active_beams.includes(source.data['beam'][i]) && + active_beam_climbing.includes(source.data['beam_climbing'][i]) && + active_sparsify.includes(source.data['sparsify_errors'][i])); + } + + if (is_active) { + const x = source.data['x'][i]; + const y = source.data['y'][i]; + + if (!(key in best_x)) { + best_x[key] = {min_x: Infinity, x: -1, y: -1, marker: 'circle'}; + best_y[key] = {min_y: Infinity, x: -1, y: -1, marker: 'circle'}; + } + if (x < best_x[key].min_x) { + best_x[key].min_x = x; best_x[key].x = x; best_x[key].y = y; best_x[key].marker = source.data['marker'][i]; + } + if (y < best_y[key].min_y) { + best_y[key].min_y = y; best_y[key].x = x; best_y[key].y = y; best_y[key].marker = source.data['marker'][i]; + } + } + } + } + + const f_x = [], f_y = [], f_marker = []; + const a_x = [], a_y = [], a_marker = []; + + for (const k in best_x) { + if (best_x[k].x !== -1) { + f_x.push(best_x[k].x); f_y.push(best_x[k].y); f_marker.push(best_x[k].marker); + } + } + for (const k in best_y) { + if (best_y[k].x !== -1) { + a_x.push(best_y[k].x); a_y.push(best_y[k].y); a_marker.push(best_y[k].marker); + } + } + + fastest_source.data = {'x': f_x, 'y': f_y, 'marker': f_marker}; + accurate_source.data = {'x': a_x, 'y': a_y, 'marker': a_marker}; + fastest_source.change.emit(); + accurate_source.change.emit(); + + if (!window.original_legend_items) { + window.original_legend_items = legend.items.slice(); + } + + const new_legend_items = []; + for (let i = 0; i < window.original_legend_items.length; i++) { + const item = window.original_legend_items[i]; + if (item.renderers && item.renderers.length > 0) { + const r_name = item.renderers[0].name; + if (r_name && r_name.startsWith("fixed_")) { + new_legend_items.push(item); + } else if (r_name) { + if (active_td.includes(r_name) && active_decoders.length > 0) { + new_legend_items.push(item); + } + } else { + new_legend_items.push(item); + } + } else { + new_legend_items.push(item); + } + } + legend.items = new_legend_items; + + source.change.emit(); + """ + + update_js = CustomJS(args=update_args, code=update_js_code) + + for ctrl in [decoder_group, mode_group, order_group, op_type_group, beam_group, beam_climbing_group, sparsify_group]: + ctrl.js_on_change('active', update_js) + p_select.js_on_change('value', update_js) + + # bind all distance checkboxes natively + for td_val, chk in d_chk_widgets.items(): + chk.js_on_change('active', update_js) + + mode_select_all.js_on_click(CustomJS(args=dict(group=mode_group), code="group.active = Array.from({length: group.labels.length}, (v, k) => k);")) + mode_clear_all.js_on_click(CustomJS(args=dict(group=mode_group), code="group.active = [];")) + order_select_all.js_on_click(CustomJS(args=dict(group=order_group), code="group.active = Array.from({length: group.labels.length}, (v, k) => k);")) + order_clear_all.js_on_click(CustomJS(args=dict(group=order_group), code="group.active = [];")) + + p = figure(width=1000, height=800, title="GARI Tradeoffs", + x_axis_type="log", y_axis_type="log", + x_axis_label="Time per round (seconds)", y_axis_label="Logical Error Rate per round", + tools="pan,wheel_zoom,box_zoom,reset,tap,save") + p.xaxis.axis_line_color = "black" + p.yaxis.axis_line_color = "black" + p.xaxis.major_tick_line_color = "black" + p.yaxis.major_tick_line_color = "black" + + scatter = p.scatter(x='x', y='y', source=source, view=view, size='size', color='color', line_color='line_color', line_width='line_width', marker='marker', alpha='alpha', nonselection_alpha=0.1) + + p.scatter(x='x', y='y', source=fastest_source, marker='marker', size=20, fill_color=None, line_color='green', line_width=4, alpha=1.0) + p.scatter(x='x', y='y', source=accurate_source, marker='marker', size=20, fill_color=None, line_color='red', line_width=4, alpha=1.0) + p.scatter(x='x', y='y', source=selected_baseline_source, marker='marker', size=15, fill_color='color', line_color='black', line_width=2, alpha=1.0) + + p.multi_line(xs='xs', ys='ys', source=line_source, color="black", line_dash="dashed", line_width=2, alpha=0.6) + p.segment(x0='x', y0='y0', x1='x', y1='y1', source=tap_error_source, color='color', line_width=2, alpha=0.8) + p.segment(x0='x', y0='y0', x1='x', y1='y1', source=hover_error_source, color='color', line_width=3, alpha=0.8) + + p.text(x='x', y='y', text='text', source=text_source, text_font_size="9pt", + x_offset=5, y_offset=5, text_baseline="bottom") + + details_box = LabelSet( + x='x', y='y', text='text', source=tap_details_source, + text_font_size="9pt", level="annotation", + x_offset='x_offset', y_offset='y_offset', + text_align='text_align', + border_line_color="#cccccc", border_line_alpha=1.0, + background_fill_color="#ffffff", background_fill_alpha=0.95 + ) + p.add_layout(details_box) + + # Dummy Glyphs for Combination Legend + unique_combos = sorted(list(set( + (source_data['type'][i], int(source_data['d'][i]), int(source_data['q'][i])) + for i in range(len(source_data['type'])) + )), key=lambda x: (x[0], x[1], x[2])) + + legend_renderers = {} + for combo in unique_combos: + c_type, d, q = combo + color = color_map.get(c_type, 'black') + marker = marker_map.get((d, q), 'circle') + short = short_names.get(c_type, c_type) + + if c_type == 'bivariatebicyclecodes': + label = f"{short}, d={d}, q={q}" + else: + label = f"{short}, d={d}" + + renderer = p.scatter(x=[float('nan')], y=[float('nan')], fill_color=color, line_color=color, marker=marker, size=10, legend_label=label, name=f"{c_type}_{d}") + legend_renderers[f"{c_type}_{d}"] = renderer + + p.scatter(x=[float('nan')], y=[float('nan')], fill_color='gray', line_color='black', line_width=1.5, marker='circle', size=12, legend_label='Baseline (Black Border)', name="fixed_baseline") + p.scatter(x=[float('nan')], y=[float('nan')], fill_color='white', line_color='green', line_width=4, size=12, legend_label='Fastest Gari', name="fixed_fastest") + p.scatter(x=[float('nan')], y=[float('nan')], fill_color='white', line_color='red', line_width=4, size=12, legend_label='Most Accurate Gari', name="fixed_accurate") + + update_js.args['legend'] = p.legend[0] + + p.add_layout(p.legend[0], 'right') + p.legend.click_policy = "hide" + + hover = HoverTool(renderers=[scatter], tooltips=[ + ("Config", "@hover_label"), + ("Decoder", "@decoder"), + ("Time", "@x"), + ("LER", "@y"), + ("Low Conf Fraction of Errors", "@low_conf_str") + ]) + p.add_tools(hover) + + hover_js = CustomJS(args=dict(source=source, hover_error_source=hover_error_source), code=""" + if (cb_data && cb_data.index && cb_data.index.indices.length > 0) { + const i = cb_data.index.indices[0]; + const x = source.data['x'][i]; + const y = source.data['y'][i]; + const err_low = source.data['ler_err_low'][i]; + const err_high = source.data['ler_err_high'][i]; + const c = source.data['color'][i]; + hover_error_source.data = { + 'x': [x], + 'y0': [y - err_low], + 'y1': [y + err_high], + 'color': [c] + }; + hover_error_source.change.emit(); + } else { + hover_error_source.data = {'x': [], 'y0': [], 'y1': [], 'color': []}; + hover_error_source.change.emit(); + } + """) + hover.callback = hover_js + + details_panel = Div( + text="

Selected Point Details

Click on a data point to see details here.

", + width=1000, + styles={"padding": "10px", "border": "1px solid #ccc", "background-color": "#f9f9f9", "margin-top": "10px"} + ) + + tap_js = CustomJS(args=dict(source=source, line_source=line_source, text_source=text_source, tap_error_source=tap_error_source, tap_details_source=tap_details_source, details_panel=details_panel, details_box=details_box, selected_baseline_source=selected_baseline_source), code=""" + const sel = source.selected.indices; + const line_xs = []; + const line_ys = []; + const text_x = []; + const text_y = []; + const texts = []; + const err_x = []; + const err_y0 = []; + const err_y1 = []; + const err_color = []; + const det_x = []; + const det_y = []; + const det_text = []; + const det_xoff = []; + const det_yoff = []; + const det_align = []; + const base_xs = []; + const base_ys = []; + const base_colors = []; + const base_markers = []; + let panel_html = "

Selected Points Details

"; + + if (sel.length > 0) { + for (let k = 0; k < sel.length; k++) { + const i = sel[k]; + + const decoder = source.data['decoder'][i]; + const ler = source.data['y'][i].toExponential(4); + const beam = source.data['beam'][i]; + const order = source.data['order'][i]; + const sp_err = source.data['sparsify_errors'][i]; + const sp_lim = source.data['sparsify_reactivate_limit'][i]; + const shots = source.data['shots'][i]; + const mode = source.data['mode'][i]; + const ler_err_low = source.data['ler_err_low'][i].toExponential(2); + const ler_err_high = source.data['ler_err_high'][i].toExponential(2); + const low_conf = source.data['low_conf_str'][i]; + const hover_lbl = source.data['hover_label'][i]; + const x1 = source.data['x'][i]; + const y1 = source.data['y'][i]; + const time_val = x1.toExponential(4); + const ler_val = y1.toExponential(4); + + panel_html += `

Point ${k+1}: ${hover_lbl}

+ `; + + const details_text = `Config: ${hover_lbl}\\nDecoder: ${decoder}\\nTime: ${time_val}s\\nLER: ${ler_val} (-${ler_err_low}, +${ler_err_high})\\nLow Conf: ${low_conf}`; + + const is_base = source.data['is_baseline'][i] === 'True'; + + let is_faster = false; + if (!is_base && source.data['speedup_str'][i] !== "") { + const x0 = source.data['base_x'][i]; + if (x1 < x0) { + is_faster = true; + } + } + + det_x.push(x1); + det_y.push(y1); + det_text.push(details_text); + det_yoff.push(12); + + if (is_faster) { + det_xoff.push(-12); + det_align.push('right'); + } else { + det_xoff.push(12); + det_align.push('left'); + } + + if (!is_base && source.data['speedup_str'][i] !== "") { + const x0 = source.data['base_x'][i]; + const y0 = source.data['base_y'][i]; + const spd = source.data['speedup_str'][i]; + const ler_str = source.data['ler_str'][i]; + const mid_x = source.data['mid_x'][i]; + const mid_y = source.data['mid_y'][i]; + + line_xs.push([x0, x1]); + line_ys.push([y0, y1]); + + text_x.push(mid_x); + text_y.push(mid_y); + texts.push(spd + "\\n" + ler_str); + + const err0_low = source.data['base_err_low'][i]; + const err0_high = source.data['base_err_high'][i]; + const err1_low = source.data['ler_err_low'][i]; + const err1_high = source.data['ler_err_high'][i]; + const c = source.data['color'][i]; + + err_x.push(x0, x1); + err_y0.push(y0 - err0_low, y1 - err1_low); + err_y1.push(y0 + err0_high, y1 + err1_high); + err_color.push(c, c); + + base_xs.push(x0); + base_ys.push(y0); + base_colors.push(source.data['color'][i]); + base_markers.push(source.data['marker'][i]); + } + } + + line_source.data = {'xs': line_xs, 'ys': line_ys}; + text_source.data = {'x': text_x, 'y': text_y, 'text': texts}; + tap_error_source.data = {'x': err_x, 'y0': err_y0, 'y1': err_y1, 'color': err_color}; + tap_details_source.data = {'x': det_x, 'y': det_y, 'text': det_text, 'x_offset': det_xoff, 'y_offset': det_yoff, 'text_align': det_align}; + selected_baseline_source.data = {'x': base_xs, 'y': base_ys, 'color': base_colors, 'marker': base_markers}; + + line_source.change.emit(); + text_source.change.emit(); + tap_error_source.change.emit(); + tap_details_source.change.emit(); + selected_baseline_source.change.emit(); + + details_panel.text = panel_html; + } else { + details_panel.text = "

Selected Point Details

Click on a data point to see details here.

"; + line_source.data = {'xs': [], 'ys': []}; + text_source.data = {'x': [], 'y': [], 'text': []}; + tap_error_source.data = {'x': [], 'y0': [], 'y1': [], 'color': []}; + tap_details_source.data = {'x': [], 'y': [], 'text': [], 'x_offset': [], 'y_offset': [], 'text_align': []}; + selected_baseline_source.data = {'x': [], 'y': [], 'color': [], 'marker': []}; + + line_source.change.emit(); + text_source.change.emit(); + tap_error_source.change.emit(); + tap_details_source.change.emit(); + selected_baseline_source.change.emit(); + } + """) + + tap = p.select(type=TapTool) + source.selected.js_on_change('indices', tap_js) + + controls = column( + p_select, + Div(text="Decoder"), + decoder_group, + Div(text="Code Types & Distances"), + accordion_container, + Div(text="Prior Modes"), + row(mode_select_all, mode_clear_all), + mode_group, + Div(text="Detector Orders"), + row(order_select_all, order_clear_all), + order_group, + Div(text="Operation Types"), + op_type_group, + Div(text="Beams"), + beam_group, + Div(text="Beam Climbing"), + beam_climbing_group, + Div(text="Sparsify Errors"), + sparsify_group, + width=350 + ) + + MODE_DESC = { + 'modeA': "ez,ex,ey Original (Keep), ez',ex' Aggregated", + 'modeD': "ez,ex,ey Free, ez',ex' Penalized", + 'modeE': "ez,ex,ey Scaled, ez',ex' Penalized", + 'modeF': "ez,ex,ey Free, ez',ex' Aggregated", + 'modeG': "ez,ex,ey Scaled to make the original columns almost free cost, ez',ex' Aggregated", + 'modeH': "ez,ex,ey Keep, ez',ex' Keep", + 'modeI': "ez,ex,ey Keep, ez',ex' Maxed", + 'modeJ': "ez,ex,ey Free, ez',ex' Maxed", + 'modeK': "ez,ex,ey Scaled, ez',ex' Maxed", + 'modeL': "ez,ex Keep, ey Free, ez',ex' Aggregated", + 'modeM': "ez,ex Keep, ey Free, ez',ex' Maxed", + 'modeN': "ez,ex,ey Keep, ez',ex' XOR Aggregated", + 'modeO': "ez,ex,ey Scaled, ez',ex' XOR Aggregated", + 'modeP': "Minimax Cost-Split (Exact mathematical split, fixes double counting)", + 'modeS': "LP Mode", + 'modeS2': "LP Mode", + 'modeSO': "LP Mode", + 'modeSO2': "LP Mode", + 'modeU': "LP Mode", + 'modeV': "LP Mode" + } + + ORDER_DESC = { + 'order1': "RealX, VirtX, RealZ, VirtZ", + 'order2': "RealX, RealZ, VirtX, VirtZ", + 'order3': "RealX, VirtZ, RealZ, VirtX", + 'order4': "RealX, RealZ, VirtZ, VirtX", + 'order4r': "VirtZ, VirtX, RealX, RealZ (Virtual First)", + 'order5': "First-Touch Interleaved Chrono", + 'order6': "Last-Touch Interleaved Chrono", + 'order7': "Chrono Real, Chrono Virt", + 'order7r': "Chrono Virt, Chrono Real (Virtual First)", + 'order7a': "Chrono Real, Chrono Virt (With ey max)", + 'order7ar': "Chrono Virt, Chrono Real (Virtual First, With ey max)", + 'order7b': "Chrono Real, Chrono Virt min (No ey)", + 'order7br': "Chrono Virt min, Chrono Real (Virtual First, No ey)", + 'order7c': "Chrono Real, Chrono Virt min (With ey)", + 'order7cr': "Chrono Virt min, Chrono Real (Virtual First, With ey)", + 'order8': "RealX, VirtX, VirtZ, RealZ", + 'order9': "RealZ, RealX, VirtZ, VirtX", + 'order9r': "VirtZ, VirtX, RealZ, RealX (Virtual First)", + 'order10': "RealZ, RealX, VirtX, VirtZ", + 'order10r': "VirtX, VirtZ, RealZ, RealX (Virtual First)", + 'Ensemble of all': "A combination of all multiple detector orders" + } + + modes_html = "" + for m in unique_modes: + desc = MODE_DESC.get(m, "Unknown Mode") + if desc == "Unknown Mode" and (m.startswith("mode") or m.startswith("LP")): + desc = "LP Mode" + modes_html += f"
  • {m}: {desc}
  • \n" + + orders_html = "" + for o in unique_orders: + desc = ORDER_DESC.get(o, "Unknown Order") + orders_html += f"
  • {o}: {desc}
  • \n" + + manual_div_str = f""" +
    + Legend & Manual
    + +
    +
    + Prior Modes: +
      + {modes_html} +
    +
    +
    + Detector Orders: +
      + {orders_html} +
    +
    +
    +
    + """ + manual_div = Div(text=manual_div_str) + + layout = column(row(controls, p), manual_div, details_panel) + + output_dir = 'plots' + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + output_file(os.path.join(output_dir, "interactive_plot.html"), title="GARI Interactive Plot") + save(layout) + +if __name__ == '__main__': + generate_interactive_plot() + print("Successfully generated plots/interactive_plot.html") diff --git a/benchmarking/gari/submit.sh b/benchmarking/gari/submit.sh new file mode 100644 index 0000000..7ae6a06 --- /dev/null +++ b/benchmarking/gari/submit.sh @@ -0,0 +1,166 @@ +#!/usr/bin/env bash + +# Exit immediately if a command exits with a non-zero status (-e), +# treat unset variables as an error (-u), and catch errors in pipes (-o pipefail). +set -euo pipefail + +mkdir -p out + +TESSERACT_BIN=./bazel-bin/src/tesseract +SIMPLEX_BIN=./bazel-bin/src/simplex +# Create a timestamp in nanoseconds for when the script starts +STARTTIME=$(($(date +%s%N))) + +COUNTER=0 + +for num in $(seq 0 20); do + for p_err in 0.001; do + # for circuit in testdata/bivariatebicyclecodes/r=6,*p=$p_err,noise=si1000,c=bivariate_bicycle_Z*.stim; do + # echo "$circuit" + # done + for circuit in testdata/colorcodes/r=7,*p=$p_err,noise=si1000,c=superdense_color_code_Z,*cz.stim; do + echo "$circuit" + done + # for circuit in testdata/surfacecodes/r=9,*p=$p_err,noise=si1000,c=surface_code_Z,*cz.stim; do + # echo "$circuit" + # done + done +done | shuf | while read circuit; do + circuit_dir=$(dirname "$circuit") + circuit_name=$(basename "$circuit" .stim) + mapping_file="$circuit_dir/gari/${circuit_name}_mapping.json" + + echo "=========================================" + echo "Running benchmark for circuit: $circuit_name" + +# Determine base degree based on the folder/filename + if [[ "$circuit" == *"surfacecodes"* ]]; then + SPARSIFY_BASE_DEGREE=3 + else + SPARSIFY_BASE_DEGREE=4 + fi + # Simplex Baseline + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$SIMPLEX_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --stats-out out/${STARTTIME}-${COUNTER}-simplex-baseline.json" + # COUNTER=$((COUNTER + 1)) + + # Baseline 1 + sbatch --partition=c2 --job-name=gari \ + --ntasks=1 \ + --mem=120gb \ + --cpus-per-task=30 \ + --time=20:00:00 \ + --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --no-revisit-dets --beam 20 --beam-climbing --num-det-orders 1 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-1det.json" + COUNTER=$((COUNTER + 1)) + + # # Baseline 2 + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --no-revisit-dets --beam 20 --beam-climbing --num-det-orders 21 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-21det.json" + # COUNTER=$((COUNTER + 1)) + + # # Baseline 3 + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --no-revisit-dets --beam 5 --beam-climbing --num-det-orders 1 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-5beam-1det.json" + # COUNTER=$((COUNTER + 1)) + + # # Baseline 4 + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --no-revisit-dets --beam 5 --beam-climbing --num-det-orders 21 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-5beam-21det.json" + # COUNTER=$((COUNTER + 1)) + + # GARI Runs + for L_type in "ogL_"; do + for mode in modeN modeQ modeR modeS modeS2 modeSO modeSO2; do + dem_file="$circuit_dir/gari/${circuit_name}_${L_type}${mode}.dem" + echo "Running GARI mode: $dem_file" + + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$SIMPLEX_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --stats-out out/${STARTTIME}-${COUNTER}-simplex-gari-${L_type}${mode}.json" + # COUNTER=$((COUNTER + 1)) + + for order in order10 all; do + echo " Running order: $order" + + sbatch --partition=c2 --job-name=gari \ + --ntasks=1 \ + --mem=120gb \ + --cpus-per-task=30 \ + --time=20:00:00 \ + --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 5 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam5.json" + COUNTER=$((COUNTER + 1)) + + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 20 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam20.json" + # COUNTER=$((COUNTER + 1)) + + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 5 --num-det-orders 1 --pqlimit 1000000 --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam5_nc.json" + # COUNTER=$((COUNTER + 1)) + + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 10 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam10.json" + # COUNTER=$((COUNTER + 1)) + + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 10 --num-det-orders 1 --pqlimit 1000000 --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam10_nc.json" + # COUNTER=$((COUNTER + 1)) + + # same with sparcifacation + # for SPARSIFY_REACTIVATE_LIMIT in 0 2 4 8 16 32 64 128 256; do + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 5 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --sparsify-errors --sparsify-base-degree $SPARSIFY_BASE_DEGREE --sparsify-reactivate-limit $SPARSIFY_REACTIVATE_LIMIT --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam5_sparsify.json" + # COUNTER=$((COUNTER + 1)) + + # sbatch --partition=c2 --job-name=gari \ + # --ntasks=1 \ + # --mem=120gb \ + # --cpus-per-task=30 \ + # --time=20:00:00 \ + # --wrap="$TESSERACT_BIN --circuit \"$circuit\" --sample-num-shots 10000 --max-errors 10 --threads 30 --beam 5 --num-det-orders 1 --pqlimit 1000000 --sparsify-errors --sparsify-base-degree $SPARSIFY_BASE_DEGREE --sparsify-reactivate-limit $SPARSIFY_REACTIVATE_LIMIT --dem \"$dem_file\" --det-mapping-file \"$mapping_file\" --custom-order \"$order\" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam5_nc_sparsify.json" + # COUNTER=$((COUNTER + 1)) + # done + + done + done + done +done diff --git a/benchmarking/gari/submit_locally.sh b/benchmarking/gari/submit_locally.sh new file mode 100644 index 0000000..54072aa --- /dev/null +++ b/benchmarking/gari/submit_locally.sh @@ -0,0 +1,82 @@ +#!/usr/bin/env bash + +# Exit immediately if a command exits with a non-zero status (-e), +# treat unset variables as an error (-u), and catch errors in pipes (-o pipefail). +set -euo pipefail + +mkdir -p out + +TESSERACT_BIN=./bazel-bin/src/tesseract +SIMPLEX_BIN=./bazel-bin/src/simplex +# Create a timestamp in nanoseconds for when the script starts +STARTTIME=$(($(date +%s%N))) + +COUNTER=0 +SHOTS=100 +SEED=2 + +for num in 0; do + for p_err in 0.001; do + # for circuit in testdata/bivariatebicyclecodes/r=6,*p=$p_err,noise=si1000,c=bivariate_bicycle_Z,*.stim; do + # echo "$circuit" + # done + for circuit in testdata/colorcodes/r=7,*p=$p_err,noise=si1000,c=superdense_color_code_Z,*.stim; do + echo "$circuit" + done + # for circuit in testdata/surfacecodes/r=7,*p=$p_err,noise=si1000,c=surface_code_Z,*.stim; do + # echo "$circuit" + # done + done +done | shuf | while read circuit; do + circuit_dir=$(dirname "$circuit") + circuit_name=$(basename "$circuit" .stim) + mapping_file="$circuit_dir/gari/${circuit_name}_mapping.json" + + echo "=========================================" + echo "Running benchmark for circuit: $circuit_name" + # Submit also one baseline job + $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --no-revisit-dets --beam 20 --beam-climbing --num-det-orders 1 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-1det.json + $SIMPLEX_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --stats-out out/${STARTTIME}-${COUNTER}-simplex-baseline.json + COUNTER=$((COUNTER + 1)) + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --no-revisit-dets --beam 20 --beam-climbing --num-det-orders 21 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-21det.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --no-revisit-dets --beam 5 --beam-climbing --num-det-orders 1 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-5beam.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --no-revisit-dets --beam 5 --beam-climbing --num-det-orders 21 --det-order-index --pqlimit 1000000 --stats-out out/${STARTTIME}-${COUNTER}-baseline-5beam.json + + #run with gari + for L_type in ""; do + for mode in modeN modeQ modeR modeS modeS2 modeSO modeSO2; do + dem_file="$circuit_dir/gari/${circuit_name}_${L_type}${mode}.dem" + echo "Running GARI mode: $dem_file" + $SIMPLEX_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --dem "$dem_file" --det-mapping-file "$mapping_file" --stats-out out/${STARTTIME}-${COUNTER}-simplex-gari-${L_type}${mode}.json + COUNTER=$((COUNTER + 1)) + for order in order10 all; do + echo " Running order: $order" + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 5 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam5.json + # # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 20 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam20.json + # # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 5 --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam5_nc.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 2 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam2.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 3 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam3.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 4 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam4.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 6 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit_beam6.json + + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 20 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 20 --num-det-orders 1 --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 20 --num-det-orders 1 --pqlimit 10000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 5 --num-det-orders 1 --pqlimit 10000 --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + # for SPARSIFY_REACTIVATE_LIMIT in 0 2 4 8 16 32 64 128; do + # echo "Limit: $SPARSIFY_REACTIVATE_LIMIT" + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 5 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --sparsify-errors --sparsify-base-degree 4 --sparsify-reactivate-limit $SPARSIFY_REACTIVATE_LIMIT --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + # # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 10 --beam-climbing --num-det-orders 1 --pqlimit 1000000 --sparsify-errors --sparsify-base-degree 4 --sparsify-reactivate-limit $SPARSIFY_REACTIVATE_LIMIT --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 5 --num-det-orders 1 --pqlimit 1000000 --sparsify-errors --sparsify-base-degree 4 --sparsify-reactivate-limit $SPARSIFY_REACTIVATE_LIMIT --dem "$dem_file" --det-mapping-file "$mapping_file" --custom-order "$order" --stats-out out/${STARTTIME}-${COUNTER}-gari-${L_type}${mode}-${order}-revisit.json + + done + done + done + done + # $TESSERACT_BIN --circuit "$circuit" --sample-num-shots $SHOTS --sample-seed $SEED --max-errors 100 --threads 32 --beam 20 --beam-climbing --num-det-orders 21 --det-order-index --pqlimit 1000000 --dem "$dem_file" --det-mapping-file "$mapping_file" --stats-out out/${STARTTIME}-${COUNTER}-gari-revisit-21det.json + +# $TESSERACT_BIN --circuit "$circuit" --sample-num-shots 1000 --sample_seed $SEED --max-errors 10 --threads 32 --no-revisit-dets --beam 20 --beam-climbing --num-det-orders 1 --det-order-index --pqlimit 1000000 --dem "$dem_file" —det-mapping-file "$mapping_file" —stats-out out/${STARTTIME}-${COUNTER}.json + + # Increment counter for every single job so JSON files don't get overwritten + COUNTER=$((COUNTER + 1)) +done diff --git a/src/py/BUILD b/src/py/BUILD index e0bf9d8..5d14057 100644 --- a/src/py/BUILD +++ b/src/py/BUILD @@ -140,3 +140,18 @@ compile_pip_requirements( src = "requirements.in", requirements_txt = "requirements_lock.txt", ) + +py_binary( + name = "gari_simulation_test", + srcs = ["gari_simulation_test.py"], + visibility = ["//visibility:public"], + deps = [ + "//src/py/_tesseract_py_util:_tesseract_py_util", + "@pypi//numpy", + "@pypi//stim", + "@pypi//scipy", + "@pypi//matplotlib", + "//src:lib_tesseract_decoder", + ], + imports = ["..", "."], +) diff --git a/src/py/_tesseract_py_util/BUILD b/src/py/_tesseract_py_util/BUILD index 59f2131..a4ccbca 100644 --- a/src/py/_tesseract_py_util/BUILD +++ b/src/py/_tesseract_py_util/BUILD @@ -1,5 +1,6 @@ load("@rules_python//python:py_test.bzl", "py_test") load("@rules_python//python:py_library.bzl", "py_library") +load("@rules_python//python:py_binary.bzl", "py_binary") py_library( name = "_tesseract_py_util", @@ -8,6 +9,19 @@ py_library( deps = [ "@pypi//stim", "@pypi//numpy", + "@pypi//scipy", + ], +) + +py_binary( + name = "gari_dem_utils", + srcs = ["gari_dem_utils.py"], + visibility = ["//visibility:public"], + deps = [ + ":_tesseract_py_util", + "@pypi//scipy", + "@pypi//numpy", + "@pypi//stim", ], ) diff --git a/src/py/_tesseract_py_util/__init__.py b/src/py/_tesseract_py_util/__init__.py index fe103fe..e4276b8 100644 --- a/src/py/_tesseract_py_util/__init__.py +++ b/src/py/_tesseract_py_util/__init__.py @@ -20,3 +20,11 @@ from _tesseract_py_util.demutil import decompose_errors from _tesseract_py_util.generalize_dem import \ generalize as regeneralize_spatial_dem + +from .gari_dem_utils import ( + get_detector_types, + dem_to_check_matrices, + matrices_to_dem, + gari_transform, + get_detector_orderings, assign_prior_weights +) diff --git a/src/py/_tesseract_py_util/gari_dem_utils.py b/src/py/_tesseract_py_util/gari_dem_utils.py new file mode 100644 index 0000000..bb9fe31 --- /dev/null +++ b/src/py/_tesseract_py_util/gari_dem_utils.py @@ -0,0 +1,821 @@ +import os +import glob +import json +import argparse +import pathlib +import stim +import numpy as np +from scipy.sparse import csc_matrix + +def get_target_path(relative_path): + workspace_root = os.environ.get("BUILD_WORKSPACE_DIRECTORY", "") + return os.path.join(workspace_root, relative_path) + +def get_detector_types(circuit: stim.Circuit): + coords = circuit.get_detector_coordinates() + n_det = circuit.num_detectors + det_types = np.zeros(n_det, dtype=int) + for i in range(n_det): + if coords.get(i)[3] >=3: + # 3 means Z detector, 1 means X detector + det_types[i] = 3 + else: + det_types[i] = 1 # default + return det_types + +def dem_to_check_matrices(dem: stim.DetectorErrorModel, allow_undecomposed_hyperedges=True): + errors = [] + observables_list = [] + priors = [] + for inst in dem.flattened(): + if inst.type == "error": + priors.append(inst.args_copy()[0]) + targets = inst.targets_copy() + dets = [t.val for t in targets if t.is_relative_detector_id()] + obs = [t.val for t in targets if t.is_logical_observable_id()] + errors.append(dets) + observables_list.append(obs) + + M = dem.num_detectors + N = len(errors) + O = dem.num_observables + + row_ind = [] + col_ind = [] + data = [] + for j, dets in enumerate(errors): + for d in dets: + row_ind.append(d) + col_ind.append(j) + data.append(1) + + H = csc_matrix((data, (row_ind, col_ind)), shape=(M, N), dtype=np.uint8) + + obs_row_ind = [] + obs_col_ind = [] + obs_data = [] + for j, obs in enumerate(observables_list): + for o in obs: + obs_row_ind.append(o) + obs_col_ind.append(j) + obs_data.append(1) + + L = csc_matrix((obs_data, (obs_row_ind, obs_col_ind)), shape=(O, N), dtype=np.uint8) + + return H, L, np.array(priors), errors + +def matrices_to_dem(H: csc_matrix, L: csc_matrix, priors: np.ndarray) -> stim.DetectorErrorModel: + dem = stim.DetectorErrorModel() + H_csc = H.tocsc() + L_csc = L.tocsc() + + for j in range(H_csc.shape[1]): + p = priors[j] + targets = [] + + # Add detectors + for i in H_csc.indices[H_csc.indptr[j]:H_csc.indptr[j+1]]: + targets.append(stim.target_relative_detector_id(int(i))) + + # Add observables + for o in L_csc.indices[L_csc.indptr[j]:L_csc.indptr[j+1]]: + targets.append(stim.target_logical_observable_id(int(o))) + + if targets: + dem.append("error", p, targets) + + return dem + + + +def gari_transform(H: csc_matrix, L: csc_matrix, det_types: np.ndarray) -> dict: + """ + Applies the Gari Transform (arXiv:2510.14060) to the check matrix. + Splits Y errors into independent X and Z components and adds virtual detectors. + """ + is_x_det = (det_types == 1) + is_z_det = (det_types == 3) + + x_orig_indices = np.where(is_x_det)[0] + z_orig_indices = np.where(is_z_det)[0] + + H_csr = H.tocsr() + hx = H_csr[is_x_det, :] + hz = H_csr[is_z_det, :] + + hx_nnz = hx.getnnz(axis=0) + hz_nnz = hz.getnnz(axis=0) + + hx_any = hx_nnz > 0 + hz_any = hz_nnz > 0 + + i_hy = np.where(hx_any & hz_any)[0] + i_hx_only = np.where(hx_any & ~hz_any)[0] + i_hz_only = np.where(~hx_any & hz_any)[0] + + hx_csc = hx.tocsc() + hz_csc = hz.tocsc() + L_csc = L.tocsc() + + dx = hx_csc[:, i_hx_only] + dz = hz_csc[:, i_hz_only] + hx_yonly = hx_csc[:, i_hy] + hz_yonly = hz_csc[:, i_hy] + + L_dx = L_csc[:, i_hx_only] + L_dz = L_csc[:, i_hz_only] + L_y = L_csc[:, i_hy] + + mx, nx = dx.shape + mz, nz = dz.shape + ny = hx_yonly.shape[1] + + def get_col_hashes(mat_csc): + hashes = [] + for j in range(mat_csc.shape[1]): + start = mat_csc.indptr[j] + end = mat_csc.indptr[j+1] + hashes.append(tuple(mat_csc.indices[start:end])) + return hashes + + dx_hashes = get_col_hashes(dx) + dz_hashes = get_col_hashes(dz) + + dx_hash_to_idx = {h: i for i, h in enumerate(dx_hashes)} + dz_hash_to_idx = {h: i for i, h in enumerate(dz_hashes)} + + hx_yonly_hashes = get_col_hashes(hx_yonly) + hz_yonly_hashes = get_col_hashes(hz_yonly) + + U_rows, U_cols = [], [] + for j, h in enumerate(hx_yonly_hashes): + i = dx_hash_to_idx.get(h, -1) + if i != -1: + U_rows.append(i) + U_cols.append(j) + U_data = np.ones(len(U_rows), dtype=np.uint8) + U = csc_matrix((U_data, (U_rows, U_cols)), shape=(nx, ny), dtype=np.uint8) + + V_rows, V_cols = [], [] + for j, h in enumerate(hz_yonly_hashes): + i = dz_hash_to_idx.get(h, -1) + if i != -1: + V_rows.append(i) + V_cols.append(j) + V_data = np.ones(len(V_rows), dtype=np.uint8) + V = csc_matrix((V_data, (V_rows, V_cols)), shape=(nz, ny), dtype=np.uint8) + + from scipy.sparse import eye, bmat + + I_nx = eye(nx, format='csc', dtype=np.uint8) + I_nz = eye(nz, format='csc', dtype=np.uint8) + + blocks = [ + [None, None, None, dx, None], + [None, None, None, None, dz ], + [I_nx, None, U, I_nx, None], + [None, I_nz, V, None, I_nz] + ] + + gari_matrix = bmat(blocks, format='csc', dtype=np.uint8) + + L_ez = L_dx.astype(np.uint8) + L_ex = L_dz.astype(np.uint8) + L_ey = L_y.astype(np.uint8) + + O_num = L.shape[0] + Z_nx = csc_matrix((O_num, nx), dtype=np.uint8) + Z_nz = csc_matrix((O_num, nz), dtype=np.uint8) + Z_ny = csc_matrix((O_num, ny), dtype=np.uint8) + + L_blocks = [ + [Z_nx, Z_nz, Z_ny, L_dx.astype(np.uint8), L_dz.astype(np.uint8)] + ] + gari_obs_matrix = bmat(L_blocks, format='csc', dtype=np.uint8) + + L_blocks_og = [ + [L_dx.astype(np.uint8), L_dz.astype(np.uint8), L_ey, Z_nx, Z_nz] + ] + gari_obs_matrix_og = bmat(L_blocks_og, format='csc', dtype=np.uint8) + + gari_structure = { + "gari_matrix": gari_matrix, + "gari_obs_matrix": gari_obs_matrix, + "gari_obs_matrix_og": gari_obs_matrix_og, + "U": U, + "V": V, + "nx_virt": nx, + "nz_virt": nz, + "nx_real": mx, + "nz_real": mz, + "num_original_detectors": H.shape[0], + "num_gari_detectors": gari_matrix.shape[0], + "i_hx_only": i_hx_only, + "i_hz_only": i_hz_only, + "i_hy": i_hy, + "dx": dx, + "dz": dz, + "hx_csc": hx_csc, + "hz_csc": hz_csc, + "x_orig_indices": x_orig_indices, + "z_orig_indices": z_orig_indices, + } + return gari_structure + + +def assign_prior_weights(gari_structure: dict, method: str, original_priors: np.ndarray) -> np.ndarray: + i_hx_only = gari_structure["i_hx_only"] + i_hz_only = gari_structure["i_hz_only"] + i_hy = gari_structure["i_hy"] + + P_dx = original_priors[i_hx_only] + P_dz = original_priors[i_hz_only] + P_y = original_priors[i_hy] + + P_ez = P_dx + P_ex = P_dz + P_ey = P_y + + U = gari_structure["U"] + V = gari_structure["V"] + nx = gari_structure["nx_virt"] + nz = gari_structure["nz_virt"] + + def cost(p): + p = np.clip(p, 1e-15, 1 - 1e-15) + return np.log((1-p)/p) + + def prob(c): + c = np.clip(c, -30, 30) + return 1 / (1 + np.exp(c)) + + P_ez_prime_agg = P_ez + U @ P_ey + P_ex_prime_agg = P_ex + V @ P_ey + + if method == "modeA": + return np.concatenate([P_ez, P_ex, P_ey, P_ez_prime_agg, P_ex_prime_agg]) + + if method == "modeB": + return np.concatenate([ + P_ez, P_ex, P_ey, + np.full_like(P_ez, 0.499), + np.full_like(P_ex, 0.499) + ]) + + lam = 0.001 + C_z_prime_agg = cost(P_ez_prime_agg) + C_x_prime_agg = cost(P_ex_prime_agg) + + if method == "modeC": + return np.concatenate([ + P_ez, P_ex, P_ey, + prob(lam * C_z_prime_agg), + prob(lam * C_x_prime_agg) + ]) + + if method == "modeD": + return np.concatenate([ + np.full_like(P_ez, 0.499), + np.full_like(P_ex, 0.499), + np.full_like(P_ey, 0.499), + P_ez, P_ex + ]) + + C_z = cost(P_ez) + C_x = cost(P_ex) + C_y = cost(P_ey) + + if method == "modeE": + return np.concatenate([ + prob(lam * C_z), + prob(lam * C_x), + prob(lam * C_y), + P_ez, P_ex + ]) + + if method == "modeF": + return np.concatenate([ + np.full_like(P_ez, 0.499), + np.full_like(P_ex, 0.499), + np.full_like(P_ey, 0.499), + P_ez_prime_agg, P_ex_prime_agg + ]) + + if method == "modeG": + return np.concatenate([ + prob(lam * C_z), + prob(lam * C_x), + prob(lam * C_y), + P_ez_prime_agg, P_ex_prime_agg + ]) + + if method == "modeH": + return np.concatenate([ + P_ez, P_ex, P_ey, + P_ez, P_ex + ]) + + U_weighted = U.multiply(P_ey) + max_ey_to_z = U_weighted.max(axis=1).toarray().flatten() + P_ez_prime_max = np.maximum(P_ez, max_ey_to_z) + + V_weighted = V.multiply(P_ey) + max_ey_to_x = V_weighted.max(axis=1).toarray().flatten() + P_ex_prime_max = np.maximum(P_ex, max_ey_to_x) + + if method == "modeI": + return np.concatenate([ + P_ez, P_ex, P_ey, + P_ez_prime_max, P_ex_prime_max + ]) + + if method == "modeJ": + return np.concatenate([ + np.full_like(P_ez, 0.499), + np.full_like(P_ex, 0.499), + np.full_like(P_ey, 0.499), + P_ez_prime_max, P_ex_prime_max + ]) + + if method == "modeK": + return np.concatenate([ + prob(lam * C_z), + prob(lam * C_x), + prob(lam * C_y), + P_ez_prime_max, P_ex_prime_max + ]) + + if method == "modeL": + return np.concatenate([ + P_ez, P_ex, np.full_like(P_ey, 0.499), + P_ez_prime_agg, P_ex_prime_agg + ]) + + if method == "modeM": + return np.concatenate([ + P_ez, P_ex, np.full_like(P_ey, 0.499), + P_ez_prime_max, P_ex_prime_max + ]) + + eps = 1e-15 + log_1m2_P_ez = np.log(np.clip(1 - 2 * P_ez, eps, 1.0)) + log_1m2_P_ex = np.log(np.clip(1 - 2 * P_ex, eps, 1.0)) + log_1m2_P_ey = np.log(np.clip(1 - 2 * P_ey, eps, 1.0)) + + P_ez_prime_xor = 0.5 * (1 - np.exp(log_1m2_P_ez + U @ log_1m2_P_ey)) + P_ex_prime_xor = 0.5 * (1 - np.exp(log_1m2_P_ex + V @ log_1m2_P_ey)) + + if method == "modeN": + return np.concatenate([ + P_ez, P_ex, P_ey, + P_ez_prime_xor, P_ex_prime_xor + ]) + + if method == "modeO": + return np.concatenate([ + prob(lam * C_z), prob(lam * C_x), prob(lam * C_y), + P_ez_prime_xor, P_ex_prime_xor + ]) + + from scipy.sparse import bmat, eye as speye, vstack, hstack + from scipy.optimize import linprog + + P_real_orig = np.concatenate([P_dx, P_dz, P_y]) + b_ub = cost(P_real_orig) + + N_real = len(b_ub) + N_virtual = nx + nz + + I_nx = speye(nx, format='csc') + I_nz = speye(nz, format='csc') + + A_ub_orig = bmat([ + [I_nx, None], + [None, I_nz], + [U.T, V.T] + ], format='csc') + + col_ones_real = np.ones((N_real, 1)) + A_upper = hstack([A_ub_orig, col_ones_real], format='csc') + + I_virt = speye(N_virtual, format='csc') + col_ones_virt = np.ones((N_virtual, 1)) + A_lower = hstack([-I_virt, col_ones_virt], format='csc') + + A_ub_new = vstack([A_upper, A_lower], format='csc') + b_ub_new = np.concatenate([b_ub, np.zeros(N_virtual)]) + bounds = [(0, None)] * (N_virtual + 1) + lambda_reg = 1e-4 + + if method == "modeQ": + c_obj_Q = np.zeros(N_virtual + 1) + c_obj_Q[-1] = -1.0 + res_lp_Q = linprog(c_obj_Q, A_ub=A_ub_new, b_ub=b_ub_new, bounds=bounds, method='highs') + if res_lp_Q.success: + g_virt_Q = res_lp_Q.x[:-1] + g_real_Q = b_ub - A_ub_orig.dot(g_virt_Q) + return np.concatenate([prob(g_real_Q), prob(g_virt_Q)]) + raise NotImplementedError("LP modeQ failed") + + if method == "modeP": + c_obj_P = np.zeros(N_real + N_virtual) + c_obj_P[-N_virtual:] = -1.0 + A_eq_P = hstack([speye(N_real, format='csc'), A_ub_orig], format='csc') + b_eq_P = b_ub + w_z, w_x, w_y = C_z, C_x, C_y + min_w_vals = [np.min(w) for w in (w_z, w_x, w_y) if len(w) > 0] + min_w = min(min_w_vals) if min_w_vals else 1e-10 + dynamic_eps = min_w / 1e5 + bounds_P = [(dynamic_eps, None)] * (N_real + N_virtual) + res_lp_P = linprog(c_obj_P, A_eq=A_eq_P, b_eq=b_eq_P, bounds=bounds_P, method='highs') + if res_lp_P.success: + return np.concatenate([prob(res_lp_P.x[:N_real]), prob(res_lp_P.x[N_real:])]) + return np.concatenate([P_ez, P_ex, P_ey, P_ez_prime_agg, P_ex_prime_agg]) + + if method == "modeR": + c_obj_R = -lambda_reg * np.ones(N_virtual + 1) + c_obj_R[-1] = -1.0 + res_lp_R = linprog(c_obj_R, A_ub=A_ub_new, b_ub=b_ub_new, bounds=bounds, method='highs') + if res_lp_R.success: + g_virt_R = res_lp_R.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_R)), prob(g_virt_R)]) + raise NotImplementedError("LP modeR failed") + + if method == "modeS": + weights_virtual = b_ub[:N_virtual] + lambda_array = lambda_reg * (weights_virtual / np.max(weights_virtual)) + c_obj_S = np.zeros(N_virtual + 1) + c_obj_S[:-1] = -lambda_array + c_obj_S[-1] = -1.0 + res_lp_S = linprog(c_obj_S, A_ub=A_ub_new, b_ub=b_ub_new, bounds=bounds, method='highs') + if res_lp_S.success: + g_virt_S = res_lp_S.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_S)), prob(g_virt_S)]) + raise NotImplementedError("LP modeS failed") + + if method == "modeS2": + W_z = b_ub[:nx] + W_x = b_ub[nx:N_virtual] + W_y = b_ub[N_virtual:] + S_z = W_z + U.dot(W_y) + S_x = W_x + V.dot(W_y) + S_virt = np.concatenate([S_z, S_x]) + lambda_array_S2 = lambda_reg * (S_virt / np.max(S_virt)) + c_obj_S2 = np.zeros(N_virtual + 1) + c_obj_S2[:-1] = -lambda_array_S2 + c_obj_S2[-1] = -1.0 + res_lp_S2 = linprog(c_obj_S2, A_ub=A_ub_new, b_ub=b_ub_new, bounds=bounds, method='highs') + if res_lp_S2.success: + g_virt_S2 = res_lp_S2.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_S2)), prob(g_virt_S2)]) + raise NotImplementedError("LP modeS2 failed") + + if method == "modeSO": + weights_virtual = b_ub[:N_virtual] + normalized_W = weights_virtual / np.max(weights_virtual) + safe_lambda_reg_SO = 0.99 / np.sum(normalized_W) + lambda_array_SO = safe_lambda_reg_SO * normalized_W + c_obj_SO = np.zeros(N_virtual + 1) + c_obj_SO[:-1] = -lambda_array_SO + c_obj_SO[-1] = -1.0 + res_lp_SO = linprog(c_obj_SO, A_ub=A_ub_new, b_ub=b_ub_new, bounds=bounds, method='highs') + if res_lp_SO.success: + g_virt_SO = res_lp_SO.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_SO)), prob(g_virt_SO)]) + raise NotImplementedError("LP modeSO failed") + + if method == "modeSO2": + W_z = b_ub[:nx] + W_x = b_ub[nx:N_virtual] + W_y = b_ub[N_virtual:] + S_z_o = W_z + U.dot(W_y) + S_x_o = W_x + V.dot(W_y) + S_virt_o = np.concatenate([S_z_o, S_x_o]) + normalized_S = S_virt_o / np.max(S_virt_o) + safe_lambda_reg_SO2 = 0.99 / np.sum(normalized_S) + lambda_array_SO2 = safe_lambda_reg_SO2 * normalized_S + c_obj_SO2 = np.zeros(N_virtual + 1) + c_obj_SO2[:-1] = -lambda_array_SO2 + c_obj_SO2[-1] = -1.0 + res_lp_SO2 = linprog(c_obj_SO2, A_ub=A_ub_new, b_ub=b_ub_new, bounds=bounds, method='highs') + if res_lp_SO2.success: + g_virt_SO2 = res_lp_SO2.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_SO2)), prob(g_virt_SO2)]) + raise NotImplementedError("LP modeSO2 failed") + + # Modes U and V + w_z, w_x, w_y = C_z, C_x, C_y + min_w_vals = [np.min(w) for w in (w_z, w_x, w_y) if len(w) > 0] + min_w = min(min_w_vals) if min_w_vals else 1e-10 + dynamic_eps = min_w / 1e5 + + A_upper_UV = hstack([A_ub_orig, np.zeros((N_real, 1))], format='csc') + A_ub_UV = vstack([A_upper_UV, A_lower], format='csc') + bounds_UV = [(0, None)] * N_virtual + [(dynamic_eps, None)] + b_ub_UV = np.concatenate([b_ub - dynamic_eps, np.zeros(N_virtual)]) + + if method == "modeU": + c_obj_U = np.zeros(N_virtual + 1) + c_obj_U[:-1] = -1.0 + c_obj_U[-1] = -1.0 + res_lp_U = linprog(c_obj_U, A_ub=A_ub_UV, b_ub=b_ub_UV, bounds=bounds_UV, method='highs') + if res_lp_U.success: + g_virt_U = res_lp_U.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_U)), prob(g_virt_U)]) + raise NotImplementedError("LP modeU failed") + + if method == "modeV": + W_z = b_ub[:nx] + W_x = b_ub[nx:N_virtual] + W_y = b_ub[N_virtual:] + S_z = W_z + U.dot(W_y) + S_x = W_x + V.dot(W_y) + S_virt = np.concatenate([S_z, S_x]) + lambda_array_V = lambda_reg * (S_virt / np.max(S_virt)) + c_obj_V = np.zeros(N_virtual + 1) + c_obj_V[:-1] = -lambda_array_V + c_obj_V[-1] = -1.0 + res_lp_V = linprog(c_obj_V, A_ub=A_ub_UV, b_ub=b_ub_UV, bounds=bounds_UV, method='highs') + if res_lp_V.success: + g_virt_V = res_lp_V.x[:-1] + return np.concatenate([prob(b_ub - A_ub_orig.dot(g_virt_V)), prob(g_virt_V)]) + raise NotImplementedError("LP modeV failed") + + raise ValueError(f"Unknown method {method}") + + +def get_detector_orderings(gari_structure: dict, det_types: np.ndarray, ordering_name: str) -> list: + hx_csc = gari_structure["hx_csc"] + hz_csc = gari_structure["hz_csc"] + U_csr = gari_structure["U"].tocsr() + V_csr = gari_structure["V"].tocsr() + + i_hx_only = gari_structure["i_hx_only"] + i_hz_only = gari_structure["i_hz_only"] + i_hy = gari_structure["i_hy"] + + x_orig_indices = gari_structure["x_orig_indices"] + z_orig_indices = gari_structure["z_orig_indices"] + + nx_virt = gari_structure["nx_virt"] + nz_virt = gari_structure["nz_virt"] + nx_real = gari_structure["nx_real"] + nz_real = gari_structure["nz_real"] + num_original_detectors = gari_structure["num_original_detectors"] + num_gari_detectors = gari_structure["num_gari_detectors"] + + # Precompute time_vx_old, time_vz_old, time_vx_new, time_vz_new, time_vx_min_old, time_vz_min_old, time_vx_min_new, time_vz_min_new + + time_vy = [max( + np.max([x_orig_indices[r] for r in hx_csc.indices[hx_csc.indptr[c]:hx_csc.indptr[c+1]]]) if hx_csc.indptr[c+1] > hx_csc.indptr[c] else 0, + np.max([z_orig_indices[r] for r in hz_csc.indices[hz_csc.indptr[c]:hz_csc.indptr[c+1]]]) if hz_csc.indptr[c+1] > hz_csc.indptr[c] else 0 + ) for c in i_hy] + + time_vx_old = [np.max([x_orig_indices[r] for r in hx_csc.indices[hx_csc.indptr[c]:hx_csc.indptr[c+1]]]) if hx_csc.indptr[c+1] > hx_csc.indptr[c] else 0 for c in i_hx_only] + time_vz_old = [np.max([z_orig_indices[r] for r in hz_csc.indices[hz_csc.indptr[c]:hz_csc.indptr[c+1]]]) if hz_csc.indptr[c+1] > hz_csc.indptr[c] else 0 for c in i_hz_only] + + time_vx_new = [max( + time_vx_old[i], + max([time_vy[k] for k in U_csr.indices[U_csr.indptr[i]:U_csr.indptr[i+1]]]) if U_csr.indptr[i+1] > U_csr.indptr[i] else 0 + ) for i, c in enumerate(i_hx_only)] + + time_vz_new = [max( + time_vz_old[i], + max([time_vy[k] for k in V_csr.indices[V_csr.indptr[i]:V_csr.indptr[i+1]]]) if V_csr.indptr[i+1] > V_csr.indptr[i] else 0 + ) for i, c in enumerate(i_hz_only)] + + time_vx_min_old = [np.min([x_orig_indices[r] for r in hx_csc.indices[hx_csc.indptr[c]:hx_csc.indptr[c+1]]]) if hx_csc.indptr[c+1] > hx_csc.indptr[c] else 0 for c in i_hx_only] + time_vz_min_old = [np.min([z_orig_indices[r] for r in hz_csc.indices[hz_csc.indptr[c]:hz_csc.indptr[c+1]]]) if hz_csc.indptr[c+1] > hz_csc.indptr[c] else 0 for c in i_hz_only] + + time_vy_min = [] + for c in i_hy: + t_hx = np.min([x_orig_indices[r] for r in hx_csc.indices[hx_csc.indptr[c]:hx_csc.indptr[c+1]]]) if hx_csc.indptr[c+1] > hx_csc.indptr[c] else float('inf') + t_hz = np.min([z_orig_indices[r] for r in hz_csc.indices[hz_csc.indptr[c]:hz_csc.indptr[c+1]]]) if hz_csc.indptr[c+1] > hz_csc.indptr[c] else float('inf') + t_min = min(t_hx, t_hz) + time_vy_min.append(t_min if t_min != float('inf') else 0) + + time_vx_min_new = [] + for i, c in enumerate(i_hx_only): + t_y = [time_vy_min[k] for k in U_csr.indices[U_csr.indptr[i]:U_csr.indptr[i+1]]] + min_y = min(t_y) if t_y else float('inf') + t_min = min(time_vx_min_old[i], min_y) + time_vx_min_new.append(t_min if t_min != float('inf') else 0) + + time_vz_min_new = [] + for i, c in enumerate(i_hz_only): + t_y = [time_vy_min[k] for k in V_csr.indices[V_csr.indptr[i]:V_csr.indptr[i+1]]] + min_y = min(t_y) if t_y else float('inf') + t_min = min(time_vz_min_old[i], min_y) + time_vz_min_new.append(t_min if t_min != float('inf') else 0) + + # Common mappings + real_x_gari = list(range(0, nx_real)) + real_z_gari = list(range(nx_real, nx_real + nz_real)) + virt_x_gari = list(range(nx_real + nz_real, nx_real + nz_real + nx_virt)) + virt_z_gari = list(range(nx_real + nz_real + nx_virt, num_gari_detectors)) + + x_map = {orig_idx: int(gari_idx) for gari_idx, orig_idx in enumerate(x_orig_indices)} + z_map = {orig_idx: int(nx_real + gari_idx) for gari_idx, orig_idx in enumerate(z_orig_indices)} + + if ordering_name == "order_1" or ordering_name == "order1": + return real_x_gari + virt_x_gari + real_z_gari + virt_z_gari + elif ordering_name == "order_2" or ordering_name == "order2": + return real_x_gari + real_z_gari + virt_x_gari + virt_z_gari + elif ordering_name == "order_3" or ordering_name == "order3": + return real_x_gari + virt_z_gari + real_z_gari + virt_x_gari + elif ordering_name == "order_4" or ordering_name == "order4": + return real_x_gari + real_z_gari + virt_z_gari + virt_x_gari + elif ordering_name == "order_8" or ordering_name == "order8": + return real_x_gari + virt_x_gari + virt_z_gari + real_z_gari + elif ordering_name == "order_9" or ordering_name == "order9": + return real_z_gari + real_x_gari + virt_z_gari + virt_x_gari + elif ordering_name == "order_10" or ordering_name == "order10": + return real_z_gari + real_x_gari + virt_x_gari + virt_z_gari + + real_gari_chronological = [] + for orig_idx in range(num_original_detectors): + if orig_idx in x_map: + real_gari_chronological.append(x_map[orig_idx]) + elif orig_idx in z_map: + real_gari_chronological.append(z_map[orig_idx]) + + if ordering_name == "order_7" or ordering_name == "order7": + virt_with_time_old = [] + for c in range(nx_virt): + virt_with_time_old.append((time_vx_old[c], nx_real + nz_real + c)) + for c in range(len(time_vz_old)): + virt_with_time_old.append((time_vz_old[c], nx_real + nz_real + nx_virt + c)) + virt_with_time_old.sort(key=lambda x: x[0]) + return real_gari_chronological + [v[1] for v in virt_with_time_old] + + elif ordering_name == "order_7a" or ordering_name == "order7a": + virt_with_time_new = [] + for c in range(nx_virt): + virt_with_time_new.append((time_vx_new[c], nx_real + nz_real + c)) + for c in range(len(time_vz_new)): + virt_with_time_new.append((time_vz_new[c], nx_real + nz_real + nx_virt + c)) + virt_with_time_new.sort(key=lambda x: x[0]) + return real_gari_chronological + [v[1] for v in virt_with_time_new] + + elif ordering_name == "order_7b" or ordering_name == "order7b": + virt_with_time_min_old = [] + for c in range(nx_virt): + virt_with_time_min_old.append((time_vx_min_old[c], nx_real + nz_real + c)) + for c in range(len(time_vz_min_old)): + virt_with_time_min_old.append((time_vz_min_old[c], nx_real + nz_real + nx_virt + c)) + virt_with_time_min_old.sort(key=lambda x: x[0]) + return real_gari_chronological + [v[1] for v in virt_with_time_min_old] + + elif ordering_name == "order_7c" or ordering_name == "order7c": + virt_with_time_min_new = [] + for c in range(nx_virt): + virt_with_time_min_new.append((time_vx_min_new[c], nx_real + nz_real + c)) + for c in range(len(time_vz_min_new)): + virt_with_time_min_new.append((time_vz_min_new[c], nx_real + nz_real + nx_virt + c)) + virt_with_time_min_new.sort(key=lambda x: x[0]) + return real_gari_chronological + [v[1] for v in virt_with_time_min_new] + + def get_order_11_variant(v_x_times, v_z_times): + all_det = [] + for orig_idx in range(num_original_detectors): + if orig_idx in x_map: + all_det.append((orig_idx, x_map[orig_idx])) + elif orig_idx in z_map: + all_det.append((orig_idx, z_map[orig_idx])) + for c in range(nx_virt): + all_det.append((v_x_times[c] + 0.1, nx_real + nz_real + c)) + for c in range(len(v_z_times)): + all_det.append((v_z_times[c] + 0.1, nx_real + nz_real + nx_virt + c)) + all_det.sort(key=lambda x: x[0]) + return [v[1] for v in all_det] + + if ordering_name == "order_11" or ordering_name == "order11": + return get_order_11_variant(time_vx_new, time_vz_new) + elif ordering_name == "order_11a" or ordering_name == "order11a": + return get_order_11_variant(time_vx_old, time_vz_old) + elif ordering_name == "order_11b" or ordering_name == "order11b": + return get_order_11_variant(time_vx_min_new, time_vz_min_new) + + def get_order_12_variant(v_x_times, v_z_times): + all_det = [] + for orig_idx in range(num_original_detectors): + if orig_idx in x_map: + all_det.append((orig_idx, x_map[orig_idx])) + elif orig_idx in z_map: + all_det.append((orig_idx, z_map[orig_idx])) + for c in range(nx_virt): + all_det.append((v_x_times[c] - 0.1, nx_real + nz_real + c)) + for c in range(len(v_z_times)): + all_det.append((v_z_times[c] - 0.1, nx_real + nz_real + nx_virt + c)) + all_det.sort(key=lambda x: x[0], reverse=True) + return [v[1] for v in all_det] + + if ordering_name == "order_12" or ordering_name == "order12": + return get_order_12_variant(time_vx_min_new, time_vz_min_new) + elif ordering_name == "order_12a" or ordering_name == "order12a": + return get_order_12_variant(time_vx_min_old, time_vz_min_old) + elif ordering_name == "order_12b" or ordering_name == "order12b": + return get_order_12_variant(time_vx_new, time_vz_new) + + raise ValueError(f"Unknown ordering {ordering_name}") + + +def process_directory(input_path): + input_path = get_target_path(input_path) + if os.path.isdir(input_path): + stim_files = glob.glob(os.path.join(input_path, "**", "*.stim"), recursive=True) + else: + if "*" not in input_path: + input_path += "*" + stim_files = [f for f in glob.glob(input_path, recursive=True) if f.endswith(".stim")] + for stim_path in stim_files: + print(f"Processing {stim_path}") + try: + if "color_code" in stim_path or "colorcodes" in stim_path: + if "superdense_color_code" not in os.path.basename(stim_path): + continue + + circuit = stim.Circuit.from_file(stim_path) + dem = circuit.detector_error_model( + decompose_errors=False, + flatten_loops=True, + ignore_decomposition_failures=True + ) + det_types = get_detector_types(circuit) + H, L, priors, errors = dem_to_check_matrices(dem) + + gari_structure = gari_transform(H, L, det_types) + + gari_matrix = gari_structure["gari_matrix"] + gari_obs_matrix = gari_structure["gari_obs_matrix"] + gari_obs_matrix_og = gari_structure["gari_obs_matrix_og"] + + stim_dir = os.path.dirname(stim_path) + stim_stem = os.path.splitext(os.path.basename(stim_path))[0] + gari_dir = os.path.join(stim_dir, "gari") + os.makedirs(gari_dir, exist_ok=True) + base_path = os.path.join(gari_dir, stim_stem) + + # modes_to_generate = [ + # "modeA", "modeB", "modeC", "modeF", "modeG", "modeH", + # "modeI", "modeJ", "modeK", "modeL", "modeM", "modeN", + # "modeO", "modeP", "modeQ", "modeR", "modeS", "modeU", + # "modeV", "modeS2", "modeSO", "modeSO2" + # ] + modes_to_generate = [ + "modeA", "modeN", + "modeO", "modeQ", "modeR", "modeS", "modeS2", "modeSO", "modeSO2" + ] + + for mode_name in modes_to_generate: + try: + priors_for_mode = assign_prior_weights(gari_structure, mode_name, priors) + dem_for_mode = matrices_to_dem(gari_matrix, gari_obs_matrix, priors_for_mode) + dem_for_mode.to_file(base_path + f"_{mode_name}.dem") + + dem_for_mode_og = matrices_to_dem(gari_matrix, gari_obs_matrix_og, priors_for_mode) + dem_for_mode_og.to_file(base_path + f"_ogL_{mode_name}.dem") + except NotImplementedError as e: + print(f"Skipping {mode_name}: {e}") + + x_orig_indices = np.where(det_types == 1)[0] + z_orig_indices = np.where(det_types == 3)[0] + + nx = len(x_orig_indices) + + mapping = [-1] * dem.num_detectors + for gari_idx, orig_idx in enumerate(x_orig_indices): + mapping[int(orig_idx)] = int(gari_idx) + for gari_idx, orig_idx in enumerate(z_orig_indices): + mapping[int(orig_idx)] = int(nx + gari_idx) + + orderings = { + "order2": get_detector_orderings(gari_structure, det_types, "order2"), + "order4": get_detector_orderings(gari_structure, det_types, "order4"), + "order7": get_detector_orderings(gari_structure, det_types, "order7"), + "order9": get_detector_orderings(gari_structure, det_types, "order9"), + "order10": get_detector_orderings(gari_structure, det_types, "order10"), + } + + det_orders_clean = {} + for k, v in orderings.items(): + det_orders_clean[k] = [int(x) for x in v] + + mapping_dict = { + "num_original_detectors": dem.num_detectors, + "mapping": mapping, + "det_orders": det_orders_clean + } + with open(base_path + "_mapping.json", "w") as f: + json.dump(mapping_dict, f, indent=2) + + print(f"Successfully processed {stim_path}") + except Exception as e: + print(f"Error processing {stim_path}: {e}") + +if __name__ == "__main__": + import argparse + parser = argparse.ArgumentParser(description="Process stim files and generate Gari DEMs.") + parser.add_argument("input_path", type=str, help="Input directory, file prefix, or glob pattern for .stim files") + args = parser.parse_args() + process_directory(args.input_path) diff --git a/src/py/gari_simulation_test.py b/src/py/gari_simulation_test.py new file mode 100644 index 0000000..da656d4 --- /dev/null +++ b/src/py/gari_simulation_test.py @@ -0,0 +1,192 @@ +import stim +import numpy as np +import matplotlib.pyplot as plt +from scipy.sparse import csc_matrix +import glob +import os +from _tesseract_py_util.gari_dem_utils import get_detector_types, dem_to_check_matrices, matrices_to_dem, gari_transform, assign_prior_weights, get_detector_orderings + +def get_target_path(relative_path): + workspace_root = os.environ.get("BUILD_WORKSPACE_DIRECTORY", "") + return os.path.join(workspace_root, relative_path) + +def testdata_one_basis_circuit(circuit_og,z_basis=True): + circuit = circuit_og.flattened().copy() + for i in range(len(circuit)-1,-1,-1): + if circuit[i].name == "DETECTOR": + args = circuit[i].gate_args_copy() + if z_basis and len(args) > 3 and args[3] <= 2: + # remove x detectors + circuit.pop(i) + + if not z_basis and len(args) > 3 and args[3] >= 3: + circuit.pop(i) + + return circuit + +def load_circuit(codename,d,r,p,obs_basis,noise='si1000'): + fname="" + if codename == "surfacecodes": + fname = f"testdata/{codename}/r={r},d={d},p={p},noise={noise},c=surface_code_{obs_basis}" + if codename == "bivariatebicyclecodes": + fname = f"testdata/{codename}/r={r},d={d},p={p},noise={noise},c=bivariate_bicycle_{obs_basis}" + if codename == "colorcodes": + fname = f"testdata/{codename}/r={r},d={d},p={p},noise={noise},c=superdense_color_code_{obs_basis}" + circuit = None + try: + fname = get_target_path(fname) + fname = glob.glob(fname + '*.stim')[0] + circuit = stim.Circuit.from_file(fname) + print("successful") + except: + raise("could not find the circuit") + + return circuit + + + +def test_gari_transform(): + print("Reading Circuit...") + import sys + if len(sys.argv) > 1: + circuit = stim.Circuit.from_file(sys.argv[1]) + else: + circuit = load_circuit("surfacecodes", d=3, r=3, p=0.001, obs_basis='Z') + # circuit = load_circuit("bivariatebicyclecodes", d=10, r=10, p=0.001, obs_basis='Z') + # circuit = load_circuit("colorcodes", d=7, r=7, p=0.001, obs_basis='Z') + + dem = circuit.detector_error_model() + print("Extracting DEM...") + dem = circuit.detector_error_model( + decompose_errors=False, + flatten_loops=True, + ignore_decomposition_failures=True + ) + # dem2 = circuit.detector_error_model( + # decompose_errors=True, + # flatten_loops=True, + # ignore_decomposition_failures=True + # ) + # print(dem) + # print("decomposed",dem2) + det_types = get_detector_types(circuit) + print(f"Total Detectors: {dem.num_detectors}") + print(f"X Detectors: {np.sum(det_types == 1)}, Z Detectors: {np.sum(det_types == 3)}") + + print("Building original Check Matrix...") + H, L, priors, errors = dem_to_check_matrices(dem) + # H, L, priors, errors = dem_to_check_matrices(dem2) + print(f"Original Check Matrix Shape: {H.shape}") + print(f"Original Observables Matrix Shape: {L.shape}") + + # print("Applying Gari Transform (matrices)...") + # H_gari, L_gari, priors_gari, L_ez_prime, L_ex_prime = gari_transform(H, L, det_types, priors, return_dem=False) + # print(f"Gari Matrix Shape: {H_gari.shape}") + # print(f"Gari Observables Matrix Shape: {L_gari.shape}") + # print(f"Gari Priors Shape: {priors_gari.shape}") + + + print("Applying Gari Transform (return_dem=False)...") + gari_structure = gari_transform(H, L, det_types) + gari_matrix = gari_structure["gari_matrix"] + gari_obs_matrix = gari_structure["gari_obs_matrix"] + gari_obs_matrix_og = gari_structure["gari_obs_matrix_og"] + + p_agg = assign_prior_weights(gari_structure, "modeA", priors) + dem_agg = matrices_to_dem(gari_matrix, gari_obs_matrix, p_agg) + + modes_to_test = { + "Mode N: ez,ex,ey Keep, ez',ex' XOR Aggregated": "modeN", + "Mode P: LP based Weight Distribution (Fixed Epsilon)": "modeP", + "Mode Q: LP based Weight Distribution (No Lambda)": "modeQ", + "Mode R: LP based Weight Distribution (Uniform Lambda)": "modeR", + "Mode S: LP based Weight Distribution (Weighted Lambda)": "modeS", + "Mode U: Max-Min Formulation (Zero-Cost Reals allowed, Lambda = 1)": "modeU", + "Mode V: Max-Min Formulation (Zero-Cost Reals allowed, Topologically Weighted Lambda)": "modeV", + "Mode S2: Max-Min Formulation (g_real >= t, Topologically Weighted Lambda)": "modeS2", + "Mode SO: Max-Min Formulation (g_real >= t, Original Weighted Lambda, Safe Lambda Reg)": "modeSO", + "Mode SO2: Max-Min Formulation (g_real >= t, Topologically Weighted Lambda, Safe Lambda Reg)": "modeSO2" + } + + prior_modes = {} + for pretty_name, mode_id in modes_to_test.items(): + try: + p_weights = assign_prior_weights(gari_structure, mode_id, priors) + prior_modes[pretty_name] = matrices_to_dem(gari_matrix, gari_obs_matrix, p_weights) + except NotImplementedError: + pass + + print(f"Gari Matrix Shape: {gari_matrix.shape}") + print(f"Gari DEM has {dem_agg.num_detectors} detectors and {dem_agg.num_errors} errors.") + orig_row_weights = np.sum(H.toarray(), axis=1) + print(f"Original Avg Row Weight: {np.mean(orig_row_weights):.2f}") + + print("Success! Gari matrix generated and verified.") + + import time + from tesseract_decoder.tesseract_sinter_compat import make_tesseract_sinter_decoders_dict + import tesseract_decoder.tesseract as tesseract + import tesseract_decoder.utils + has_tesseract = True + + if has_tesseract: + print("\n--- Running Tesseract Decoder Comparison ---") + + num_shots = 100 + sampler = circuit.compile_detector_sampler(seed=0) + shots, obs_shots = sampler.sample(shots=num_shots, separate_observables=True) + + sinter_decoders = make_tesseract_sinter_decoders_dict() + short_beam_decoder_obj = sinter_decoders["tesseract-long-beam"] + base_config = short_beam_decoder_obj.compile_decoder_for_dem(dem=dem).decoder.config + base_config.det_orders = [list(range(dem.num_detectors))] + + print("\nCompiling Tesseract for Original DEM (Forced 1 Order)...") + decoder_orig = tesseract.TesseractDecoder(base_config) + + print("Decoding Original DEM (1 Order)...") + start_time = time.time() + predicted_obs_orig = decoder_orig.decode_batch(shots) + time_orig = time.time() - start_time + + correct_orig = np.sum(np.all(predicted_obs_orig == obs_shots, axis=1)) + print(f"Original DEM (1 Order): {correct_orig}/{num_shots} correct, Time: {time_orig:.4f}s") + + orders = { + "order2": get_detector_orderings(gari_structure, det_types, "order2"), + "order4": get_detector_orderings(gari_structure, det_types, "order4"), + "order7": get_detector_orderings(gari_structure, det_types, "order7"), + "order9": get_detector_orderings(gari_structure, det_types, "order9"), + "order10": get_detector_orderings(gari_structure, det_types, "order10"), + } + + is_x_det = (det_types == 1) + is_z_det = (det_types == 3) + num_virtual = dem_agg.num_detectors - dem.num_detectors + x_shots = shots[:, is_x_det] + z_shots = shots[:, is_z_det] + virtual_shots = np.zeros((num_shots, num_virtual), dtype=bool) + gari_shots = np.concatenate([x_shots, z_shots, virtual_shots], axis=1) + + print("\nTesting Priority Modes on Gari DEM ...") + + for mode_name, target_dem in prior_modes.items(): + print(f"\n>> {mode_name}") + for name, order in orders.items(): + base_gari_config = short_beam_decoder_obj.compile_decoder_for_dem(dem=target_dem).decoder.config + base_gari_config.det_orders = [order] + base_gari_config.no_revisit_dets = False + + decoder_gari = tesseract.TesseractDecoder(base_gari_config) + + start_time = time.time() + predicted_obs_gari = decoder_gari.decode_batch(gari_shots) + time_gari = time.time() - start_time + correct_gari = np.sum(np.all(predicted_obs_gari == obs_shots, axis=1)) + print(f"{name[:30]:30s} | {correct_gari}/{num_shots} correct | {time_gari:.4f}s") + +if __name__ == "__main__": + # plot_gari_structure() + test_gari_transform() + + diff --git a/src/simplex_main.cc b/src/simplex_main.cc index 95437ae..2937829 100644 --- a/src/simplex_main.cc +++ b/src/simplex_main.cc @@ -28,6 +28,7 @@ struct Args { std::string circuit_path; std::string dem_path; bool no_merge_errors = false; + std::string det_mapping_file; // Sampling options size_t sample_num_shots = 0; @@ -138,6 +139,18 @@ struct Args { void extract(SimplexConfig& config, std::vector& shots, std::unique_ptr& writer) { + std::vector det_mapping; + uint64_t num_original_detectors = 0; + if (!det_mapping_file.empty()) { + std::ifstream f(det_mapping_file); + if (!f) { + throw std::invalid_argument("Could not open the mapping file: " + det_mapping_file); + } + nlohmann::json j = nlohmann::json::parse(f); + num_original_detectors = j.at("num_original_detectors").get(); + det_mapping = j.at("mapping").get>(); + } + // Get a circuit, if available stim::Circuit circuit; if (!circuit_path.empty()) { @@ -182,7 +195,7 @@ struct Args { shots[k].obs_mask = obs_T[k]; for (size_t d = 0; d < num_detectors; d++) { if (dets[d][k]) { - shots[k].hits.push_back(d); + shots[k].hits.push_back(!det_mapping.empty() ? det_mapping[d] : d); } } } @@ -195,14 +208,20 @@ struct Args { throw std::invalid_argument("Could not open the file: " + in_fname); } stim::FileFormatData shots_in_format = stim::format_name_to_enum_map().at(in_format); + size_t num_dets = det_mapping.empty() ? config.dem.count_detectors() : num_original_detectors; auto reader = stim::MeasureRecordReader::make( - shots_file, shots_in_format.id, 0, config.dem.count_detectors(), + shots_file, shots_in_format.id, 0, num_dets, append_observables * config.dem.count_observables()); // Load the shots from a file stim::SparseShot sparse_shot; sparse_shot.clear(); while (reader->start_and_read_entire_record(sparse_shot)) { + if (!det_mapping.empty()) { + for (auto& hit : sparse_shot.hits) { + hit = det_mapping[hit]; + } + } shots.push_back(sparse_shot); sparse_shot.clear(); } @@ -277,6 +296,7 @@ int main(int argc, char* argv[]) { program.add_argument("--no-merge-errors") .help("If provided, will not merge identical error mechanisms.") .store_into(args.no_merge_errors); + program.add_argument("--det-mapping-file").help("JSON file containing detector mapping").default_value(std::string("")).store_into(args.det_mapping_file); program.add_argument("--sample-num-shots") .help( "If provided, will sample the requested number of shots from the " diff --git a/src/tesseract_main.cc b/src/tesseract_main.cc index c54f82b..237f838 100644 --- a/src/tesseract_main.cc +++ b/src/tesseract_main.cc @@ -33,6 +33,8 @@ struct Args { std::string circuit_path; std::string dem_path; bool no_merge_errors = false; + std::string det_mapping_file; + std::string custom_order; // Manifold orientation options uint64_t det_order_seed; @@ -179,6 +181,35 @@ struct Args { void extract(TesseractConfig& config, std::vector& shots, std::unique_ptr& writer) { + std::vector det_mapping; + uint64_t num_original_detectors = 0; + std::vector> custom_det_orders; + if (!det_mapping_file.empty()) { + std::ifstream f(det_mapping_file); + if (!f) { + throw std::invalid_argument("Could not open the mapping file: " + det_mapping_file); + } + nlohmann::json j = nlohmann::json::parse(f); + num_original_detectors = j.at("num_original_detectors").get(); + det_mapping = j.at("mapping").get>(); + + if (!custom_order.empty() && j.contains("det_orders")) { + if (custom_order == "all") { + for (auto& [key, value] : j.at("det_orders").items()) { + custom_det_orders.push_back(value.get>()); + } + } else { + if (j.at("det_orders").contains(custom_order)) { + custom_det_orders.push_back(j.at("det_orders").at(custom_order).get>()); + } else { + throw std::invalid_argument("Mapping file does not contain custom order: " + custom_order); + } + } + } + } else if (!custom_order.empty()) { + throw std::invalid_argument("--custom-order requires a --det-mapping-file containing the orders."); + } + // Get a circuit, if available stim::Circuit circuit; if (!circuit_path.empty()) { @@ -225,15 +256,25 @@ struct Args { std::cout << ")" << std::endl; } } - DetOrder order = DetOrder::DetIndex; - if (det_order_bfs) { - order = DetOrder::DetBFS; - } else if (det_order_index) { - order = DetOrder::DetIndex; - } else if (det_order_coordinate) { - order = DetOrder::DetCoordinate; + if (!custom_det_orders.empty()) { + for (const auto& order : custom_det_orders) { + if (order.size() != config.dem.count_detectors()) { + throw std::invalid_argument("Custom detector order size does not match DEM detector count."); + } + } + config.det_orders = custom_det_orders; + num_det_orders = custom_det_orders.size(); + } else { + DetOrder order = DetOrder::DetIndex; + if (det_order_bfs) { + order = DetOrder::DetBFS; + } else if (det_order_index) { + order = DetOrder::DetIndex; + } else if (det_order_coordinate) { + order = DetOrder::DetCoordinate; + } + config.det_orders = build_det_orders(config.dem, num_det_orders, order, det_order_seed); } - config.det_orders = build_det_orders(config.dem, num_det_orders, order, det_order_seed); } if (sample_num_shots > 0) { @@ -248,7 +289,7 @@ struct Args { shots[k].obs_mask = obs_T[k]; for (size_t d = 0; d < num_detectors; d++) { if (dets[d][k]) { - shots[k].hits.push_back(d); + shots[k].hits.push_back(!det_mapping.empty() ? det_mapping[d] : d); } } } @@ -261,14 +302,20 @@ struct Args { throw std::invalid_argument("Could not open the file: " + in_fname); } stim::FileFormatData shots_in_format = stim::format_name_to_enum_map().at(in_format); + size_t num_dets = det_mapping.empty() ? config.dem.count_detectors() : num_original_detectors; auto reader = stim::MeasureRecordReader::make( - shots_file, shots_in_format.id, 0, config.dem.count_detectors(), + shots_file, shots_in_format.id, 0, num_dets, append_observables * config.dem.count_observables()); // Load the shots from a file stim::SparseShot sparse_shot; sparse_shot.clear(); while (reader->start_and_read_entire_record(sparse_shot)) { + if (!det_mapping.empty()) { + for (auto& hit : sparse_shot.hits) { + hit = det_mapping[hit]; + } + } shots.push_back(sparse_shot); sparse_shot.clear(); } @@ -347,6 +394,8 @@ int main(int argc, char* argv[]) { Args args; program.add_argument("--circuit").help("Stim circuit file path").store_into(args.circuit_path); program.add_argument("--dem").help("Stim dem file path").store_into(args.dem_path); + program.add_argument("--det-mapping-file").help("JSON file containing detector mapping").default_value(std::string("")).store_into(args.det_mapping_file); + program.add_argument("--custom-order").help("Specific detector order to use from mapping JSON (e.g. 'order5' or 'all')").default_value(std::string("")).store_into(args.custom_order); program.add_argument("--no-merge-errors") .help("If provided, will not merge identical error mechanisms.") .store_into(args.no_merge_errors); @@ -650,6 +699,8 @@ int main(int argc, char* argv[]) { nlohmann::json stats_json = { {"circuit_path", args.circuit_path}, {"dem_path", args.dem_path}, + {"custom_order", args.custom_order}, + {"det_mapping_file", args.det_mapping_file}, {"max_errors", args.max_errors}, {"sample_seed", args.sample_seed},