diff --git a/model/density_evolution.py b/model/density_evolution.py index c1ea23f..480f196 100644 --- a/model/density_evolution.py +++ b/model/density_evolution.py @@ -702,9 +702,240 @@ def verify_matrix(H_base, z=32): # ============================================================================= -# CLI placeholder (will be extended in later tasks) +# FER Validation # ============================================================================= +def validate_matrix(H_base, lam_s_points, n_frames=200, lam_b=0.1, max_iter=30): + """ + Run full Monte Carlo FER simulation for a base matrix. + + Returns dict mapping lam_s -> {fer, ber, avg_iter}. + """ + from ldpc_analysis import generic_decode, peg_encode + + z = Z + k = z + H_full = _build_full_h_from_base(H_base, z) + + results = {} + for lam_s in lam_s_points: + bit_errors = 0 + frame_errors = 0 + total_bits = 0 + total_iter = 0 + + for _ in range(n_frames): + info = np.random.randint(0, 2, k).astype(np.int8) + try: + codeword = peg_encode(info, H_base, H_full, z=z) + except Exception: + frame_errors += 1 + total_bits += k + continue + + from ldpc_sim import poisson_channel, quantize_llr + llr_float, _ = poisson_channel(codeword, lam_s, lam_b) + llr_q = quantize_llr(llr_float) + + decoded, converged, iters, sw = generic_decode( + llr_q, H_base, z=z, max_iter=max_iter, q_bits=Q_BITS + ) + total_iter += iters + + errs = np.sum(decoded != info) + bit_errors += errs + total_bits += k + if errs > 0: + frame_errors += 1 + + ber = bit_errors / total_bits if total_bits > 0 else 0 + fer = frame_errors / n_frames + avg_iter = total_iter / n_frames + + results[lam_s] = { + 'fer': float(fer), + 'ber': float(ber), + 'avg_iter': float(avg_iter), + } + + return results + + +def run_full_pipeline(lam_b=0.1, seed=42): + """ + Full pipeline: optimize -> construct -> validate. + """ + np.random.seed(seed) + + print("=" * 70) + print("DENSITY EVOLUTION FULL PIPELINE") + print("=" * 70) + + # Step 1: Compute threshold for reference matrices + print("\n--- Step 1: Reference thresholds ---") + from ldpc_analysis import build_peg_matrix + ref_matrices = { + 'Original staircase [7,2,2,2,2,2,2,1]': H_BASE, + 'PEG ring [7,3,3,3,2,2,2,2]': build_peg_matrix(z=32)[0], + } + ref_thresholds = {} + for name, H_b in ref_matrices.items(): + profile = make_profile(H_b) + threshold = compute_threshold(profile, lam_b=lam_b, z_pop=20000, tol=0.25) + ref_thresholds[name] = threshold + print(f" {name}: threshold = {threshold:.2f} photons/slot") + + # Step 2: Run optimizer + print("\n--- Step 2: Degree distribution optimization ---") + opt_results = optimize_degree_distribution( + m_base=7, lam_b=lam_b, top_k=5, + z_pop_coarse=10000, z_pop_fine=20000, tol=0.25 + ) + + if not opt_results: + print("ERROR: No valid distributions found") + return None + + best_degrees, best_threshold = opt_results[0] + print(f"\nBest distribution: {best_degrees}") + print(f"Best threshold: {best_threshold:.2f} photons/slot") + + # Step 3: Construct matrix + print("\n--- Step 3: Matrix construction ---") + H_opt, girth = construct_base_matrix(best_degrees, z=32, n_trials=2000) + checks = verify_matrix(H_opt, z=32) + print(f" Constructed matrix: girth={girth}, rank={checks['actual_rank']}/{checks['expected_rank']}") + print(f" Encodable: {checks['encodable']}, Full rank: {checks['full_rank']}") + print(f" Base matrix:\n{H_opt}") + + # Step 4: FER validation + print("\n--- Step 4: FER validation ---") + lam_s_points = [2.0, 3.0, 4.0, 5.0, 7.0, 10.0] + + print(f"\n{'lam_s':>8s} {'Optimized':>12s} {'Original':>12s} {'PEG ring':>12s}") + print("-" * 50) + + fer_opt = validate_matrix(H_opt, lam_s_points, n_frames=200, lam_b=lam_b) + fer_orig = validate_matrix(H_BASE, lam_s_points, n_frames=200, lam_b=lam_b) + fer_peg = validate_matrix(ref_matrices['PEG ring [7,3,3,3,2,2,2,2]'], + lam_s_points, n_frames=200, lam_b=lam_b) + + for lam_s in lam_s_points: + f_opt = fer_opt[lam_s]['fer'] + f_orig = fer_orig[lam_s]['fer'] + f_peg = fer_peg[lam_s]['fer'] + print(f"{lam_s:8.1f} {f_opt:12.3f} {f_orig:12.3f} {f_peg:12.3f}") + + # Save results + output = { + 'reference_thresholds': ref_thresholds, + 'optimizer_results': [(d, t) for d, t in opt_results], + 'best_degrees': best_degrees, + 'best_threshold': float(best_threshold), + 'constructed_matrix': H_opt.tolist(), + 'matrix_checks': {k: (v if not isinstance(v, (np.integer, np.floating, np.bool_)) else v.item()) + for k, v in checks.items()}, + 'fer_comparison': { + 'lam_s_points': lam_s_points, + 'optimized': {str(k): v for k, v in fer_opt.items()}, + 'original': {str(k): v for k, v in fer_orig.items()}, + 'peg_ring': {str(k): v for k, v in fer_peg.items()}, + }, + } + + 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.json') + with open(out_path, 'w') as f: + json.dump(output, f, indent=2, default=str) + print(f"\nResults saved to {out_path}") + + return output + + +# ============================================================================= +# CLI +# ============================================================================= + +def main(): + parser = argparse.ArgumentParser( + description='Density Evolution Optimizer for QC-LDPC Codes', + ) + subparsers = parser.add_subparsers(dest='command') + + # threshold + p_thresh = subparsers.add_parser('threshold', help='Compute threshold for known matrices') + p_thresh.add_argument('--z-pop', type=int, default=20000) + p_thresh.add_argument('--tol', type=float, default=0.25) + p_thresh.add_argument('--seed', type=int, default=42) + + # optimize + p_opt = subparsers.add_parser('optimize', help='Run optimizer, print top-10') + p_opt.add_argument('--top-k', type=int, default=10) + p_opt.add_argument('--seed', type=int, default=42) + + # construct + p_con = subparsers.add_parser('construct', help='Build matrix for given degrees') + p_con.add_argument('degrees', type=str, help='Comma-separated degrees, e.g. 7,3,3,3,2,2,2,2') + p_con.add_argument('--n-trials', type=int, default=2000) + p_con.add_argument('--seed', type=int, default=42) + + # validate + p_val = subparsers.add_parser('validate', help='FER comparison') + p_val.add_argument('--n-frames', type=int, default=200) + p_val.add_argument('--seed', type=int, default=42) + + # full + p_full = subparsers.add_parser('full', help='Full pipeline') + p_full.add_argument('--seed', type=int, default=42) + + args = parser.parse_args() + + if args.command == 'threshold': + np.random.seed(args.seed) + from ldpc_analysis import build_peg_matrix + matrices = { + 'Original staircase [7,2,2,2,2,2,2,1]': H_BASE, + 'PEG ring [7,3,3,3,2,2,2,2]': build_peg_matrix(z=32)[0], + } + print("DE Threshold Computation") + print("-" * 50) + for name, H_b in matrices.items(): + profile = make_profile(H_b) + threshold = compute_threshold(profile, lam_b=0.1, z_pop=args.z_pop, tol=args.tol) + print(f" {name}: {threshold:.2f} photons/slot") + + elif args.command == 'optimize': + np.random.seed(args.seed) + optimize_degree_distribution(m_base=7, lam_b=0.1, top_k=args.top_k) + + elif args.command == 'construct': + np.random.seed(args.seed) + degrees = [int(x) for x in args.degrees.split(',')] + H_base, girth = construct_base_matrix(degrees, z=32, n_trials=args.n_trials) + checks = verify_matrix(H_base, z=32) + print(f"Constructed matrix for degrees {degrees}:") + print(f" Girth: {girth}") + print(f" Rank: {checks['actual_rank']}/{checks['expected_rank']}") + print(f" Encodable: {checks['encodable']}") + print(f"Base matrix:") + print(H_base) + + elif args.command == 'validate': + np.random.seed(args.seed) + lam_s_points = [2.0, 3.0, 4.0, 5.0, 7.0, 10.0] + print("FER Validation: Original staircase") + results = validate_matrix(H_BASE, lam_s_points, n_frames=args.n_frames) + for lam_s in lam_s_points: + r = results[lam_s] + print(f" lam_s={lam_s:.1f}: FER={r['fer']:.3f}, BER={r['ber']:.6f}") + + elif args.command == 'full': + run_full_pipeline(seed=args.seed) + + else: + parser.print_help() + + if __name__ == '__main__': - print("Density Evolution Optimizer") - print("Run with subcommands: threshold, optimize, construct, validate, full") + main() diff --git a/model/test_density_evolution.py b/model/test_density_evolution.py index 86f1ef3..80375e7 100644 --- a/model/test_density_evolution.py +++ b/model/test_density_evolution.py @@ -175,3 +175,27 @@ class TestPEGBaseMatrixConstructor: target = [7, 3, 3, 3, 2, 2, 2, 2] H_base, girth = construct_base_matrix(target, z=32, n_trials=500) assert girth >= 4, f"Girth {girth} should be >= 4" + + +class TestFERValidationAndCLI: + """Tests for FER validation and CLI.""" + + def test_validate_returns_results(self): + """validate_matrix should return FER results dict.""" + from density_evolution import validate_matrix + from ldpc_sim import H_BASE + np.random.seed(42) + results = validate_matrix(H_BASE, lam_s_points=[10.0], n_frames=10, lam_b=0.1) + assert 10.0 in results, f"Expected key 10.0, got {list(results.keys())}" + assert 'fer' in results[10.0] + assert 0.0 <= results[10.0]['fer'] <= 1.0 + + def test_cli_threshold(self): + """CLI threshold subcommand should exit 0.""" + import subprocess + result = subprocess.run( + ['python3', 'model/density_evolution.py', 'threshold', '--z-pop', '5000', '--tol', '1.0'], + capture_output=True, text=True, timeout=120, + ) + assert result.returncode == 0, f"CLI failed: {result.stderr}" + assert 'threshold' in result.stdout.lower() or 'photon' in result.stdout.lower()