simulation_vs_theta_voi_unaware_no_loc.py

import time

import numpy as np
from concurrent.futures import ProcessPoolExecutor
from scipy.stats import chi2
import pandas as pd

import push_success as push
import pull_optimization as pull

import common as cmn
from grid_search_position import grid_search_position
from fim_evaluation import eval_peb_ub


def run_episode(episode_idx: int,
                theta: float, epsilon: float,
                num_frame: int = int(1e2),
                num_push_res: int = cmn.P, num_pull_res: int = cmn.Q,
                num_sym: int = cmn.T, num_sc: int = cmn.F,
                num_ue: int = cmn.U, snr: float = cmn.SNR,
                num_ap: int = cmn.L, num_ant: int = cmn.M,
                ap_pos: np.ndarray = cmn.AP_POS, ap_ant_pos: np.ndarray = cmn.AP_ANT_POS,
                side: float = cmn.ROOM_SIDE,
                voi_dist: str = 'uniform',
                data_size_dist: str = 'geometric',
                data_size_mean: float = cmn.DATA_MEAN,
                grid_test_bool: bool = False, verbose: bool = False,
                debug_mode: bool = False):

    assert epsilon > 0, "Epsilon must be > 0"
    assert 0 <= theta <= 1, "Theta must be within 0 and 1"

    rng = np.random.default_rng(episode_idx)

    ### Episode results variables
    voi_tot = np.zeros(num_frame)
    res_employed = np.zeros(num_frame)
    com_vs_loc = np.zeros(num_frame)
    push_access_rate = np.ones(num_frame)
    pull_access_rate = np.ones(num_frame)

    # Initialization
    ue_pos = 2 * 0.99 * side * (rng.random((num_ue, 2)) - 0.5)

    for t in cmn.std_bar(range(num_frame)):
        # Reset use of resources
        Qu_loc = np.zeros(num_ue, dtype=np.int32)
        Qu_com = np.zeros(num_ue, dtype=np.int32)

        # Update UE position
        ue_pos = cmn.move_ues(ue_pos)

        ## CHANNEL (ONLY PUSH)
        # Compute the large-scale fading for the channel
        beta = cmn.eval_fspl(ue_pos, ap_pos)

        ## Energy detector for estimating beta (ignoring the real Gaussian cross term)
        # Chi-2 noise at the ED (scaled for CE and averaged for the OFDM symbols)
        w_ed = chi2(num_sym * num_sc).rvs((num_ue, num_ap)) / snr / num_sym / num_sc / num_ant
        # Energy detector for a single RE (U x L)
        beta_hat = beta + w_ed

        # Compute the value of information
        value_quantile = rng.random(num_ue)

        # Draw VoI from correct distribution
        if (voi_dist == 'uniform'):
            voi = value_quantile
        # Draw observation size from correct distribution in bits (the geometric gives the bytes)
        if (data_size_dist == 'geometric'):
            data_size = 8 * rng.geometric(1 / data_size_mean, size = num_ue)

        ## Push-based subframe
        # Take active users
        active = np.nonzero(value_quantile > theta)[0]
        # sanity check
        if len(active) == 0:
            continue
        # Simulate push access
        pushed = push.simulate_push(num_push_res, active)
        if verbose:
            print('p', pushed)
        # Compute access rate
        if len(pushed) == 0:
            push_access_rate[t] = 0
            continue
        else:
            push_access_rate[t] = len(pushed) / len(active)

        # All the pushed nodes inform the CPU of the data_size
        # Compute minimum no of resources for communication
        se = cmn.effective_se(beta_hat[pushed])
        Qu_com[pushed] = np.ceil(data_size[pushed] / se / (num_sym - 1) / num_sc).astype(np.int32)

        peb_ub = np.zeros(len(pushed))
        if grid_test_bool:
            xu_worst = np.zeros((len(pushed), 2))
            for i, u in enumerate(pushed):
                peb_ub[i], xu_worst[i] = grid_search_position(beta_hat[u], ap_pos, ap_ant_pos, snr)
        else:
            # When the ED noise is low the worst position is on the edge
            xu_worst = 0.99 * side * np.ones(2)
            for i, u in enumerate(pushed):
                peb_ub[i] = eval_peb_ub(xu_worst, ap_pos, ap_ant_pos, beta_hat[u], snr)

        # Compute minimum no of resources for localization
        Qu_loc[pushed] = 0

        # Print resources for tuning
        if verbose:
            print(f"Qu_com: {Qu_com[pushed]}")
            print(f"Qu_loc: {Qu_loc[pushed]}")
        com_vs_loc[t] = np.sum(Qu_com[pushed] > Qu_loc[pushed]) / len(pushed)

        # Run pull-based subframe scheduler
        scheduled = pushed[pull.optimize_pull(np.ones(len(pushed)), Qu_com[pushed], Qu_loc[pushed], num_pull_res)]

        # Compute access rate
        pull_access_rate[t] = len(scheduled) / len(pushed)

        # Compute overall VoI and scheduled resources
        voi_tot[t] = np.sum(voi[scheduled])
        res_employed[t] = np.sum(np.maximum(Qu_com[scheduled], Qu_loc[scheduled]))


    # Mean results
    avg_com_vs_loc = np.mean(com_vs_loc)
    avg_res_employed = np.mean(res_employed)
    avg_voi_tot = np.mean(voi_tot)
    avg_pull_access_rate = np.mean(pull_access_rate)
    avg_push_access_rate = np.mean(push_access_rate)

    if verbose:
        print(f"Push access rate: {avg_push_access_rate * 100.0:.2f}%")
        print(f"Pull access rate: {avg_pull_access_rate * 100.0:.2f}%")
        print(f"% Comm REs > Loc REs (per UE): {avg_com_vs_loc * 100.0:.2f}%")
        print(f"No. of Pull REs employed: {avg_res_employed} over {num_pull_res}")
        print(f"Avg. tot. VoI (per frame): {avg_voi_tot:.3f}")

    return avg_com_vs_loc, avg_res_employed, avg_voi_tot, avg_pull_access_rate, avg_push_access_rate



if __name__ == '__main__':
    # Parse arguments, if any
    parallel, data_folder, _, overwrite = cmn.common_parser()

    ### MAIN SIMULATION PARAMETERS ###
    theta_vec = cmn.theta_vec
    epsilon = cmn.std_epsilon
    num_frame = cmn.num_frame

    # Check if files exist and load it if there
    prefix = 'voi_unaware_no_loc_vs_theta'
    data_shape = (len(theta_vec), len(cmn.results_label))
    avg_results, filename = cmn.check_data(data_shape, prefix, data_folder, overwrite_flag=overwrite)

    # Start evaluating
    for m, theta in enumerate(theta_vec):
        # Check if data is there
        if overwrite or np.any(np.isnan(avg_results[m])):

            args = (theta, epsilon, num_frame)

            start_time = time.time()
            if parallel:
                with ProcessPoolExecutor() as executor:
                    futures = [executor.submit(run_episode, ep, *args) for ep in range(cmn.E)]
                    results = [f.result() for f in futures]
            else:
                results = []
                for ep in range(cmn.E):
                    print(f'\tEpisode: {ep:02d}/{cmn.E - 1:02d}')
                    results.append(run_episode(ep, *args))

            # Average the results
            avg_results[m] = np.mean(np.array(results), axis=0)

            # Save data (doing it every time is redundant, but it is safe)
            avg_df = pd.DataFrame(avg_results, columns=cmn.results_label)
            avg_df.insert(0, 'Theta', theta_vec)
            avg_df.to_csv(filename, index=False)

            # Print time
            elapsed = time.time() - start_time
            print(f"\t...done in {elapsed:.3f} seconds")

        else:
            print("\t...already done!")
            continue