# -*- coding: utf-8 -*-
"""
Created on Fri Oct 20 11:54:42 2017

@author: florent
http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html
"""


import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as spl
import sklearn.mixture  

#Ces deux lignes permettent de regler le nombre de decimales dans les "print" de tableaux np:
float_formatter = lambda x: "%.4f" % x
np.set_printoptions(formatter={'float_kind':float_formatter})

#Parametres fixes par l'utilisateur:
rangeK=range(1,10) # valeurs de K testees pour le nombre de classes (dans la selection de modele)

#Creation du tableau de donnees "data" a partir du fichier:
#f=open("densitesOs.txt","r")
#datatxt=f.readlines()
#data=[]
#for line in datatxt:
#    line_of_data=[]
#    line = line.strip() #Returns a copy of the string line with possible leading and trailing whitespace removed
#    for number in line.split(): #(line.split() est la liste des str obtenus en scidant line aux espaces)
#        line_of_data.append(float(number))
#    data.append(np.array(line_of_data))
#Beaucoup, beaucoup plus simple :
data=np.loadtxt("densitesOs.txt")
n=len(data)
d=len(data[0])


#Fonction utile pour materialiser la covariance et la moyenne d'un echantillon gaussien:
def coord_elliposide(mu=[0,0],C=np.eye(2),u=.5,m=1000):
    r=-2*np.log(u)
    A=spl.sqrtm(C)
    a,b,c=A[0,0],A[0,1],A[1,1]
    t=np.linspace(0,2*np.pi,m)
    co,si=r*np.cos(t),r*np.sin(t)
    x=mu[0]+a*co+b*si
    y=mu[1]+b*co+c*si
    return (x,y)

#Affichage des donnees:
plt.figure()
data=np.array(data)
plt.plot(data[:,0],data[:,1],"r.")
plt.title("Les "+str(n)+" points du data set" )
plt.show()


lowest_bic = np.infty
bic=[]
for k in rangeK:
    gmm = sklearn.mixture.GaussianMixture(n_components=k, covariance_type='full')
    gmm.fit(data)
    bic.append(gmm.bic(data))
    if bic[-1] < lowest_bic:
        lowest_bic = bic[-1]
        best_gmm = gmm
        
l=np.argmin(bic)
Koptimal=rangeK[l]
alpha_em,mu_em,covariance_em=best_gmm.weights_, best_gmm.means_, best_gmm.covariances_ 
Classes = best_gmm.predict(data)

#Affichage des donnees colorees en fonction de leur classe telle qu'estimee par la fct EM:
plt.figure()
for k in xrange(Koptimal):
    datak=np.array([data[i] for i in np.where(Classes == k)[0]])
    plt.plot(datak[:,0],datak[:,1],".")
    (x,y)=coord_elliposide(mu_em[k],covariance_em[k])
    plt.plot(x,y,"r")
    plt.plot(mu_em[k][0],mu_em[k][1],"r+")        
plt.title("Les "+str(n)+" points du data set clusterises par EM" )
plt.show()

#Affichage de la fonction dont le "K" optimal est l'argmin
plt.figure()
plt.plot(rangeK,bic, "r")
plt.title("Model selection")
plt.show()


#Affichage des parametres du modele tels qu'estimes par la fct EM:
print("\nValeur de K :\n")
print(Koptimal)

print("\nEstimation de alpha :\n")
print("%.2f, "*Koptimal % tuple(alpha_em))

print("\nEstimation de mu ")
for k in xrange(Koptimal):
    print("\nClasse "+str(k)+" :")
    print("%.2f, "*d % tuple(mu_em[k]))

print("\nEstimation de la covariance :")
for k in xrange(Koptimal):
    print("\nClasse "+str(k)+" :")
    print(covariance_em[k])