text: ########################################
  #   data_fig05_barplot_cellfree.py
  #
  #   Description. Script used to generate data for Fig. 5 of the paper regarding
  #   the cell-free 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.
  #   
  #   Please, make sure that you have the files produced by:
  #       
  #       - lookup_fig05_06_delta.py
  #       - lookup_fig05_06_best_pair.py
  #
  ########################################
  import numpy as np
  
  import time
  
  ########################################
  # Preamble
  ########################################
  np.random.seed(42)
  
  ########################################
  # System parameters
  ########################################
  
  # Define number of APs
  L = 64
  
  # Define number of antennas per AP
  N = 8
  
  # UL transmit power
  p = 100
  
  # DL transmit power per AP
  ql = 200/L
  
  # Define noise power
  sigma2 = 1
  
  # Number of RA pilot signals
  taup = 5
  
  ########################################
  # SELECTION
  ########################################
  
  # Choose the estimator
  estimator = "est1"
  estimator = "est2"
  estimator = "est3"
  
  ########################################
  # Lookup table
  ########################################
  
  # Load best pair look up table
  load = np.load("lookup/lookup_fig05_best_pair_" + estimator + ".npz", allow_pickle=True)
  best_pair_lookup = load["best_pair"]
  best_pair_lookup = best_pair_lookup.item()
  
  # Load possible values of delta for Estimator 3
  if estimator == "est3":
  
      load = np.load("lookup/lookup_fig05_06_delta.npz", allow_pickle=True)
      delta_lookup = load["delta"]
      delta_lookup = delta_lookup.item()
  
  ########################################
  # Geometry
  ########################################
  
  # Define square length
  squareLength = 400
  
  # Create square grid of APs
  APperdim = int(np.sqrt(L))
  APpositions = np.linspace(squareLength/APperdim, squareLength, APperdim) - squareLength/APperdim/2
  APpositions = APpositions + 1j*APpositions[:, None]
  APpositions = APpositions.reshape(L)
  
  ########################################
  # 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 -- cell-free")
  print("\testimator: " + estimator)
  print("\tN = " + str(N))
  print("--------------------------------------------------\n")
  
  # Store total time
  total_time = time.time()
  
  # Store enumeration of L
  enumerationL = np.arange(L)
  
  # Prepare to save NMSE stats
  nmse = np.zeros((3, collisions.size))
  
  
  #####
  
  
  # Generate noise realizations at APs
  n_ = np.sqrt(sigma2/2)*(np.random.randn(numsetups, N, L, numchannel) + 1j*np.random.randn(numsetups, N, L, 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
      UElocations = squareLength*(np.random.rand(numsetups, collisionSize) + 1j*np.random.rand(numsetups, collisionSize))
  
      # Compute UEs distances to each AP
      UEdistances = np.abs(UElocations[:, :, np.newaxis] - APpositions)
  
      # 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 each AP equipped with N antennas
      Gnorm_ = np.sqrt(1/2)*(np.random.randn(numsetups, N, collisionSize, L, numchannel) + 1j*np.random.randn(numsetups, N, collisionSize, L, numchannel))
  
      # Compute channel matrix
      G_ = np.sqrt(channel_gains[:, None, :, :, None]) * Gnorm_
  
      # Compute received signal according to Eq. (4)
      Yt_ = np.sqrt(p * taup) * G_.sum(axis=2) +  n_
  
      # Store l2-norms of Yt
      Yt_norms = np.linalg.norm(Yt_, axis=1)
  
      # Obtain pilot activity vector according to Eq. (8)
      atilde_t = (1/N) * Yt_norms**2
      atilde_t[atilde_t < sigma2] = 0.0
  
      # Extract current Csize and Lmax 
      (Csize, Lmax) = best_pair_lookup[(collisionSize, N)]
  
      # Prepare to save inner NMSE
      nmse_in = np.zeros((numsetups, collisionSize, numchannel))
  
  
      # Go through all setups
      for ss in range(numsetups):
  
          # Go through all channel realizations
          for ch in range(numchannel):
  
              # Obtain set of pilot-serving APs (Definition 2)
              Pcal = np.argsort(atilde_t[ss, :, ch])[-Lmax:]
              Pcal = np.delete(Pcal, atilde_t[ss, Pcal, ch] == 0)
  
  
              #####
              # SUCRe - step 2
              #####
  
  
              if estimator == 'est3':
  
                  # Denominator according to Eqs. (34) and (35)
                  den = np.sqrt(N * (atilde_t[ss, :, ch] - sigma2).sum())
  
                  # Compute precoded DL signal according to Eq. (35)
                  Vt_ = np.sqrt(ql) * (Yt_[ss][:, Pcal, ch] / den)
  
              else:
  
                  # Compute precoded DL signal according to Eq. (10)
                  Vt_ = np.sqrt(ql) * (Yt_[ss][:, Pcal, ch] / Yt_norms[ss, Pcal, ch][None, :])
  
              # Compute true total UL signal power of colliding UEs 
              # according to Eq. (16)
              alpha_true = p * taup * channel_gains[ss, :, Pcal].sum()
  
              # Go through all colliding UEs
              for k in range(collisionSize):
  
                  # Compute received DL signal at UE k according to Eq. 
                  # (12)
                  z_k = np.sqrt(taup) * (G_[ss][:, k, Pcal, ch].conj() * Vt_).sum() + eta[ss, k, ch]
  
                  # Obtain set of nearby APs of UE k (Definition 1)
                  Ccal_k = np.argsort(ql * channel_gains[ss, k, :])[-Csize:]
  
                  # Obtain natural set of nearby APs of UE k (Definition 1)
                  checkCcal_k = enumerationL[ql * channel_gains[ss, k, :] > sigma2]
  
                  if len(checkCcal_k) == 0:
                      checkCcal_k = np.array([np.argmax(ql * channel_gains[ss, k, :])])
  
                  if len(Ccal_k) > len(checkCcal_k):
                      Ccal_k = checkCcal_k
  
  
                  #####
                  # Estimation
                  #####
  
  
                  # Compute constants
                  cte = z_k.real/np.sqrt(N)
                  num = np.sqrt(ql * p) * taup * channel_gains[ss, k, 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':
  
                      # Define compensation factor in Eq. (39)
                      delta = delta_lookup[(collisionSize, N, Lmax)]
  
                      # 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 * channel_gains[ss, k, Ccal_k].sum()
  
                  # Avoiding underestimation
                  if alphahat < gamma:
                      alphahat = gamma
  
                  # Get and store inner loop stats
                  nmse_in[ss, k, ch] = (np.abs(alphahat - alpha_true)**2)/(alpha_true**2)
  
      # Average out channel realizations
      nmse_in = nmse_in.mean(axis=-1)
  
      # Save NMSE 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_cellfree_' + estimator + '.npz',
      nmse=nmse
  )
  
  print("the data has been saved in the /data folder.\n")
  
  print("------------------- all done :) ------------------")