data_fig07_08_cellular.py

########################################
#   data_fig07_08_cellular.py
#
#   Description. Script used to generate data for Figs. 7 and 8 of the paper
#   concerning cellular curves (Ce-SUCRe).
#
#   Author. @victorcroisfelt
#
#   Date. December 27, 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.
########################################
import numpy as np

import time

from settings_fig07_08 import *

########################################
# Simulation
########################################
print("--------------------------------------------------")
print("Data Figs 07 & 08: CeSUCRe -- ANAA and TCP")
print("--------------------------------------------------\n")

# Store total time 
total_time = time.time()

# Prepare to save simulation results
finalWaitingTimes = np.zeros((K0values.size, maxAttempts))
avg_activePilots = np.zeros((K0values.size))


#####


# Generate noise realizations at the BS
n_ = np.sqrt(sigma2/2)*(np.random.randn(M, taup, numRAblocks) + 1j*np.random.randn(M, taup, numRAblocks))


# 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}/{K0values.size-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)

    # 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):


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


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

        # Compute UEs distances to BS
        newUEdistances = abs(BSposition - newUElocations)

        # 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
        betas = np.concatenate((newBetas, np.array(waitingBetas)))

        # Compute number of UEs that will send pilots
        numberOfAccessingUEs = len(betas)

        # 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))

        # 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))

        # Randomize which pilot each UE chose
        pilotSelections = accessAttempt*np.random.randint(1, taup+1, size=numberOfAccessingUEs)
        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))

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

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

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

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

            # Update number of active pilots
            avg_activePilots[kk] += len(activePilots)

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

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

                # Obtain collision size
                collisionSize = len(UEindices)

                # Compute received signal (equivalent to Eq. (4))
                yt = np.sqrt(p * taup) * (G[:, UEindices]).sum(axis=1) + n_[:, tt, r]

                # Compute precoded DL signal (equiavalent to Eq. (10))
                v_ = np.sqrt(q) * (yt / np.linalg.norm(yt))

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

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

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


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


                    # Compute constants
                    cte = z_k.real/np.sqrt(M)
                    num = np.sqrt(q * p) * taup * betas[UEindices[k]]

                    # Compute estimate
                    alphahat = ((num/cte)**2) - sigma2

                    # Compute own total UL signal power (equivalent to Eq. 
                    #(15))
                    gamma = p * taup * betas[UEindices[k]]

                    # Avoiding underestimatio
                    if alphahat < gamma:
                        alphaest = gamma

                    # Apply retransmission decision rule -- Eq. (17)
                    retransmit[k] = gamma > alphahat/2

                # Check if only one UE has decided to retransmit pilot t and
                # then admit the UE for data transmission and store the 
                # number of access attempts that the UE made
                if sum(retransmit) == 1:
                    successfulAttempt[UEindices[retransmit]] = True
                    finalWaitingTimes[kk, int(accessAttempts[UEindices[retransmit]] - 1)] += 1

            # 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 active pilots
avg_activePilots *= 1/numRAblocks

# 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)
tcp = anaa * (1 + taup) * q * avg_activePilots

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_cellular.npz',
    K0values=K0values,
    anaa=anaa,
    tcp=tcp
)

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

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