% Generate a compact Markdown simulation report.

thisdir = fileparts(mfilename('fullpath'));
rootdir = fullfile(thisdir, '..');
addpath(genpath(fullfile(rootdir, 'src')));

rng(7);
report_path = fullfile(rootdir, 'docs', 'SIMULATION_REPORT.md');
fid = fopen(report_path, 'w');
if fid < 0
    error('generate_simulation_report:open_failed', 'Could not open report file.');
end

cleanup = onCleanup(@() fclose(fid));

fprintf(fid, '# Simulation Report\n\n');
fprintf(fid, 'Generated by `octave --no-gui examples/generate_simulation_report.m`.\n\n');
fprintf(fid, 'This report uses intentionally small sample counts so it can run quickly in CI or a local terminal. Exact recovery tests live in `tests/run_all_tests.m`.\n\n');
fprintf(fid, '## Reproducibility Metadata\n\n');
fprintf(fid, '| Field | Value |\n');
fprintf(fid, '| --- | --- |\n');
report_timestamp = getenv('QEC_REPORT_TIMESTAMP');
if isempty(report_timestamp)
    report_timestamp = 'reproducible';
end
fprintf(fid, '| Generated at | `%s` |\n', report_timestamp);
fprintf(fid, '| Octave version | `%s` |\n', version());
fprintf(fid, '| RNG seed | `7` |\n');
fprintf(fid, '| Surface-3 cache version | `%s` |\n\n', surface3_cache_version());

fprintf(fid, '## Code Summary\n\n');
fprintf(fid, '| Code | Physical qubits | Distance / model | Correctness coverage |\n');
fprintf(fid, '| --- | ---: | --- | --- |\n');
fprintf(fid, '| Bit-flip repetition | 3 | Distance 3 for X errors | Corrects one X error; Monte Carlo X-noise sweep |\n');
fprintf(fid, '| Phase-flip repetition | 3 | Distance 3 for Z errors | Corrects one Z error |\n');
fprintf(fid, '| 5-qubit perfect | 5 | Distance 3 stabilizer code | Corrects every single-qubit X/Y/Z error |\n');
fprintf(fid, '| Steane | 7 | Distance 3 CSS code | Corrects every single-qubit X/Y/Z error |\n');
fprintf(fid, '| Shor | 9 | Distance 3 concatenated-style code | Corrects every single-qubit X/Y/Z error |\n');
fprintf(fid, '| Bacon-Shor | 9 | Distance 3 subsystem Pauli-frame model | Corrects every single-qubit X/Y/Z error by logical residual parity |\n');
fprintf(fid, '| Surface-3 prototype | 9 data qubits | X/Z/Pauli code-capacity model | Minimum-weight lookup corrects single X, Z, and Y errors |\n');
fprintf(fid, '| Generic surface layout | odd d-by-d data grids | Code-capacity Pauli model | Exact d=3 lookup, bounded d=5 lookup, bounded graph baseline, and peeling fallback for decoder comparisons |\n\n');

fprintf(fid, '## Depolarizing Noise Snapshot\n\n');
fprintf(fid, 'Independent per-qubit depolarizing noise; logical failure is estimated by state fidelity after recovery.\n\n');
fprintf(fid, '| Code | p | Trials | Logical failure |\n');
fprintf(fid, '| --- | ---: | ---: | ---: |\n');
codes = {'bitflip', 'phaseflip', 'fivequbit', 'steane'};
labels = {'Bit-flip', 'Phase-flip', '5-qubit', 'Steane'};
ps = [0.02 0.08 0.14];
N = 40;
for c = 1:numel(codes)
    [~, fail] = sweep_depolarizing_logical_error(codes{c}, ps, N);
    for i = 1:numel(ps)
        fprintf(fid, '| %s | %.2f | %d | %.3f |\n', labels{c}, ps(i), N, fail(i));
    end
end

fprintf(fid, '\n## Noisy Syndrome Rounds\n\n');
fprintf(fid, '3-qubit bit-flip code with data error probability 0.08 and syndrome readout error probability 0.08.\n\n');
rounds = [1 3 5];
Nsyn = 400;
fail_rounds = sweep_noisy_syndrome_bitflip(0.08, 0.08, rounds, Nsyn);
fprintf(fid, '| Rounds | Trials | Logical failure |\n');
fprintf(fid, '| ---: | ---: | ---: |\n');
for i = 1:numel(rounds)
    fprintf(fid, '| %d | %d | %.3f |\n', rounds(i), Nsyn, fail_rounds(i));
end

fprintf(fid, '\n## Bacon-Shor Pauli-Frame Snapshot\n\n');
fprintf(fid, '3x3 Bacon-Shor subsystem-code Pauli-frame model under independent depolarizing noise.\n\n');
fprintf(fid, '| p | Trials | Logical failure |\n');
fprintf(fid, '| ---: | ---: | ---: |\n');
Nbacon = 200;
for i = 1:numel(ps)
    failed = 0;
    for trial = 1:Nbacon
        out = simulate_bacon_shor_pauli_once(ps(i));
        failed = failed + (~out.success);
    end
    fprintf(fid, '| %.2f | %d | %.3f |\n', ps(i), Nbacon, failed / Nbacon);
end

fprintf(fid, '\n## Surface-3 Prototype\n\n');
fprintf(fid, 'Small surface-code model with cached minimum-weight correction lookup for X-only, Z-only, and Pauli noise.\n\n');
ps_surface = [0.02 0.08 0.14];
Nsurf = 500;
[~, fail_surface, fail_surface_z, fail_surface_pauli, meta_surface] = surface3_cached_channel_sweep(ps_surface, Nsurf, 7, false);
fprintf(fid, '| Error probability | Trials | X-only failure | Z-only failure | Pauli failure |\n');
fprintf(fid, '| ---: | ---: | ---: | ---: | ---: |\n');
for i = 1:numel(ps_surface)
    fprintf(fid, '| %.2f | %d | %.3f | %.3f | %.3f |\n', ps_surface(i), Nsurf, fail_surface(i), fail_surface_z(i), fail_surface_pauli(i));
end
fprintf(fid, '\nSurface-3 channel cache version: `%s`.\n', meta_surface.cache_version);

fprintf(fid, '\n## Surface-3 Noisy Syndrome Rounds\n\n');
fprintf(fid, 'Classical readout noise on syndrome bits with majority vote across repeated rounds.\n\n');
rounds_surface = [1 3 5];
NsurfSyn = 300;
[~, fail_syn_x, fail_syn_z, meta_surface_syn] = surface3_cached_noisy_syndrome_sweep(0.08, 0.08, rounds_surface, NsurfSyn, 7, false);
fprintf(fid, '| Rounds | Trials | X-channel failure | Z-channel failure |\n');
fprintf(fid, '| ---: | ---: | ---: | ---: |\n');
for i = 1:numel(rounds_surface)
    fprintf(fid, '| %d | %d | %.3f | %.3f |\n', rounds_surface(i), NsurfSyn, fail_syn_x(i), fail_syn_z(i));
end
fprintf(fid, '\nSurface-3 noisy-syndrome cache version: `%s`.\n', meta_surface_syn.cache_version);

fprintf(fid, '\n## Generic Surface-Layout Distance Benchmark\n\n');
fprintf(fid, 'Variable-distance d-by-d surface-layout model under independent Pauli noise. The table compares the default lookup-plus-peeling decoder against MWPM-style matching and bounded syndrome-graph baselines.\n\n');
distances_surface = [3 5 7];
ps_distance = [0 0.03 0.06 0.09];
Ndistance = 8;
distance_decoders = {'min_weight', 'mwpm', 'graph_matching'};
fprintf(fid, '| Decoder | Error probability | Trials | d=3 failure | d=5 failure | d=7 failure |\n');
fprintf(fid, '| --- | ---: | ---: | ---: | ---: | ---: |\n');
for decoder_idx = 1:numel(distance_decoders)
    distance_results = sweep_surface_distance_logical_error(distances_surface, ps_distance, Ndistance, 11, distance_decoders{decoder_idx});
    for i = 1:numel(ps_distance)
        fprintf(fid, '| %s | %.2f | %d | %.3f | %.3f | %.3f |\n', ...
                distance_decoders{decoder_idx}, ps_distance(i), Ndistance, ...
                distance_results.logical_error(1, i), distance_results.logical_error(2, i), ...
                distance_results.logical_error(3, i));
    end
end

fprintf(fid, '\n## Generated Figures\n\n');
fprintf(fid, '- `images/bitflip_logical_error_vs_physical_error.png`\n');
fprintf(fid, '- `images/bitflip_syndrome_distribution.png`\n');
fprintf(fid, '- `images/bitflip_decoder_confusion_matrix.png`\n');
fprintf(fid, '- `images/bitflip_error_weight_distribution.png`\n');
fprintf(fid, '- `images/bitflip_monte_carlo_demo.png`\n');
fprintf(fid, '- `images/qec_recovery_failure_by_error_weight.png`\n');
fprintf(fid, '- `images/qec_depolarizing_logical_error_comparison.png`\n');
fprintf(fid, '- `images/bitflip_noisy_syndrome_rounds.png`\n');
fprintf(fid, '- `images/surface3_logical_error_vs_x_error.png`\n');
fprintf(fid, '- `images/surface3_channel_logical_error_comparison.png`\n');
fprintf(fid, '- `images/surface3_noisy_syndrome_rounds.png`\n');
fprintf(fid, '- `images/surface_distance_logical_error_scaling.png`\n');

clear cleanup;
disp(['Wrote ', report_path]);