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.where(y >= 0, 0, 1)

    s = np.concatenate([y, np.array([0])])
    for i in range(args.homotopy_iter):
        x_hat = np.where(s[:-1] >= 0, 0, 1)
        if not np.any(np.mod(H @ x_hat, 2)):
            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

    num_frames = 0

    bit_errors = 0
    frame_errors = 0
    decoding_failures = 0

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

    for _ in tqdm(range(args.max_frames)):
        Eb_N0_lin = 10**(Eb_N0 / 10)
        N0 = 1 / (2 * k / n * Eb_N0_lin)

        u = np.random.randint(2, size=k)
        # u = np.zeros(shape=k).astype(np.int32)
        c = np.array(GF(u) @ G)
        x = 1 - 2*c

        y = x + np.sqrt(N0) * np.random.normal(size=n)

        homotopy.update_received(y)

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

        if np.any(np.mod(H @ c_hat, 2)):
            tracker.set_sigma(-1 * args.sigma)
            c_hat = decode(tracker, y, H, args)
            tracker.set_sigma(args.sigma)

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

        bit_errors += np.sum(c_hat != c)
        frame_errors += np.any(c_hat != c)
        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) -> pd.DataFrame:
    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)
        print(pd.DataFrame({"SNR": Eb_N0_list[:len(FERs)], "FER": FERs, "BER": BERs,
                            "DFR": DFRs, "frame_errors": frame_errors_list}))

    return pd.DataFrame({"SNR": Eb_N0_list, "FER": FERs, "BER": BERs,
                         "DFR": DFRs, "frame_errors": 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])
    df = simulate_error_rates(H, SNRs, simulation_args)

    print(df)


if __name__ == "__main__":
    main()