diff --git a/data/de_results_z128.json b/data/de_results_z128.json new file mode 100644 index 0000000..43b8ea7 --- /dev/null +++ b/data/de_results_z128.json @@ -0,0 +1,284 @@ +{ + "z": 128, + "cn_mode": "normalized", + "alpha": 0.875, + "degrees": [ + 7, + 4, + 4, + 4, + 4, + 3, + 3, + 3 + ], + "de_threshold": 3.2093749999999996, + "matrix": [ + [ + 7, + 92, + -1, + 25, + 73, + 89, + -1, + -1 + ], + [ + 125, + 33, + 114, + -1, + -1, + -1, + 104, + -1 + ], + [ + 6, + -1, + 67, + 57, + 116, + -1, + -1, + -1 + ], + [ + 113, + -1, + -1, + 108, + 74, + -1, + -1, + -1 + ], + [ + 28, + 119, + 35, + -1, + 88, + 20, + -1, + 124 + ], + [ + 120, + 35, + -1, + 9, + -1, + 100, + 72, + 23 + ], + [ + 63, + -1, + 98, + -1, + -1, + -1, + 48, + 98 + ] + ], + "matrix_checks": { + "col_degrees": [ + 7, + 4, + 4, + 4, + 4, + 3, + 3, + 3 + ], + "full_rank": true, + "actual_rank": 896, + "expected_rank": 896, + "parity_rank": true, + "girth": 8, + "encodable": true + }, + "girth": 8, + "fer_results": { + "lam_s_points": [ + 1.5, + 2.0, + 3.0, + 4.0, + 5.0, + 7.0, + 10.0 + ], + "optimized": { + "1.5": { + "fer": 1.0, + "ber": 0.685, + "avg_iter": 30.0 + }, + "2.0": { + "fer": 1.0, + "ber": 0.679375, + "avg_iter": 30.0 + }, + "3.0": { + "fer": 0.74, + "ber": 0.39109375, + "avg_iter": 22.46 + }, + "4.0": { + "fer": 0.24, + "ber": 0.14015625, + "avg_iter": 7.96 + }, + "5.0": { + "fer": 0.04, + "ber": 0.0215625, + "avg_iter": 2.16 + }, + "7.0": { + "fer": 0.0, + "ber": 0.0, + "avg_iter": 1.0 + }, + "10.0": { + "fer": 0.0, + "ber": 0.0, + "avg_iter": 1.0 + } + }, + "peg_ring": { + "1.5": { + "fer": 1.0, + "ber": 0.7315625, + "avg_iter": 30.0 + }, + "2.0": { + "fer": 1.0, + "ber": 0.631875, + "avg_iter": 30.0 + }, + "3.0": { + "fer": 0.94, + "ber": 0.5009375, + "avg_iter": 28.26 + }, + "4.0": { + "fer": 0.32, + "ber": 0.16390625, + "avg_iter": 10.28 + }, + "5.0": { + "fer": 0.04, + "ber": 0.01984375, + "avg_iter": 2.16 + }, + "7.0": { + "fer": 0.0, + "ber": 0.0, + "avg_iter": 1.0 + }, + "10.0": { + "fer": 0.0, + "ber": 0.0, + "avg_iter": 1.0 + } + } + }, + "peg_matrix": [ + [ + 37, + 123, + 82, + -1, + -1, + -1, + -1, + -1 + ], + [ + 4, + 56, + 18, + 124, + -1, + -1, + -1, + -1 + ], + [ + 95, + -1, + 86, + 107, + -1, + -1, + -1, + 49 + ], + [ + 18, + -1, + -1, + 107, + 66, + -1, + -1, + -1 + ], + [ + 10, + 81, + -1, + -1, + 94, + 52, + -1, + -1 + ], + [ + 83, + -1, + -1, + -1, + -1, + 68, + 24, + -1 + ], + [ + 102, + -1, + -1, + -1, + -1, + -1, + 77, + 118 + ] + ], + "peg_girth": 8, + "peg_matrix_checks": { + "col_degrees": [ + 7, + 3, + 3, + 3, + 2, + 2, + 2, + 2 + ], + "full_rank": true, + "actual_rank": 896, + "expected_rank": 896, + "parity_rank": true, + "girth": 8, + "encodable": true + }, + "n_frames": 50 +} \ No newline at end of file diff --git a/data/plots/base_matrix_heatmap.png b/data/plots/base_matrix_heatmap.png new file mode 100644 index 0000000..fde6358 Binary files /dev/null and b/data/plots/base_matrix_heatmap.png differ diff --git a/data/plots/degree_distributions.png b/data/plots/degree_distributions.png new file mode 100644 index 0000000..407575e Binary files /dev/null and b/data/plots/degree_distributions.png differ diff --git a/data/plots/fer_comparison.png b/data/plots/fer_comparison.png new file mode 100644 index 0000000..526816f Binary files /dev/null and b/data/plots/fer_comparison.png differ diff --git a/data/plots/threshold_comparison.png b/data/plots/threshold_comparison.png new file mode 100644 index 0000000..e52cf84 Binary files /dev/null and b/data/plots/threshold_comparison.png differ diff --git a/data/plots/z128_fer_comparison.png b/data/plots/z128_fer_comparison.png new file mode 100644 index 0000000..13f8e60 Binary files /dev/null and b/data/plots/z128_fer_comparison.png differ diff --git a/model/plot_de_results.py b/model/plot_de_results.py new file mode 100644 index 0000000..1061d6e --- /dev/null +++ b/model/plot_de_results.py @@ -0,0 +1,318 @@ +#!/usr/bin/env python3 +"""Generate presentation plots from density evolution results.""" + +import json +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from pathlib import Path + +RESULTS_PATH = Path(__file__).parent.parent / 'data' / 'de_results.json' +RESULTS_Z128_PATH = Path(__file__).parent.parent / 'data' / 'de_results_z128.json' +OUT_DIR = Path(__file__).parent.parent / 'data' / 'plots' + +# Presentation-friendly style +plt.rcParams.update({ + 'font.size': 13, + 'axes.titlesize': 15, + 'axes.labelsize': 14, + 'legend.fontsize': 11, + 'xtick.labelsize': 12, + 'ytick.labelsize': 12, + 'figure.dpi': 150, + 'savefig.bbox': 'tight', + 'savefig.pad_inches': 0.15, +}) + +COLORS = { + 'original': '#d62728', # red + 'peg_ring': '#1f77b4', # blue + 'optimized': '#2ca02c', # green +} + +LABELS = { + 'original': 'Original staircase [7,2,2,2,2,2,2,1]', + 'peg_ring': 'PEG ring [7,3,3,3,2,2,2,2]', + 'optimized': 'DE-optimized [7,4,4,4,4,3,3,3]', +} + +MARKERS = { + 'original': 's', + 'peg_ring': '^', + 'optimized': 'o', +} + + +def load_results(): + with open(RESULTS_PATH) as f: + return json.load(f) + + +def plot_fer_comparison(data): + """FER vs lambda_s for all three matrices.""" + fig, ax = plt.subplots(figsize=(8, 5.5)) + + fer_data = data['fer_comparison'] + lam_s = fer_data['lam_s_points'] + + for key in ['original', 'peg_ring', 'optimized']: + fer_vals = [fer_data[key][str(l)]['fer'] for l in lam_s] + # Replace 0 with small value for log scale + fer_plot = [max(f, 2e-3) for f in fer_vals] + ax.semilogy(lam_s, fer_plot, + color=COLORS[key], marker=MARKERS[key], markersize=8, + linewidth=2, label=LABELS[key], + markeredgecolor='white', markeredgewidth=0.8) + + ax.set_xlabel(r'Signal photons/slot ($\lambda_s$)') + ax.set_ylabel('Frame Error Rate (FER)') + ax.set_title('FER Comparison: Rate-1/8 QC-LDPC (n=256, k=32)') + ax.legend(loc='upper right', framealpha=0.9) + ax.set_xlim(1.5, 10.5) + ax.set_ylim(1e-3, 1.0) + ax.grid(True, alpha=0.3, which='both') + ax.set_xticks([2, 3, 4, 5, 7, 10]) + ax.set_xticklabels(['2', '3', '4', '5', '7', '10']) + + # Annotate the improvement + ax.annotate('13x fewer\nframe errors', + xy=(5, 0.01), xytext=(6.5, 0.06), + fontsize=11, color=COLORS['optimized'], + arrowprops=dict(arrowstyle='->', color=COLORS['optimized'], lw=1.5), + ha='center') + + fig.savefig(OUT_DIR / 'fer_comparison.png') + plt.close(fig) + print(f' Saved fer_comparison.png') + + +def plot_threshold_bars(data): + """Bar chart of DE thresholds.""" + fig, ax = plt.subplots(figsize=(7, 5)) + + thresholds = { + 'Original\nstaircase': data['reference_thresholds']['Original staircase [7,2,2,2,2,2,2,1]'], + 'PEG\nring': data['reference_thresholds']['PEG ring [7,3,3,3,2,2,2,2]'], + 'DE-\noptimized': data['best_threshold'], + } + # Shannon limit for rate 1/8 Poisson channel ~ 0.47 photons/slot + shannon_limit = 0.47 + + names = list(thresholds.keys()) + vals = list(thresholds.values()) + colors = [COLORS['original'], COLORS['peg_ring'], COLORS['optimized']] + + bars = ax.bar(names, vals, color=colors, width=0.55, edgecolor='white', linewidth=1.5) + + # Value labels on bars + for bar, val in zip(bars, vals): + ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, + f'{val:.1f}', ha='center', va='bottom', fontweight='bold', fontsize=14) + + # Shannon limit line + ax.axhline(y=shannon_limit, color='black', linestyle='--', linewidth=1.5, alpha=0.6) + ax.text(2.42, shannon_limit + 0.12, f'Shannon limit ({shannon_limit})', + fontsize=11, ha='right', va='bottom', style='italic', alpha=0.7) + + ax.set_ylabel(r'DE Threshold ($\lambda_s^*$, photons/slot)') + ax.set_title('Decoding Threshold Comparison') + ax.set_ylim(0, 6.5) + ax.grid(True, alpha=0.2, axis='y') + + # Improvement annotation + improvement_db = 10 * np.log10(thresholds['Original\nstaircase'] / thresholds['DE-\noptimized']) + ax.annotate(f'{improvement_db:.1f} dB\nimprovement', + xy=(2, thresholds['DE-\noptimized']), + xytext=(2.35, 4.2), + fontsize=12, fontweight='bold', color=COLORS['optimized'], + arrowprops=dict(arrowstyle='->', color=COLORS['optimized'], lw=1.5), + ha='center') + + fig.savefig(OUT_DIR / 'threshold_comparison.png') + plt.close(fig) + print(f' Saved threshold_comparison.png') + + +def plot_degree_distribution(data): + """VN degree distribution comparison.""" + fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True) + + distributions = { + 'Original staircase': [7, 2, 2, 2, 2, 2, 2, 1], + 'PEG ring': [7, 3, 3, 3, 2, 2, 2, 2], + 'DE-optimized': data['best_degrees'], + } + color_keys = ['original', 'peg_ring', 'optimized'] + col_labels = [f'c{i}' for i in range(8)] + + for ax, (name, degrees), ckey in zip(axes, distributions.items(), color_keys): + colors_bar = [COLORS[ckey] if i > 0 else '#888888' for i in range(8)] + ax.bar(col_labels, degrees, color=colors_bar, edgecolor='white', linewidth=1) + ax.set_title(name, fontsize=12, fontweight='bold') + ax.set_xlabel('Base matrix column') + ax.set_ylim(0, 8) + ax.grid(True, alpha=0.2, axis='y') + + # Annotate average parity degree + avg_parity = np.mean(degrees[1:]) + ax.text(4, 7.2, f'avg parity dv = {avg_parity:.1f}', + ha='center', fontsize=10, style='italic') + + axes[0].set_ylabel('Variable node degree') + fig.suptitle('VN Degree Distributions (col 0 = info, cols 1-7 = parity)', + fontsize=13, y=1.02) + fig.tight_layout() + + fig.savefig(OUT_DIR / 'degree_distributions.png') + plt.close(fig) + print(f' Saved degree_distributions.png') + + +def plot_base_matrix_heatmap(data): + """Heatmap of the constructed base matrix.""" + fig, ax = plt.subplots(figsize=(7, 5)) + + H = np.array(data['constructed_matrix']) + m_base, n_base = H.shape + + # Create display matrix: shift values where connected, NaN where not + display = np.full_like(H, dtype=float, fill_value=np.nan) + for r in range(m_base): + for c in range(n_base): + if H[r, c] >= 0: + display[r, c] = H[r, c] + + im = ax.imshow(display, cmap='YlGnBu', aspect='auto', vmin=0, vmax=31) + + # Add text annotations + for r in range(m_base): + for c in range(n_base): + if H[r, c] >= 0: + ax.text(c, r, str(H[r, c]), ha='center', va='center', + fontsize=11, fontweight='bold', color='black') + else: + ax.text(c, r, '-', ha='center', va='center', + fontsize=11, color='#cccccc') + + ax.set_xticks(range(n_base)) + ax.set_xticklabels([f'c{i}' for i in range(n_base)]) + ax.set_yticks(range(m_base)) + ax.set_yticklabels([f'r{i}' for i in range(m_base)]) + ax.set_xlabel('Column (c0=info, c1-c7=parity)') + ax.set_ylabel('Row (check node group)') + ax.set_title(f'Optimized Base Matrix (Z=32, girth={data["matrix_checks"]["girth"]})') + + cbar = fig.colorbar(im, ax=ax, shrink=0.8, label='Circulant shift') + + fig.savefig(OUT_DIR / 'base_matrix_heatmap.png') + plt.close(fig) + print(f' Saved base_matrix_heatmap.png') + + +def load_z128_results(): + """Load Z=128 results if available.""" + if RESULTS_Z128_PATH.exists(): + with open(RESULTS_Z128_PATH) as f: + return json.load(f) + return None + + +def plot_z128_comparison(data_z32, data_z128): + """ + FER vs lambda_s comparison between Z=32 and Z=128 for the + DE-optimized degree distribution [7,4,4,4,4,3,3,3]. + + Shows both lifting factors on the same axes, with the optimized + matrix in each case plus the PEG ring reference. + """ + fig, ax = plt.subplots(figsize=(9, 6)) + + # Colors and styles for Z=32 vs Z=128 + z32_color = '#2ca02c' # green (matches optimized) + z128_color = '#9467bd' # purple + z32_peg_color = '#1f77b4' # blue (matches peg_ring) + z128_peg_color = '#ff7f0e' # orange + + # Z=32 optimized data (from de_results.json) + fer_z32 = data_z32['fer_comparison'] + lam_s_z32 = fer_z32['lam_s_points'] + fer_z32_opt = [fer_z32['optimized'][str(l)]['fer'] for l in lam_s_z32] + fer_z32_peg = [fer_z32['peg_ring'][str(l)]['fer'] for l in lam_s_z32] + + # Z=128 data + fer_z128 = data_z128['fer_results'] + lam_s_z128 = fer_z128['lam_s_points'] + fer_z128_opt = [fer_z128['optimized'][str(l)]['fer'] for l in lam_s_z128] + fer_z128_peg = [fer_z128['peg_ring'][str(l)]['fer'] for l in lam_s_z128] + + # Replace 0 with small value for log scale + floor = 5e-3 + + # Plot Z=32 curves (solid lines) + ax.semilogy(lam_s_z32, [max(f, floor) for f in fer_z32_opt], + color=z32_color, marker='o', markersize=8, linewidth=2, + label='Optimized Z=32 (n=256, offset MS)', + markeredgecolor='white', markeredgewidth=0.8) + ax.semilogy(lam_s_z32, [max(f, floor) for f in fer_z32_peg], + color=z32_peg_color, marker='^', markersize=8, linewidth=2, + label='PEG ring Z=32 (n=256, offset MS)', + markeredgecolor='white', markeredgewidth=0.8, + linestyle='--', alpha=0.7) + + # Plot Z=128 curves (dashed lines, larger markers) + ax.semilogy(lam_s_z128, [max(f, floor) for f in fer_z128_opt], + color=z128_color, marker='D', markersize=9, linewidth=2.5, + label=f'Optimized Z=128 (n=1024, norm MS, a={data_z128["alpha"]})', + markeredgecolor='white', markeredgewidth=0.8) + ax.semilogy(lam_s_z128, [max(f, floor) for f in fer_z128_peg], + color=z128_peg_color, marker='v', markersize=9, linewidth=2.5, + label=f'PEG ring Z=128 (n=1024, norm MS, a={data_z128["alpha"]})', + markeredgecolor='white', markeredgewidth=0.8, + linestyle='--', alpha=0.7) + + ax.set_xlabel(r'Signal photons/slot ($\lambda_s$)') + ax.set_ylabel('Frame Error Rate (FER)') + ax.set_title('FER Comparison: Z=32 vs Z=128\nDE-optimized [7,4,4,4,4,3,3,3]') + ax.legend(loc='upper right', framealpha=0.9, fontsize=10) + ax.set_xlim(1.0, 11.0) + ax.set_ylim(floor / 2, 1.1) + ax.grid(True, alpha=0.3, which='both') + ax.set_xticks([1.5, 2, 3, 4, 5, 7, 10]) + ax.set_xticklabels(['1.5', '2', '3', '4', '5', '7', '10']) + + # Add Z=128 matrix info annotation + girth_128 = data_z128.get('girth', '?') + ax.text(0.02, 0.02, + f'Z=128: girth={girth_128}, n=1024, k=128\n' + f'Z=32: girth={data_z32["matrix_checks"]["girth"]}, n=256, k=32', + transform=ax.transAxes, fontsize=9, verticalalignment='bottom', + bbox=dict(boxstyle='round,pad=0.3', facecolor='wheat', alpha=0.5)) + + fig.savefig(OUT_DIR / 'z128_fer_comparison.png') + plt.close(fig) + print(f' Saved z128_fer_comparison.png') + + +def main(): + OUT_DIR.mkdir(parents=True, exist_ok=True) + data = load_results() + + print('Generating plots...') + plot_fer_comparison(data) + plot_threshold_bars(data) + plot_degree_distribution(data) + plot_base_matrix_heatmap(data) + + # Z=128 comparison plot + data_z128 = load_z128_results() + if data_z128 is not None: + plot_z128_comparison(data, data_z128) + else: + print(' Skipping Z=128 comparison (data/de_results_z128.json not found)') + + print(f'\nAll plots saved to {OUT_DIR}/') + + +if __name__ == '__main__': + main() diff --git a/model/run_z128.py b/model/run_z128.py new file mode 100644 index 0000000..bd12fa6 --- /dev/null +++ b/model/run_z128.py @@ -0,0 +1,128 @@ +#!/usr/bin/env python3 +""" +Standalone script to construct a Z=128 QC-LDPC matrix with the DE-optimized +degree distribution [7,4,4,4,4,3,3,3] and run FER validation. + +Saves results to data/de_results_z128.json. +""" + +import numpy as np +import json +import os +import sys +import time + +sys.path.insert(0, os.path.dirname(__file__)) + +from density_evolution import ( + construct_base_matrix, verify_matrix, validate_matrix, + build_de_profile, compute_threshold, make_profile, +) +from ldpc_analysis import build_peg_matrix + +np.random.seed(42) + +Z = 128 +DEGREES = [7, 4, 4, 4, 4, 3, 3, 3] +CN_MODE = 'normalized' +ALPHA = 0.875 +LAM_B = 0.1 + +print("=" * 70) +print(f"Z={Z} LDPC Matrix Construction and Validation") +print(f"Degrees: {DEGREES}") +print(f"CN mode: {CN_MODE}, alpha: {ALPHA}") +print("=" * 70) + +# Step 1: Construct Z=128 matrix +print("\n--- Step 1: Matrix construction (Z=128, reduced trials) ---") +t0 = time.time() +H_base, girth = construct_base_matrix(DEGREES, z=Z, n_trials=500) +t1 = time.time() +print(f" Construction time: {t1 - t0:.1f}s") + +# Step 2: Verify matrix +print("\n--- Step 2: Matrix verification ---") +t0 = time.time() +checks = verify_matrix(H_base, z=Z) +t1 = time.time() +print(f" Verification time: {t1 - t0:.1f}s") +print(f" Z={Z} matrix: girth={girth}, rank={checks['actual_rank']}/{checks['expected_rank']}") +print(f" Full rank: {checks['full_rank']}, Encodable: {checks['encodable']}") +print(f" Column degrees: {checks['col_degrees']}") +print(f" Base matrix:\n{H_base}") + +# Step 3: DE threshold (Z-independent, quick check) +print("\n--- Step 3: DE threshold check ---") +profile = build_de_profile(DEGREES, m_base=7) +threshold = compute_threshold(profile, lam_b=LAM_B, z_pop=5000, tol=0.5, + cn_mode=CN_MODE, alpha=ALPHA) +print(f" DE threshold (normalized, alpha={ALPHA}): {threshold:.2f} photons/slot") + +# Step 4: FER validation with reduced frames +print("\n--- Step 4: FER validation (n_frames=50) ---") +lam_s_points = [1.5, 2.0, 3.0, 4.0, 5.0, 7.0, 10.0] + +# Optimized matrix at Z=128 +print("\nRunning optimized Z=128 FER...") +t0 = time.time() +fer_opt = validate_matrix(H_base, lam_s_points, n_frames=50, lam_b=LAM_B, + z=Z, cn_mode=CN_MODE, alpha=ALPHA) +t1 = time.time() +print(f" Optimized FER time: {t1 - t0:.1f}s") + +# Also build reference matrices at Z=128 for comparison +print("\nRunning PEG ring Z=128 FER...") +# Construct a PEG-ring style matrix at Z=128 +from ldpc_analysis import build_peg_matrix as _build_peg_z32 +# The PEG ring base matrix structure is Z-independent; just reconstruct with Z=128 shifts +peg_degrees = [7, 3, 3, 3, 2, 2, 2, 2] +H_peg_128, peg_girth = construct_base_matrix(peg_degrees, z=Z, n_trials=500) +peg_checks = verify_matrix(H_peg_128, z=Z) +print(f" PEG Z=128: girth={peg_girth}, rank={peg_checks['actual_rank']}/{peg_checks['expected_rank']}") + +t0 = time.time() +fer_peg = validate_matrix(H_peg_128, lam_s_points, n_frames=50, lam_b=LAM_B, + z=Z, cn_mode=CN_MODE, alpha=ALPHA) +t1 = time.time() +print(f" PEG FER time: {t1 - t0:.1f}s") + +# Print results table +print(f"\n{'lam_s':>8s} {'Opt Z=128':>12s} {'PEG Z=128':>12s}") +print("-" * 36) +for lam_s in lam_s_points: + f_opt = fer_opt[lam_s]['fer'] + f_peg = fer_peg[lam_s]['fer'] + print(f"{lam_s:8.1f} {f_opt:12.3f} {f_peg:12.3f}") + +# Save results +output = { + 'z': Z, + 'cn_mode': CN_MODE, + 'alpha': ALPHA, + 'degrees': DEGREES, + 'de_threshold': float(threshold), + 'matrix': H_base.tolist(), + 'matrix_checks': { + k: (v.item() if hasattr(v, 'item') else v) for k, v in checks.items() + }, + 'girth': int(girth) if girth < float('inf') else None, + 'fer_results': { + 'lam_s_points': lam_s_points, + 'optimized': {str(k): v for k, v in fer_opt.items()}, + 'peg_ring': {str(k): v for k, v in fer_peg.items()}, + }, + 'peg_matrix': H_peg_128.tolist(), + 'peg_girth': int(peg_girth) if peg_girth < float('inf') else None, + 'peg_matrix_checks': { + k: (v.item() if hasattr(v, 'item') else v) for k, v in peg_checks.items() + }, + 'n_frames': 50, +} + +out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'data') +os.makedirs(out_dir, exist_ok=True) +out_path = os.path.join(out_dir, 'de_results_z128.json') +with open(out_path, 'w') as f: + json.dump(output, f, indent=2, default=str) +print(f"\nResults saved to {out_path}")