un llyod max qui marche

455-Codage_Sources/algo_code/llyod_max.py

@@ -1,58 +1,53 @@
 #!/usr/bin/env python3
 import numpy as np
-from sipy import integrate
-from scipy import norm
-
-M = 8
-X = np.random.normal(0,1,1000)
+from scipy import integrate
+from scipy.stats import norm
+import matplotlib.pyplot as plt
 
 def ddp(x):
     mean = 0,
     sigma = 1
     return norm.pdf(x,mean,sigma)
 
-def init_thres_vec(M,X):
-    step = (np.max(X)-np.min(X))/M
-    thres_intervals = np.array([])
-    mid = np.mean(X)
-    for i in range(int(M/2)):
-        thres_intervals = np.append(thres_vec,mid+(i+1)*step)
-        thres_intervals = np.insert(thtres_vec,0,mid-(1+1)*step)
-    return thres_intervals
-
-def quant(x,thres,intervals):
-    thres= np.append(thres, np.inf)
-    thres= np.insert(thres, 0, -np.inf)
-    x_hat_q = np.zeros(np.shape(x))
-    for i in range(len(thres)-1):
-        if i == 0:
-            x_hat_q = np.where(np.logical_and(x > thres[i], x <= thres[i+1]),
-                                  np.full(np.size(x_hat_q), intervals[i]), x_hat_q)
-        elif i == range(len(thres))[-1]-1:
-            x_hat_q = np.where(np.logical_and(x > thres[i], x <= thres[i+1]),
-                               np.full(np.size(x_hat_q), intervals[i]), x_hat_q)
-        else:
-            x_hat_q = np.where(np.logical_and(x > thres[i], x < thres[i+1]),
-                               np.full(np.size(x_hat_q), intervals[i]), x_hat_q)
-    return x_hat_q
-
-
-def LlyodMax(X,intervals, max_iter=1000,eps=1e-5):
-    err_min = np.inf
-    for i in range(max_iter):
-        for j in range(len(x_hat_q)):
-            centroids[i] = integrate.quad(lambda x : x*ddp(x),
-                                        intervals[j],intervals[j+1])[0]/
-                         integrate.quad(lambda x : ddp(x),
-                                        intervals[j],intervals[j+1])[0]
-            intervals = 0.5*(centroids[1:]+centroids[:-1])
-            x_hat = quant(X,centroids,intervals)
-            err = np.linalg.norm(X-x_hat)
-            if err < err_min:
-                err_min =err
-                intervals_min = intervals
-                centroids_min = centroids
-        if err_min< 1e-5:
+def quant(centroids, X):
+    bornes = (centroids[:-1]+centroids[1:])/2
+    bornes = np.insert(bornes,0,-np.inf)
+    bornes = np.append(bornes,np.inf)
+    xquant =np.zeros(len(X))
+    for k in range(len(X)):
+        for i in range(len(bornes)):
+            if X[k]>=bornes[i] and X[k] <bornes[i+1]:
+                xquant[k] = centroids[i]
+    return xquant
+def llyodMax(X,M,maxiter=1000,eps=1e-6):
+    #répartition uniforme des bornes
+    step = (np.max(X)-np.min(X))/(M-2)
+    Xmin = np.min(X)
+    bornes = np.array([i*step+Xmin for i in range(M-1)])
+    bornes = np.insert(bornes,0,-np.inf)
+    bornes = np.append(bornes,np.inf)
+    centroids = np.zeros(M)
+    for  it in range(maxiter):
+        old_centroids = centroids.copy()
+        for i in range(M): 
+            centroids[i] = integrate.quad(lambda x: x*ddp(x),bornes[i],bornes[i+1])[0]\
+                         /integrate.quad(lambda x: ddp(x),bornes[i],bornes[i+1])[0]
+        bornes[1:-1] = (centroids[:-1]+centroids[1:])/2
+        err = np.linalg.norm(centroids-old_centroids)
+        print(err)
+        if err < eps :
             break
-    best_x_hat = quant(X,centroids_min,intervals_min)
-    return best_x_hat
+    return centroids    
+
+M = 4
+X = np.random.normal(0,1,1000)
+centroids = llyodMax(X,M)
+bornes = (centroids[:-1]+centroids[1:])/2
+bornes = np.insert(bornes,0,-np.inf)
+bornes = np.append(bornes,np.inf)
+
+print(centroids, bornes)
+plt.figure()
+plt.plot(X)
+plt.plot(quant(bornes,X))
+plt.show()