Commit 9ecc180b by JRock007

### Translate Python scripts to English

parent 510e8e95
 import math import matplotlib.pyplot as plt import numpy as np import Modele_N_recepteurs as NOMAStat import Theorie_N_recepteurs as NOMATh def compLogLog(xlog=True, ylog=True): """ Graph le Bit Error Rate en fonction de P1/sigma**2 """ N = 2 Pmax = 10000 g = 1 it = 2000 P2 = 0.8 * Pmax P1 = Pmax - P2 P = [P1, P2] Erreurs = [] Probas = [] Sigmas = np.linspace(math.sqrt(P1 / 50), math.sqrt(P1), 50) for sigma in Sigmas: probas = NOMATh.theorie(g, sigma, P, N) erreurs = NOMAStat.stats(g, sigma, P, N, it) Probas.append(probas[0]) Erreurs.append(erreurs[0]) fig = plt.figure() ax = fig.add_subplot(1, 1, 1) if xlog: ax.set_xscale("log") plt.xlabel("RSB (log)") else: plt.xlabel("RSB") if ylog: ax.set_yscale("log") plt.ylabel("BER (log)") else: plt.ylabel("BER") RSB = P1 / np.array(np.square(Sigmas)) ax.plot(RSB, Probas, label="Taux d'erreur théorique") ax.plot(RSB, Erreurs, label="Taux d'erreur statistique", marker="x", ls="") plt.legend() plt.show() def compRepartitionPuissance(): N = 2 sigma = 0.25 / math.sqrt(2) # https://www.gnuradio.org/doc/doxygen/classgr_1_1analog_1_1noise__source__c.html#a6f77085e298f93b770ab0693342d1dbb g = 1 Pmax = 1 it = 2048 Erreurs = [] Probas = [] Puissances = np.linspace(Pmax / 2, 0.99 * Pmax, 50) for P2 in Puissances: P = [Pmax - P2, P2] probas = NOMATh.theorie(g, sigma, P, N) erreurs = NOMAStat.stats(g, sigma, P, N, it) Probas.append(sum(probas)) Erreurs.append(sum(erreurs)) fig = plt.figure() ax = fig.add_subplot(1, 1, 1) plt.xlabel("P2/Pmax") plt.ylabel("BER total") P = Puissances / Pmax ax.plot(P, Probas, label="Taux d'erreur théorique") ax.plot(P, Erreurs, label="Taux d'erreur statistique", marker="x", ls="") plt.legend() plt.show() def compDiffCarres(loop_range=5, it=1000): """ Mesure la différence entre le modèle théorique et le modèle statistique """ print("Étude comparative de la théorie et du modèle") for i in range(loop_range): # Build a new model and calculate therory N = np.random.randint(2, 6) sigma = np.sqrt(np.random.randint(0, 50) / 10) g = np.random.rand() P = [np.random.uniform(1, 10)**k for k in range(N)] print("**********************************") print("N =", N) print("Puissances =", P) print("Atténuation g =", g) print("Bruit de variance σ² =", sigma**2) print("Simulation de n = {} itérations".format(it)) erreurs = NOMAStat.stats(g, sigma, P, N, it) probas = NOMATh.theorie(g, sigma, P, N) diff = [(erreurs[i] - probas[i])**2 for i in range(N)] print(diff) def plotTh2users(P, sigma, g): probas = [] puissances = np.array(P) for p in puissances: probas.append(NOMATh.theorie(g, sigma, p, 2)[0]) plt.plot(puissances[:, -1], probas) plt.show() if __name__ == '__main__': # compRepartitionPuissance() P1 = np.arange(0.88, 0.999, 0.001) P = [(1 - p, p) for p in P1] sigma = 0.075 g = 1 plotTh2users(P, sigma, g)
 ... ... @@ -2,15 +2,16 @@ import math import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import numpy as np import Theorie_N_recepteurs as NOMATh import Theory_N_receivers as NOMATh # Experimental parameters N = 3 Pmax = 1 sigma = 0.05 g = 1 it = 10000 # Experimental measures P1 = [0., 0., 0., 0.05, 0.1, 0., 0.05, 0.1, 0.15, 0., 0.05, 0.1, 0.15, 0.2, 0., 0.05, 0.1, 0.15, 0.2, 0.25, 0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0., 0.05, 0.1, 0.15, 0.2, 0., 0.05, 0.1, 0.] P2 = [0., 0.05, 0.1, 0.1, 0.1, 0.15, 0.15, 0.15, 0.15, 0.2, 0.2, 0.2, 0.2, 0.2, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.4, 0.4, 0.4, 0.4, 0.4, 0.45, 0.45, 0.45, 0.5] BER3 = [0, 0, 0, 0, 0, 0, 0, 0.004, 0.045, 0, 0.0005, 0.035, 0.16, 0.24, 0, 0.01, 0.14, 0.23, 0.24, 0.25, 0, 0.075, 0.22, 0.25, 0.25, 0.25, 0.25, 0, 0.19, 0.245, 0.26, 0.25, 0.25, 0.26, 0.02, 0.25, 0.25, 0.255, 0.25, 0.05, 0.24, 0.25, 0.25] ... ... @@ -20,30 +21,39 @@ BER1 = [0.5, 0.5, 0.5, 0.05, 0.25, 0.5, 0.023, 0.1, 0.3, 0.5, 0.011, 0.07, 0.3, def main(): """ Graph le BER théorique et le BER expérimental Graph the theoritical error probability and experimental BER """ # List of the sum of BERs for each parameter value BERTotalList = [] BERTotalThList = [] # List of the BER for 3rd user (highest power) for each parameter value BER3List = [] BER3ThList = [] # List of the BER for 2nd user for each parameter value BER2List = [] BER2ThList = [] # List of the BER for 1rst user (lowest power) for each parameter value BER1List = [] BER1ThList = [] for i in range(len(P1)): # Extract power frome experimental results p1 = P1[i] p2 = P2[i] p3 = Pmax - p2 - p1 # Same with BER ber3 = BER3[i] ber2 = BER2[i] ber1 = BER1[i] ber1th, ber2th, ber3th = NOMATh.theorie(g, sigma, [p1, p2, p3], N) # Now, computer the theoritical BER given the experimental conditions ber1th, ber2th, ber3th = NOMATh.theory(g, sigma, [p1, p2, p3], N) # Append all those values to the lists BERTotalList.append(ber1 + ber2 + ber3) BERTotalThList.append(ber1th + ber2th + ber3th) ... ... @@ -56,12 +66,13 @@ def main(): BER1List.append(ber1) BER1ThList.append(ber1th) # Plot in 3D the values of BER plt.figure() ax = plt.axes(projection="3d") ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("Somme des BER") ax.set_title("BER total à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = {} et σ = {}".format(N, Pmax, g, sigma)) ax.set_zlabel("Sum of BER") ax.set_title("Evolution of the total BER for N={} users, depending on the power, with Pmax = {}, g = {} and σ = {}".format(N, Pmax, g, sigma)) ax.plot(P1, P2, BERTotalThList, markerfacecolor="r", marker="o", markersize=5, linestyle="None") ax.plot(P1, P2, BERTotalList, markerfacecolor="b", marker="x", markersize=5, linestyle="None") plt.show(True) ... ... @@ -71,7 +82,7 @@ def main(): ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("BER3") ax.set_title("BER de l'utilisateur 3 à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = {} et σ = {}".format(N, Pmax, g, sigma)) ax.set_title("Evolution of the 3rd user's BER for N={} users, depending on the power, with Pmax = {}, g = {} and σ = {}".format(N, Pmax, g, sigma)) ax.plot(P1, P2, BER3List, markerfacecolor="r", marker="o", markersize=5, linestyle="None") ax.plot(P1, P2, BER3ThList, markerfacecolor="b", marker="x", markersize=5, linestyle="None") plt.show(True) ... ... @@ -81,7 +92,7 @@ def main(): ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("BER2") ax.set_title("BER de l'utilisateur 2 à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = {} et σ = {}".format(N, Pmax, g, sigma)) ax.set_title("Evolution of the 2nd user's BER for N={} users, depending on the power, with Pmax = {}, g = {} and σ = {}".format(N, Pmax, g, sigma)) ax.plot(P1, P2, BER2List, markerfacecolor="r", marker="o", markersize=5, linestyle="None") ax.plot(P1, P2, BER2ThList, markerfacecolor="b", marker="x", markersize=5, linestyle="None") plt.show(True) ... ... @@ -91,10 +102,11 @@ def main(): ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("BER1") ax.set_title("BER de l'utilisateur 1 à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = {} et σ = {}".format(N, Pmax, g, sigma)) ax.set_title("Evolution of the 1rst user's BER for N={} users, depending on the power, with Pmax = {}, g = {} and σ = {}".format(N, Pmax, g, sigma)) ax.plot(P1, P2, BER1List, markerfacecolor="r", marker="o", markersize=5, linestyle="None") ax.plot(P1, P2, BER1ThList, markerfacecolor="b", marker="x", markersize=5, linestyle="None") plt.show(True) if __name__ == '__main__': main()
 import math import matplotlib.pyplot as plt import numpy as np import Model_N_receivers as NOMAStat import Theory_N_receivers as NOMATh def compLogLog(N=2, n=1, xlog=True, ylog=True, it=10000): """ Graph of the Bit Error Rate depending on the noise value for a given user N: number of users n: user for which the SNR and BER are computed it: number of iterations for the statistical approach """ Pmax = 10000 g = 1 # List of each user's power # Each user takes up 85% of the remaining power (up to the last-but-one user) P = [] for i in range(N - 1): P = [0.85 * (Pmax - sum(P))] + P # The last uses takes up the whole remaining power P = [Pmax - sum(P)] + P # The power of the user for which the graph is displayed Pn = P[n - 1] # Find, depending on the value of sigma for the noise, the BER for the nth-user Errors = [] Probas = [] Sigmas = np.linspace(math.sqrt(Pn / 50), math.sqrt(Pn), 50) for sigma in Sigmas: probas = NOMATh.theory(g, sigma, P, N) errors = NOMAStat.stats(g, sigma, P, N, it) Probas.append(probas[n - 1]) Errors.append(errors[n - 1]) # Draw the graph fig = plt.figure() ax = fig.add_subplot(1, 1, 1) # Set wether the scales should be logarithmic or not if xlog: ax.set_xscale("log") plt.xlabel("SNR user {} (log)".format(n)) else: plt.xlabel("SNR user {}".format(n)) if ylog: ax.set_yscale("log") plt.ylabel("BER user {} (log)".format(n)) else: plt.ylabel("BER user {}".format(n)) # Convert the sigmas for the noise to an RSB value RSB = Pn / np.array(np.square(Sigmas)) ax.plot(RSB, Probas, label="Theoretical error rate") ax.plot(RSB, Errors, label="Monte-Carlo error rate", marker="x", ls="") plt.legend() plt.show() def compPowerRepartitions(N=2, sigma=0.25, g=1, Pmax=1, it=2000): """ Plot the total BER for different repartitions of power """ # Convert the passed sigma to a value comparable to the one in GNURadio # See https://www.gnuradio.org/doc/doxygen/classgr_1_1analog_1_1noise__source__c.html#a6f77085e298f93b770ab0693342d1dbb sigma /= math.sqrt(2) # Compute which values of power repartition should be tested Powers = np.linspace(Pmax / 2, 0.99 * Pmax, 50) # Compute the BER and theoritical error probability for each possibility Errors = [] Probas = [] for P2 in Powers: P = [Pmax - P2, P2] probas = NOMATh.theory(g, sigma, P, N) errors = NOMAStat.stats(g, sigma, P, N, it) Probas.append(sum(probas)) Errors.append(sum(errors)) # Plot the result fig = plt.figure() ax = fig.add_subplot(1, 1, 1) plt.xlabel("P2/Pmax") plt.ylabel("Sum of BER") P = Powers / Pmax ax.plot(P, Probas, label="Theoretical error probability") ax.plot(P, Errors, label="Statistical error rate", marker="x", ls="") plt.legend() plt.show() if __name__ == '__main__': # compLogLog(N=3, n=1) compPowerRepartitions()
 import Theorie_N_recepteurs as NOMATh import Theory_N_receivers as NOMATh import numpy as np import time from scipy.optimize import minimize ... ... @@ -31,55 +31,64 @@ https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html https://stackoverflow.com/a/52003654/3624018 """ N = 4 Pmax = 10 g = [1, 0.5, 0.25, 0.1] sigma = [0.1, 0.2, 0.3, 0.4] N = 3 Pmax = 10000 g = [1, 0.8, 0.65] sigma = [5, 10, 15] print("N = {}, Pmax = {}, g = {} and σ = {}".format(N, Pmax, g, sigma)) def graph_2D(n=100, block=False): """ Graph le min pour N = 2 Graph the optimum for N = 2 """ # On génère les puissances possibles et on calcule f(P) # Generate all possible power repartitions and compute f(P) # If P1 > Pmax / 2, then P2 < P1. This cannot be, so we stop at Pmax / 2 X = np.linspace(0, Pmax / 2, n) Y = [f([x, Pmax - x]) for x in X] # On affiche les points # Plot those values ax = plt.axes() ax.plot(X, Y) ax.set_xlabel("P1") ax.set_ylabel("Somme des BER") ax.set_title("BER total à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = [{:.3g}, {:.3g}] et σ = [{:.3g}, {:.3g}]".format(N, Pmax, g[0], g[1], sigma[0], sigma[1])) ax.set_ylabel("Sum of BER") ax.set_title("Evolution of the sum of BER for N={} users depending on the power, with Pmax = {}, g = [{:.3g}, {:.3g}] et σ = [{:.3g}, {:.3g}]".format(N, Pmax, g[0], g[1], sigma[0], sigma[1])) # On trouve le min et on l'affiche # Find the optimum opti = np.argmin(Y) P = [X[opti], Pmax - X[opti]] BERmin = Y[opti] BER1, BER2 = f2(P) print("Graphically, min was found for P1 = [{:.2g}, {:.2g}] with BER={:.3g}".format(P[0], P[1], BERmin)) # Display the min we found ax.plot([P[0]], [BERmin], markerfacecolor="r", marker="o", markersize=5) ax.text(P[0] - 0.2, BERmin + 0.02, "P=({:.2g}, {:.2g}), BER={:.2g}".format(P[0], P[1], BERmin), fontsize=12, horizontalalignment="left", verticalalignment="center") ax.text(P[0], BERmin - 0.01, "P=({:.2g}, {:.2g}), BER1={:.2g}, BER2={:.2g}".format(P[0], P[1], BER1, BER2), fontsize=11, horizontalalignment="center", verticalalignment="center") # Afficher la fenêtre sans bloquer l'interface # Display the output plt.show(block) # Renvoyer la valeur trouvée # Return the value return P, BERmin def graph_3D(step=Pmax / 100, block=False): """ Graph le min pour N = 3 Graph the optimum for N = 3 """ X = [] Y = [] Z = [] # On calcule f(P) pour les valeurs valides de P # Compute f(P) for all possible power repartitions # The constraints are Pi > 0 and Pi+1 < Pi < Pi-1 # If P2 > Pmax / 2, then P3 < P2, which is impossible. So P2 stops at Pmax / 2 P2 = 0 while P2 <= Pmax / 2: P1 = 0 # We want to stop before the sum of BER is greater than Pmax # Since P3 >= P2, the condition "P2 <= (Pmax - P2 - P1)" must be verified # so that, at worst, P3 = P2 and we don't use more than the allocated power while P1 <= P2 and P2 <= (Pmax - P2 - P1): P = [P1, P2, Pmax - P2 - P1] ... ... @@ -91,37 +100,39 @@ def graph_3D(step=Pmax / 100, block=False): P2 += step # On transforme en array numpy # Convert the lists to numpy arrays X = np.array(X) Y = np.array(Y) Z = np.array(Z) # On affiche tous les points # Display all the datapoints ax = plt.axes(projection="3d") ax.scatter3D(X, Y, Z, c="blue") ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("Somme des BER") ax.set_title("BER total à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = [{:.3g}, {:.3g}, {:.3g}] et σ = [{:.3g}, {:.3g}, {:.3g}]".format(N, Pmax, g[0], g[1], g[2], sigma[0], sigma[1], sigma[2])) ax.set_zlabel("Sum of BER") ax.set_title("Evolution of the sum of BER for N={} users depending on the power, with Pmax = {}, g = [{:.3g}, {:.3g}, {:.3g}] and σ = [{:.3g}, {:.3g}, {:.3g}]".format(N, Pmax, g[0], g[1], g[2], sigma[0], sigma[1], sigma[2])) # On trouve le min et on l'affiche # Find the optimum opti = np.argmin(Z) P = [X[opti], Y[opti], Pmax - X[opti] - Y[opti]] BERmin = Z[opti] print("Graphically, a min was found for P = [{:.2g}, {:.2g}, {:.2g}] with BER={:.3g}".format(P[0], P[1], P[2], BERmin)) # Display the min we found ax.plot([P[0]], [P[1]], [BERmin], markerfacecolor="r", marker="o", markersize=10) ax.text(P[0], P[1], BERmin - 0.02, "P=({:.2g}, {:.2g}, {:.2g}), BER={:.2g}".format(P[0], P[1], P[2], BERmin), fontsize=8, horizontalalignment="center", verticalalignment="center") ax.text(P[0], P[1], BERmin - 0.05, "P=({:.2g}, {:.2g}, {:.2g}), BER={:.2g}".format(P[0], P[1], P[2], BERmin), fontsize=8, horizontalalignment="center", verticalalignment="center") # Afficher la fenêtre sans bloquer l'interface # Display the output plt.show(block) # Renvoyer la valeur trouvée # Return the value return P, BERmin def graph_4D(step=Pmax / 100, block=False): """ Graph le min pour N = 4 Graph the optimum for N = 4 """ X = [] Y = [] ... ... @@ -129,12 +140,19 @@ def graph_4D(step=Pmax / 100, block=False): C = [] ax = plt.axes(projection="3d") # On calcule f(P) pour les valeurs valides de P # Compute f(P) for all possible power repartitions # The constraints are Pi > 0 and Pi+1 < Pi < Pi-1 # If P3 > Pmax / 2, then P4 < P3, which is impossible. So P3 stops at Pmax / 2 P3 = P2 = P1 = 0 while P3 <= Pmax / 2: P2 = 0 # We want to stop before the sum of BER is greater than Pmax # Since P4 >= P3, the condition "P3 <= (Pmax P3 - P2 - P1)" must be verified # so that, at worst, P4 = P3 and we don't use more than the allocated power while P2 <= P3 and 2 * P3 + P2 + P1 <= Pmax: P1 = 0 # Same as before, but with P1 and P2 while P1 <= P2 and 2 * P3 + P2 + P1 <= Pmax: P = [P1, P2, P3, Pmax - P3 - P2 - P1] ... ... @@ -147,17 +165,18 @@ def graph_4D(step=Pmax / 100, block=False): P2 += step P3 += step # On transforme en array numpy # Convert the lists to numpy arrays X = np.array(X) Y = np.array(Y) Z = np.array(Z) C = np.array(C) # On trouve le min et le max # Find the min and max values for BER (aka the color) Cmax = max(C) + 0.1 Cmin = min(C) - 0.1 # On affiche toutes les valeurs # Display all the datapoints # Use colors for the 4th dimension for i in range(len(X)): x = X[i] y = Y[i] ... ... @@ -165,73 +184,101 @@ def graph_4D(step=Pmax / 100, block=False): color = [0, 0, (Cmax - C[i]) / (Cmax - Cmin)] ax.plot([x], [y], [z], c=color, marker=".", markersize=10, alpha=0.7) ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("Sum of BER") ax.set_title("Evolution of the sum of BER for N={} users depending on the power, with Pmax = {}, g = [{:.3g}, {:.3g}, {:.3g}, {:.3g}] et σ = [{:.3g}, {:.3g}, {:.3g}, {:.3g}]".format(N, Pmax, g[0], g[1], g[2], g[3], sigma[0], sigma[1], sigma[2], sigma[3])) print("Graphically, a min was found for P = [{:.2g}, {:.2g}, {:.2g}, {:.2g}] with BER={:.3g}".format(P[0], P[1], P[2], P[3], BERmin)) # Display a colorbar sm = cm.ScalarMappable(cmap=cm.Blues, norm=plt.Normalize(vmin=0, vmax=1)) sm._A = [] cbar = plt.colorbar(sm) cbar.set_label("P3", rotation=0) # On trouve le min et on l'affiche # Find the optimum opti = np.argmin(Z) P = [X[opti], Y[opti], C[opti], Pmax - X[opti] - Y[opti] - C[opti]] BERmin = Z[opti] ax.set_xlabel("P1") ax.set_ylabel("P2") ax.set_zlabel("Somme des BER") ax.set_title("BER total à N={} utilisateurs en fonction des puissances avec Pmax = {}, g = [{:.3g}, {:.3g}, {:.3g}, {:.3g}] et σ = [{:.3g}, {:.3g}, {:.3g}, {:.3g}]".format(N, Pmax, g[0], g[1], g[2], g[3], sigma[0], sigma[1], sigma[2], sigma[3])) print("Graphically, a min was found for P = [{:.2g}, {:.2g}, {:.2g}, {:.2g}] with BER={:.3g}".format(P[0], P[1], P[2], P[3], BERmin)) # Display the min we found ax.plot([P[0]], [P[1]], [BERmin], markerfacecolor="r", marker=".", markersize=15) ax.text(P[0], P[1], BERmin - 0.05, "P=({:.2g}, {:.2g}, {:.2g}, {:.2g}), BER={:.2g}".format(P[0], P[1], P[2], P[3], BERmin), fontsize=8, horizontalalignment="center", verticalalignment="center") # Afficher la fenêtre sans bloquer l'interface # Display the output plt.show(block) # Renvoyer la valeur trouvée # Return the value return P, BERmin def evaluate_min(step=Pmax / 100): """ Find the min for more than 4 users """ global nbrSteps, BERmin, Pmin # Init default values Pmin = [0] * (N - 1) BERmin = f(Pmin + [Pmax]) # Init global variable to count the number of steps this took nbrSteps = 1 def should_continue(P, i, step): """ Returns True if another value of power for i should be tested """ # Make sure we don't use more than the max allocated power if sum(P) + P[-1] + step > Pmax: return False # Make sure the power remains smaller than the one of the previous user if i == len(P) - 1: # This is the first-but-one user (the last user's power is not in P since their power can be deduced) # Just don't take-up more than half of the total power return P[i] <= Pmax / 2 - step else: # Make sure that, after adding step, the power is still smaller than the one of the previous user return P[i] <= P[i + 1] - step def rec(P, i, step): """ Recursive method to compute all the BER values """ global nbrSteps, BERmin, Pmin # Make sure the user is valid if i < 0: return