homotopy-continuation-channel-decoding/python/examples/simulate_error_rate.py

import typing
import numpy as np
import galois
import argparse
from dataclasses import dataclass
from tqdm import tqdm
import pandas as pd

# autopep8: off
import sys
import os
sys.path.append(f"{os.path.dirname(os.path.abspath(__file__))}/../")

from hccd import utility, homotopy_generator, path_tracker
# autopep8: on


@dataclass
class SimulationArgs:
    euler_step_size: float
    euler_max_tries: int
    newton_max_iter: int
    newton_convergence_threshold: float
    sigma: int
    homotopy_iter: int

    max_frames: int
    target_frame_errors: int


def decode(tracker, y, H, args: SimulationArgs) -> np.ndarray:
    x_hat = np.mod(np.round(y), 2).astype('int32')

    s = np.concatenate([y, np.array([0])])
    for i in range(args.homotopy_iter):
        x_hat = np.mod(np.round(s[:-1]), 2).astype('int32')
        if not np.any(np.mod(H @ x_hat, 2)):
            return x_hat

        # if s[-1] > 1.5:
        #     return x_hat

        try:
            s = tracker.step(s)
        except:
            return x_hat

    return x_hat


def simulate_error_rates_for_SNR(H, Eb_N0, args: SimulationArgs) -> typing.Tuple[float, float, float, int]:
    GF = galois.GF(2)
    H_GF = GF(H)
    G = H_GF.null_space()
    k, n = G.shape

    homotopy = homotopy_generator.HomotopyGenerator(H)

    # print(f"G: {homotopy.G}")
    # print(f"F: {homotopy.F}")
    # print(f"H: {homotopy.H}")
    # print(f"DH: {homotopy.DH}")

    tracker = path_tracker.PathTracker(homotopy, args.euler_step_size, args.euler_max_tries,
                                       args.newton_max_iter, args.newton_convergence_threshold, args.sigma)

    num_frames = 0

    bit_errors = 0
    frame_errors = 0
    decoding_failures = 0

    for _ in tqdm(range(args.max_frames)):
        Eb_N0_lin = 10**(Eb_N0 / 10)
        N0 = 1 / (2 * k / n * Eb_N0_lin)
        y = np.zeros(n) + np.sqrt(N0) * np.random.normal(size=n)

        x_hat = decode(tracker, y, H, args)

        bit_errors += np.sum(x_hat != np.zeros(n))
        frame_errors += np.any(x_hat != np.zeros(n))

        if np.any(np.mod(H @ x_hat, 2)):
            decoding_failures += 1

        num_frames += 1

        if frame_errors >= args.target_frame_errors:
            break

    BER = bit_errors / (num_frames * n)
    FER = frame_errors / num_frames
    DFR = decoding_failures / num_frames

    return FER, BER, DFR, frame_errors


def simulate_error_rates(H, Eb_N0_list, args: SimulationArgs) -> typing.Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    FERs = []
    BERs = []
    DFRs = []
    frame_errors_list = []
    for SNR in Eb_N0_list:
        FER, BER, DFR, num_frames = simulate_error_rates_for_SNR(H, SNR, args)
        FERs.append(FER)
        BERs.append(BER)
        DFRs.append(DFR)
        frame_errors_list.append(num_frames)

    return np.array(FERs), np.array(BERs), np.array(DFRs), np.array(frame_errors_list)


def main():
    # Parse command line arguments

    parser = argparse.ArgumentParser()

    parser.add_argument("-c", "--code", type=str, required=True,
                        help="Path to the alist file containing the parity check matrix of the code")
    parser.add_argument("--snr", type=float, required=True, nargs=3,
                        metavar=('start', 'stop', 'step'),
                        help="SNRs to simulate (Eb/N0) in dB")

    parser.add_argument("--max-frames", type=int, default=int(1e6),
                        help="Maximum number of frames to simulate")
    parser.add_argument("--target-frame-errors", type=int, default=200,
                        help="Number of frame errors after which to stop the simulation")

    parser.add_argument("--euler-step-size", type=float,
                        default=0.05, help="Step size for Euler predictor")
    parser.add_argument("--euler-max-tries", type=int, default=5,
                        help="Maximum number of tries for Euler predictor")
    parser.add_argument("--newton-max-iter", type=int, default=5,
                        help="Maximum number of iterations for Newton corrector")
    parser.add_argument("--newton-convergence-threshold", type=float,
                        default=0.01, help="Convergence threshold for Newton corrector")
    parser.add_argument("-s", "--sigma", type=int, default=1,
                        help="Direction in which the path is traced")
    parser.add_argument("-n", "--homotopy-iter", type=int, default=20,
                        help="Number of iterations of the homotopy continuation method to perform for each decoding")

    args = parser.parse_args()

    # TODO: Name this section properly
    # Do stuff

    H = utility.read_alist_file(args.code)

    simulation_args = SimulationArgs(
        euler_step_size=args.euler_step_size,
        euler_max_tries=args.euler_max_tries,
        newton_max_iter=args.newton_max_iter,
        newton_convergence_threshold=args.newton_convergence_threshold,
        sigma=args.sigma,
        homotopy_iter=args.homotopy_iter,
        max_frames=args.max_frames,
        target_frame_errors=args.target_frame_errors)

    SNRs = np.arange(args.snr[0], args.snr[1], args.snr[2])
    FERs, BERs, DFRs, frame_errors_list = simulate_error_rates(
        H, SNRs, simulation_args)

    df = pd.DataFrame({"SNR": SNRs, "FER": FERs, "BER": BERs, "DFR": DFRs, "frame_errors": frame_errors_list})
    print(df)

if __name__ == "__main__":
    main()