Commit 15d8bc99 by JRock007

### Removed unused python files

parent 25072c50
 import Theorie_N_recepteurs as NOMATh import numpy as np import time from scipy.optimize import minimize from scipy.optimize import basinhopping from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import matplotlib.cm as cm """ Minimize f(P) Subject to: ∀ i ∈ ⟦1, N⟧, Pi ⩾ 0 ∀ (i, j) ∈ ⟦1, N⟧², i ≤ j ⇒ Pi ≤ Pj ΣPi ≤ Pmax It can be written as: min(f(P)) under the constraint that AX ⩾ B and ∀ i ∈ ⟦1, N⟧, Pi ⩾ 0 A = (-1, 1, 0, ..., 0, 0) ( 0, -1, 1, ..., 0, 0) ... ( 0, 0, 0, ..., -1, 1) (-1, -1, -1, ..., -1, -1) X = t(P1, P2, ..., PN) B = t(0, 0, ..., 0, -Pmax) 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] print("N = {}, Pmax = {}, g = {} and σ = {}".format(N, Pmax, g, sigma)) def graph_2D(n=100, block=False): """ Graph le min pour N = 2 """ # On génère les puissances possibles et on calcule f(P) X = np.linspace(0, Pmax / 2, n) Y = [f([x, Pmax - x]) for x in X] # On affiche les points 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])) # On trouve le min et on l'affiche opti = np.argmin(Y) P = [X[opti], Pmax - X[opti]] BERmin = Y[opti] print("Graphically, min was found for P1 = [{:.2g}, {:.2g}] with BER={:.3g}".format(P[0], P[1], BERmin)) 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") # Afficher la fenêtre sans bloquer l'interface plt.show(block) # Renvoyer la valeur trouvée return P, BERmin def graph_3D(step=Pmax / 100, block=False): """ Graph le min pour N = 3 """ X = [] Y = [] Z = [] # On calcule f(P) pour les valeurs valides de P P2 = 0 while P2 <= Pmax / 2: P1 = 0 while P1 <= P2 and P2 <= (Pmax - P2 - P1): P = [P1, P2, Pmax - P2 - P1] X.append(P1) Y.append(P2) Z.append(f(P)) P1 += step P2 += step # On transforme en array numpy X = np.array(X) Y = np.array(Y) Z = np.array(Z) # On affiche tous les points 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])) # On trouve le min et on l'affiche 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)) 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") # Afficher la fenêtre sans bloquer l'interface plt.show(block) # Renvoyer la valeur trouvée return P, BERmin def graph_4D(step=Pmax / 100, block=False): """ Graph le min pour N = 4 """ X = [] Y = [] Z = [] C = [] ax = plt.axes(projection="3d") # On calcule f(P) pour les valeurs valides de P P3 = P2 = P1 = 0 while P3 <= Pmax / 2: P2 = 0 while P2 <= P3 and 2 * P3 + P2 + P1 <= Pmax: P1 = 0 while P1 <= P2 and 2 * P3 + P2 + P1 <= Pmax: P = [P1, P2, P3, Pmax - P3 - P2 - P1] X.append(P1) Y.append(P2) C.append(P3) Z.append(f(P)) P1 += step P2 += step P3 += step # On transforme en array numpy X = np.array(X) Y = np.array(Y) Z = np.array(Z) C = np.array(C) # On trouve le min et le max Cmax = max(C) + 0.1 Cmin = min(C) - 0.1 # On affiche toutes les valeurs for i in range(len(X)): x = X[i] y = Y[i] z = Z[i] color = [0, 0, (Cmax - C[i]) / (Cmax - Cmin)] ax.plot([x], [y], [z], c=color, marker=".", markersize=10, alpha=0.7) 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 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)) 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 plt.show(block) # Renvoyer la valeur trouvée return P, BERmin def evaluate_min(step=Pmax / 100): global nbrSteps, BERmin, Pmin Pmin = [0] * (N - 1) BERmin = f(Pmin + [Pmax]) nbrSteps = 1 def should_continue(P, i, step): if sum(P) + P[-1] + step > Pmax: return False if i == len(P) - 1: return P[i] <= Pmax / 2 - step else: return P[i] <= P[i + 1] - step def rec(P, i, step): global nbrSteps, BERmin, Pmin if i < 0: return while should_continue(P, i, step): P[i] += step nbrSteps += 1 BER = f(P + [Pmax - sum(P)]) if BER < BERmin: BERmin = BER Pmin = P.copy() rec(P.copy(), i - 1, step) rec(Pmin.copy(), N - 2, step) print("Iterativally, a min was found for P = {} in {} iterations".format(Pmin + [Pmax - sum(Pmin)], nbrSteps)) return Pmin + [Pmax - sum(Pmin)], BERmin def global_min(x0, bounds, constraints=None, niter=100): """ Trouve le min global avec basinhopping """ # Stepsize must be big enough so it doesn't get stuck in a local min ret = basinhopping(f, x0=x0, niter=niter, stepsize=Pmax / 10, minimizer_kwargs={"bounds": bounds, "constraints": constraints}) print("The global min was found for P =", ret.x) print("BER =", f(ret.x)) def local_min(x0, bounds, constraints=None, method="SLSQP"): """ Trouve un min loca avec minimize """ if method == "COBYLA": # Cobylla can't handle bounds ret = minimize(f, x0=x0, method=method, constraints=constraints) else: ret = minimize(f, x0=x0, bounds=bounds, method=method, constraints=constraints) print("A local min was found for P =", ret.x) print("BER =", f(ret.x)) def f(P): """ Fonction à minimiser Entrée : Liste les puissances Sortie : float """ # Chaque utilisateur ayant une valeur de g différente, on recalcule la proba à chaque fois probas = [] for i in range(N): proba = NOMATh.theorie(g[i], sigma[i], P, N, nuser=i + 1) # Une probabilité "None" indique un problème de paramètres, on renvoie +inf if proba is None: return np.inf else: probas.append(proba) return sum(probas) def constraint(P): """ Fonction permettant de tester si les contraintes sont vérifiées Entrée : liste des puissances Sortie : Matrice dont toutes les valeurs doivent être positives ou nulles """ R = np.matmul(A, P.T) return R - B # Création de la matrice A A = np.zeros((N, N)) for i in range(N - 1): A[i, i] = -1 A[i, i + 1] = 1 A[-1, :] = -np.ones((1, N)) # Création de la matrice B B = np.zeros(N) B[-1] = -Pmax # Dictionaire permettant d'indiquer que c'est la fonction "constraint" qui calcule les contraintes cons = {"type": "ineq", "fun": constraint} # ineq: the constraint function result is to be non-negative # Création des brones sur chaque puissance (les puissances doivent être entre 0 et +∞) # Note : pas forcément utile car f(P) renvoie +∞ si une puissance est négative bounds = [(0, np.inf)] * N # Première méthode : Évaluation de la fonction en plusieur points et résolution graphique puis recherche d'un minimum local print("----- Finding graphical minimum -----") t0 = time.time() if N == 2: Pmin, BERmin = graph_2D() elif N == 3: Pmin, BERmin = graph_3D(step=Pmax / 100) elif N == 4: Pmin, BERmin = graph_4D(step=Pmax / 50) else: Pmin, BERmin = evaluate_min(step=Pmax / 25) print("Finding local min starting with this point...") local_min(Pmin, bounds, constraints=cons) print("Took {:.2g}s".format(time.time() - t0)) # Deuxième méthode : recherche d'un minimum global print("----- Finding global minimum -----") x0 = np.zeros(N) niter = 20 print("P0 =", x0) t0 = time.time() global_min(x0, bounds, constraints=cons, niter=niter) print("Took {:.2g}s for {} iterations".format(time.time() - t0, niter)) # Troisème méthode : recherche d'un minimum local à partir d'un point arbitraire print("----- Finding local minimum with arbitrary starting point -----") # Le choix des conditions initiales est important x0 = np.array([2**i for i in range(N)], dtype="float64") x0 *= Pmax / sum(x0) print("P0 =", x0) t0 = time.time() local_min(x0, bounds, constraints=cons) print("Took {:.2g}s".format(time.time() - t0)) # Afficher la fenêtre pyplot et empêcher Python de fermer si c'est utile if N < 5: plt.show(True)
 import math import matplotlib.pyplot as plt import numpy as np import Theorie_N_recepteurs as NOMATh def comp_log_log(xlog=True, ylog=False): """ Graph """ xlog = True ylog = False # Init variables N = 2 Pmax = 1 g = 1 Sigmas = np.linspace(math.sqrt(Pmax / 50), math.sqrt(Pmax), 1000) RSB = Pmax / np.array(np.square(Pmax / Sigmas)) # We look for a minimum / RSB function Opt_power = [] # P1 opt w/ respect to sigma Opt_ber = [] # ber for this P1 opt w/ respect to sigma for sigma in Sigmas: # Probas1 = [] # Probas2 = [] Sum = [] # Total BER for given sigma, in 0 -> 50% Pmax range for i in range(50): P1 = i * Pmax / 100 P2 = Pmax - P1 P = [P1, P2] probas = NOMATh.theorie(g, sigma, P, N) # Probas1.append(probas[0]) # Probas2.append(probas[1]) Sum.append(probas[0] + probas[1]) minima = Sum.index(min(Sum)) Opt_ber.append(min(Sum) * 100) Opt_power.append(minima) fig = plt.figure() ax = fig.add_subplot(1, 1, 1) if xlog: ax.set_xscale("log") plt.xlabel("RSB (log)") else: plt.xlabel("P1 / Pmax") if ylog: ax.set_yscale("log") plt.ylabel("BER (log)") else: plt.ylabel("BER") ax.plot(RSB, Opt_power, label="Optimum") ax.plot(RSB, Opt_ber, label="Ber") plt.legend() plt.show() if __name__ == '__main__': comp_log_log() \ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment