data_fig07_08_cellfree.py

########################################
#   data_fig07_08_bcf.py
#
#   Description. Script used to generate data for Figs. 7 and 8 of the paper
#   concerning cellfree curves (CF-SUCRe). You should choose the estimator and 
#   the practical way of finding the best number of nearby APs (fixed and 
#   flexible).
#
#   Author. @victorcroisfelt
#
#   Date. December 28, 2021
#
#   This code is part of the code package used to generate the numeric results
#   of the paper:
#
#   Croisfelt, V., Abrão, T., and Marinello, J. C., “User-Centric Perspective in
#   Random Access Cell-Free Aided by Spatial Separability”, arXiv e-prints, 2021.
#
#   Available on:
#
#                   https://arxiv.org/abs/2107.10294
#
#   Comment. You need to run:
#
#       - plot_fig07c_anaa_lower.py
#       - plot_fig07d_anaa_practical.py
#       - plot_fig08_tcp.py
#
#   to actually plot the figure using the data generated by this script.
#   Please, make sure that you have the files produced by:
#       
#       - lookup_fig07_08_delta.py
#       - lookup_fig07_08_Lmax_lower.py
#       - lookup_fig07_08_Lmax_practical.py
#
########################################
import numpy as np

import time

from settings_fig07_08 import *

########################################
# SELECTION
########################################

# Choose estimator
#estimator = 'est1'
estimator = 'est2'
#estimator = 'est3'

# Choose method of selection of number of nearby APs (Section V.E-1))
method = 'fixed'
#method = 'flexible'

# Choose bound
#bound = 'lower'
bound = 'practical' 

########################################
# Lookup tables
########################################

# Get lookup tables from the according files
if bound == 'lower':

    load = np.load("lookup/lookup_fig07_08_Lmax_" + estimator + "_lower.npz", allow_pickle=True)
    Lmax_lookup = load["best_Lmax"]
    Lmax_lookup = Lmax_lookup.item()

    # Load possible values of delta for Estimator 3
    if estimator == "est3":

        load = np.load("lookup/lookup_fig07_08_delta.npz", allow_pickle=True)
        delta_lookup = load["delta"]
        delta_lookup = delta_lookup.item()

if bound == 'practical':

    load = np.load("lookup/lookup_fig07_08_practical.npz", allow_pickle=True)
    Lmax_lookup = load["Lmax_practical"]

    # Load possible values of delta for Estimator 3
    if estimator == "est3":
        delta_lookup = load["delta_practical"]

########################################
# Simulation
########################################
print("--------------------------------------------------")
print("Data Figs 07 & 08: CF-SUCRe -- ANAA and TCP")
print("\t Estimator: " + estimator)
print("\t Method: " + method)
print("\t Bound: " + bound) 
print("--------------------------------------------------")

# Store total time 
total_time = time.time()

# Store enumeration of L
enumerationL = np.arange(L)

# Prepare to save final waiting times
finalWaitingTimes = np.zeros((K0values.size, maxAttempts))

# Prepare to save average number of operative APs
avg_operativeAPs = np.empty((K0values.size, numRAblocks))
avg_operativeAPs[:] = np.nan

# Prepare to save average number of active pilots per operative AP
avg_activePilots_perAP = np.empty((K0values.size, numRAblocks))
avg_activePilots_perAP[:] = np.nan

# More one thing...
if estimator == "est3":

    # Prepare to save average efficient power per operative AP
    avg_effpwr_perAP = np.empty((K0values.size, numRAblocks))
    avg_effpwr_perAP[:] = np.nan

#####


# Go through all different number of inactive UEs
for kk, K0 in enumerate(K0values):

    # Storing time
    timer_start = time.time()

    # Print current data point
    print(f"\tinactive UEs: {kk}/{len(K0values)-1}")

    # Generate the number of UEs that wish to access the network (for the first
    # time) in each of the RA blocks
    newUEs = np.random.binomial(K0, pA, size=numRAblocks)

    # Extract current Lmax
    if bound == 'practical':
        Lmax = Lmax_lookup[kk]

    # Initiate memories to store set of UEs that have failed to access the
    # network
    waitingTime = [] # Contains number of access attempts
    waitingBetas = [] # Average channel gains of UEs

    # Go through all RA blocks
    for r in range(numRAblocks):


        # Generate noise realizations at APs
        n_ = np.sqrt(sigma2/2)*(np.random.randn(N, L, taup) + 1j*np.random.randn(N, L, taup))


        #####
        # Generating UEs
        #####


        # Generate UEs locations
        newUElocations = squareLength*(np.random.rand(newUEs[r]) + 1j*np.random.rand(newUEs[r]))

        # Compute UEs distances to each AP
        newUEdistances = abs(APpositions - newUElocations[:, None])

        # Compute average channel gains according to Eq. (1)
        newBetas = 10**((94.0 - 30.5 - 36.7 * np.log10(np.sqrt(newUEdistances**2 + 10**2)))/10)


        #####
        # Preparing for RA Attempt
        #####


        # Combine the new UEs with the ones that have made previous access
        # attempts
        if len(np.array(waitingBetas)) == 0:
            betas = newBetas
        else:
            betas = np.concatenate((newBetas, np.array(waitingBetas)))

        # Compute number of UEs that will send pilots
        numberOfAccessingUEs = betas.shape[0]

        # Randomize if each of the UEs that retransmit pilots should really send
        # a pilot in this RA block. One means retransmit and zero means do not
        # retransmit in this block
        shouldWaitingUsersRetransmit = np.random.binomial(1, tryAgainProb, size=len(waitingTime)).astype(int)

        # Create a list of the UEs that will send pilots (all new UEs transmit
        # pilots)
        accessAttempt = np.concatenate((np.ones(newUEs[r], dtype=np.uint), shouldWaitingUsersRetransmit)).astype(int)

        # Randomize which pilot each UE chose
        pilotSelections = accessAttempt*np.random.randint(1, taup+1, size=numberOfAccessingUEs).astype(int)
        pilotSelections += -1

        # Count the number of pilots that each of the UEs will have transmitted,
        # after this block
        accessAttempts = np.concatenate((np.ones(newUEs[r], dtype=np.uint), waitingTime + shouldWaitingUsersRetransmit)).astype(int)


        # Check existence of transmission
        if len(accessAttempts) != 0:

            # Generate channel matrix at each AP equipped with N antennas
            G = np.sqrt(betas[None, :, :]/2)*(np.random.randn(N, numberOfAccessingUEs, L) + 1j*np.random.randn(N, numberOfAccessingUEs, L))

            # Prepare a list of UEs that succeed in the random access
            successfulAttempt = np.zeros(numberOfAccessingUEs, dtype=bool)

            # Get list of active pilots
            activePilots = np.unique(pilotSelections)

            
            # Prepare to save local results
            list_of_operative_APs = []
            list_of_active_pilots_per_opAP = np.zeros((L))


            if estimator == "est3":
                effpwr_per_pilot = np.empty((taup))
                effpwr_per_pilot[:] = np.nan

            

            # Go through all active RA pilots
            for t in activePilots:

                # Extract UEs that transmit pilot t
                UEindices = np.where(pilotSelections == t)[0]

                # Obtain collision size
                collisionSize = len(UEindices)

                # Obtain best Lmax
                if bound == 'lower':
                    Lmax = Lmax_lookup.get(collisionSize, L)

                # Compute received signal according to Eq. (4)
                Yt = np.sqrt(p * taup)*(G[:, UEindices, :]).sum(axis=1) + n_[:, :, t]

                # Store l2-norms of Yt
                Yt_norms = np.linalg.norm(Yt, axis=0)

                # Obtain pilot activity matrix according to Eq. (8)
                Atilde = (1/N) * Yt_norms**2
                Atilde[Atilde < sigma2] = 0.0

                # Obtain set of pilot-serving APs (Definition 2)
                Pcal = np.argsort(Atilde)[-Lmax:]
                Pcal = np.delete(Pcal, Atilde[Pcal] == 0)

                # Store local results
                list_of_operative_APs.append(set(Pcal))
                list_of_active_pilots_per_opAP[Pcal] += 1


                #####
                # SUCRe: step 2
                #####


                # Check estimator
                if estimator == 'est3':

                    # Denominator according to Eq. (34) and (35)
                    den = np.sqrt(N *  np.maximum(Atilde - sigma2, np.zeros((Atilde.shape))).sum())

                    # Compute precoded DL signal according to Eq. (35)
                    Vt = np.sqrt(ql) * (Yt[:, Pcal] / den)

                    # Store effective power per pilot
                    effpwr_per_pilot[t] = ql * ((Yt_norms[Pcal]**2).mean()) / (den**2)

                else:

                    # Compute precoded DL signal according to Eq. (10)
                    Vt = np.sqrt(ql) * (Yt[:, Pcal] / Yt_norms[Pcal])


                #####
                # SUCRe: step 3
                #####


                # Prepare a list of UEs that decide to retransmit pilot t
                retransmit = np.zeros((collisionSize), dtype=bool)

                # Prepare to save a list of all the natural sets of nearby APs 
                # from colliding UEs
                checkCcal = []


                # Go through all colliding UEs
                for k in range(collisionSize):

                    # Compute received DL signal at UE k according to Eq. (12)
                    noise = np.sqrt(sigma2/2)*(np.random.randn() + 1j*np.random.randn())
                    z_k = np.sqrt(taup) * (G[:, UEindices[k], Pcal].conj() * Vt).sum() + noise

                    # Obtain natural set of nearby APs of UE k (Definition 1)
                    checkCcal_k = enumerationL[ql*betas[UEindices[k], :] > sigma2]


                    if len(checkCcal_k) == 0:
                        checkCcal_k = np.array([np.argmax(ql * betas[UEindices[k], :])])

                    # Store checkCcal_k
                    checkCcal.append(set(checkCcal_k))


                    # Check method to determine size of Ccal
                    if method == 'fixed':


                        #####
                        # Estimation
                        #####


                        # Compute constants
                        cte = z_k.real/np.sqrt(N)
                        num = np.sqrt(ql * p) * taup * betas[UEindices[k], checkCcal_k]

                        if estimator == 'est1':

                            # Compute estimate according to Eq. (28)
                            alphahat = ((num.sum()/cte)**2) - sigma2

                        elif estimator == 'est2':

                            num23 = num**(2/3)
                            cte2 = (num23.sum()/cte)**2

                            # Compute estimate according to Eq. (32)
                            alphahat = (cte2 * num23 - sigma2).sum()

                        elif estimator == 'est3':

                            # Define compensation factor in Eq. (39)
                            if bound == 'lower':
                                delta = delta_lookup.get((collisionSize, Lmax), delta_lookup.get((50, Lmax)))
                            else:
                                delta = delta_lookup[kk]

                            # Compute new constant according to Eq. (38)
                            underline_cte = delta * (z_k.real - sigma2)/np.sqrt(N)

                            # Compute estimate according to Eq. (40)
                            alphahat = (num.sum() / underline_cte)**2

                        # Compute own total UL signal power in Eq. (15)
                        gamma = p * taup * betas[UEindices[k], checkCcal_k].sum()

                        # Avoiding underestimation
                        if alphahat < gamma:
                            alphahat = gamma


                        #####
                        # Decision
                        #####


                        # Apply the retransmission decision rule
                        retransmit[k] = gamma > alphahat/2


                    elif method == 'flexible':


                        # Compute unchanged constant
                        if estimator == 'est3':

                            # Define compensation factor in Eq. (39)
                            if bound == 'lower':
                                delta = delta_lookup.get((collisionSize, Lmax), delta_lookup.get((50, Lmax)))
                            else:
                                delta = delta_lookup[kk]

                            # Compute new constant according to Eq. (38)
                            cte = delta * (z_k.real - sigma2)/np.sqrt(N)

                        else:

                            cte = z_k.real/np.sqrt(N)

                        # Define aux variable
                        aux = int(checkCcal_k.size)

                        while aux > 0:

                            # Obtain current Ccal_k
                            current_Ccal_k = checkCcal_k[-aux:]


                            #####
                            # Estimation
                            #####


                            # Compute new numerator
                            num = np.sqrt(ql * p) * taup * betas[UEindices[k], current_Ccal_k]

                            if estimator == 'est1':

                                # Compute estimate according to Eq. (28)
                                alphahat = ((num.sum()/cte)**2) - sigma2

                            elif estimator == 'est2':

                                num23 = num**(2/3)
                                cte2 = (num23.sum()/cte)**2

                                # Compute estimate according to Eq. (32)
                                alphahat = (cte2 * num23 - sigma2).sum()

                            elif estimator == 'est3':

                                # Compute estimate according to Eq. (40)
                                alphahat = (num.sum() / cte)**2

                            # Compute own total UL signal power in Eq. (15)
                            gamma = p * taup * betas[UEindices[k], current_Ccal_k].sum()

                            # Avoiding underestimation
                            if alphahat < gamma:
                                alphahat = gamma

                            #####
                            # Decision
                            #####

                            # Evaluate decision
                            if gamma > alphahat/2:
                                retransmit[k] = 1
                                break

                            # Update aux variable
                            aux -= 1


                #####
                # Evaluate succesful attempts (!!!)
                #####


                # Single transmission
                if sum(retransmit) == 1:
                    successfulAttempt[UEindices[retransmit]] = True
                    finalWaitingTimes[kk, int(accessAttempts[UEindices[retransmit]] - 1)] += 1

                # Collision case
                elif sum(retransmit) > 1:

                    # Go through all colliding UEs
                    for k in range(collisionSize):

                        # Extract current natural set of nearby APs of UE k
                        checkCcal_k = checkCcal[k]

                        # Compute intersection with Pcal
                        int_Pcal_checkCcal_k = set.intersection(set(Pcal), checkCcal_k)

                        # Spatial Separability (Definition 3): the following ifs
                        # are performed according to the defition of spatial
                        # separability

                        # Check first condition
                        if len(int_Pcal_checkCcal_k) > 0:

                            # Obtain union of natural sets excluding the one from
                            # UE k
                            checkCcal__k = set.union(*(checkCcal[:k] + checkCcal[k+1:]))

                            # Compute intersection with Pcal
                            int_Pcal_checkCcal__k = set.intersection(set(Pcal), checkCcal__k)

                            # Check second condition
                            if len(int_Pcal_checkCcal_k - int_Pcal_checkCcal__k) != 0:

                                # Succesful attempt: UE k is spatially separable
                                successfulAttempt[UEindices[k]] = True
                                finalWaitingTimes[kk, int(accessAttempts[UEindices[k]] - 1)] += 1

            # Store results

            # Operative APs
            if len(list_of_operative_APs) != 0:
                avg_operativeAPs[kk, r] = len(set.union(*list_of_operative_APs))

            # Active pilots
            list_of_active_pilots_per_opAP[list_of_active_pilots_per_opAP == 0.0] = np.nan
            avg_activePilots_perAP[kk, r] = np.nanmean(list_of_active_pilots_per_opAP)

            # Effective power
            if estimator == "est3":
                avg_effpwr_perAP[kk, r] = np.nanmean(effpwr_per_pilot)

            # Determine which of the UEs have failed too many times with their 
            # access attempts and will give up
            giveUp = (accessAttempts[successfulAttempt == 0] == maxAttempts)
            finalWaitingTimes[kk, -1] += sum(giveUp)

            # Keep the important parameters for all the UEs that failed to
            # access the network and did not give up
            mask_remaining = np.logical_and((successfulAttempt == 0), (accessAttempts < maxAttempts))

            waitingTime = accessAttempts[mask_remaining]
            waitingBetas = betas[mask_remaining]

    print("\t[|U|] elapsed " + str(np.round(time.time() - timer_start, 4)) + " seconds.\n")

# Compute average number of operative APs
final_avg_operativeAPs = np.nanmean(avg_operativeAPs, axis=-1)

# Compute average number of active pilots
final_avg_activePilots_perAP = np.nanmean(avg_activePilots_perAP, axis=-1)

# Compute ANAA
anaa = (np.arange(1, maxAttempts+1)[np.newaxis, :] * (finalWaitingTimes/np.sum(finalWaitingTimes, axis=-1)[:, np.newaxis])).sum(axis=-1)

# Compute TCP in Eq. (43)
if estimator == "est3":

    # Average effective power
    final_effpwr = np.nanmean(avg_effpwr_perAP, axis=-1)

    # Compute TCP
    tcp = anaa * (1 + taup) * final_effpwr * final_avg_activePilots_perAP * final_avg_operativeAPs

else:

    # Compute TCP
    tcp = anaa * (1 + taup) * ql * final_avg_activePilots_perAP * final_avg_operativeAPs

print("total simulation time was " + str(np.round(time.time() - total_time, 4)) + " seconds.\n")

print("wait for data saving...\n")

# Save simulation results
np.savez('data/fig07_08_cellfree_' + estimator + '_' + method + '_' + bound + '.npz',
    K0values=K0values,
    anaa=anaa,
    tcp=tcp
)

print("your data has been saved in the /data folder.\n")

print("------------------- all done :) ------------------")