data_fig05_barplot_cellular.py

########################################
#   data_fig05_barplot_cellular.py
#
#   Description. Script used to generate data for Fig. 5 of the paper regarding
#   the cellular curves. You should choose the estimator among the three avai-
#   lable.
#
#   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_fig05_barplot.py
#
#   to actually plot the figure using the data generated by this script.
########################################
import numpy as np

import time

########################################
# Preamble
########################################
np.random.seed(42)

########################################
# System parameters
########################################

# Define number of BS antennas
M = 64

# UL transmit power
p = 100

# DL transmit power
q = 200

# Define noise power
sigma2 = 1

# Number of RA pilot signals
taup = 5

########################################
# Geometry
########################################

# Define square length
squareLength = 400

# Define BS position
BSposition = (squareLength/2)*(1 + 1j)

########################################
# Simulation parameters
########################################

# Set the number of setups
numsetups = 100

# Set the number of channel realizations
numchannel = 100

# Range of collision sizes
collisions = np.arange(1, 11)

########################################
# Simulation
########################################
print("--------------------------------------------------")
print("Data Fig 05: barplot -- cellular")
print("--------------------------------------------------\n")

# Store total time
total_time = time.time()

# Prepare to save NMSE stats
nmse = np.zeros((3, collisions.size))


#####


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

# Generate noise realization at UEs
eta = np.sqrt(sigma2/2)*(np.random.randn(numsetups, collisions.max(), numchannel) + 1j*np.random.randn(numsetups, collisions.max(), numchannel))


# Go through all collision sizes
for cs, collisionSize in enumerate(collisions):

    # Storing time
    timer_start = time.time()

    # Print current data point
    print(f"\tcollision: {cs}/{collisions.size-1}")


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

    # Generate UEs locations
    UEpositions = squareLength*(np.random.rand(numsetups, collisionSize) + 1j*np.random.rand(numsetups, collisionSize))

    # Compute UEs distances to the BS
    UEdistances = np.abs(BSposition - UEpositions)

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

    # Generate normalized channel matrix for the BS equipped with M antennas
    Gnorm_ = np.sqrt(1/2)*(np.random.randn(numsetups, M, collisionSize, numchannel) + 1j*np.random.randn(numsetups, M, collisionSize, numchannel))

    # Compute channel matrix
    G_ = np.sqrt(channel_gains[:, None, :, None]) * Gnorm_

    # Compute received signal (equivalent to Eq. (4))
    Yt_ = np.sqrt(p * taup) * G_.sum(axis=2) +  n_

    # Compute precoded DL signal (equiavalent to Eq. (10))
    Vt_ = np.sqrt(q) * (Yt_ / np.linalg.norm(Yt_, axis=1)[:, None, :])

    # Compute true alpha
    alpha_true = p * taup * channel_gains.sum(axis=-1)

    # Prepare to save inner NMSE
    nmse_in = np.zeros((numsetups, collisionSize))

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

        # Compute received DL signal at UE k (equivalent to Eq. (12))
        z_k = np.sqrt(taup) * (G_[:, :, k, :].conj() * Vt_).sum(axis=1) + eta[:, k, :]


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


        # Compute constants
        den = z_k.real/np.sqrt(M)
        num = np.sqrt(q * p) * taup * channel_gains[:, k]

        # Compute estimate
        alphahat = ((num[:, None]/den)**2) - sigma2

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

        # Avoiding underestimation
        for ch in range(numchannel):

            mask = alphahat[:, ch] <= gamma
            alphahat[mask, ch] = gamma[mask]

        # Compute stats
        nmse_in[:, k] = np.mean((np.abs(alphahat - alpha_true[:, None])**2), axis=-1)/(alpha_true**2)

    # Save outer loop stats
    nmse[:, cs] = np.stack((np.percentile(nmse_in, 25), np.median(nmse_in), np.percentile(nmse_in, 75)))

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

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/fig05_barplot_cellular.npz',
    nmse=nmse
)

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

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