diff --git a/455-Codage_Sources/algo_code/llyod_max.py b/455-Codage_Sources/algo_code/llyod_max.py index f34eb3d..0c9d1a9 100755 --- a/455-Codage_Sources/algo_code/llyod_max.py +++ b/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]