Added back Optimization Python script

Python/Recherche_Optimum.py

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