import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from sklearn import mixture
import math


n_samples = 300
dimension_de_x =2
k=2

def Expectation (data, sigma, mu, phi, i, j):

    #numerateur
    determinant = np.linalg.det(sigma[j])
    determinant = np.absolute(determinant)
    dataligne = [data[i,]]
    dataligne = np.matrix(np.array(dataligne))
    dataligne = dataligne - np.matrix(np.array([mu[j,]]))
    transpose = dataligne.transpose()

    num = (1/((2*math.pi)*math.sqrt(determinant)))
    num = num*math.exp(-((dataligne)*np.linalg.inv(sigma[j])*transpose))*phi[0,j]
    #print("le numerateur est",num)


    #denominateur
    determinant = np.linalg.det(sigma[0])
    determinant = np.absolute(determinant)
    den1 = (1/((2*math.pi)*math.sqrt(determinant)))
    dataligne = [data[i,]]
    dataligne = np.matrix(np.array(dataligne))
    dataligne = dataligne - np.matrix(np.array([mu[0,]]))
    transpose = dataligne.transpose()
    den1 = den1*math.exp(-((dataligne)*np.linalg.inv(sigma[0])*transpose))*phi[0,0]

    determinant = np.linalg.det(sigma[1])
    determinant = np.absolute(determinant)
    den2 = (1/((2*math.pi)*math.sqrt(determinant)))
    dataligne = [data[i,]]
    dataligne = np.matrix(np.array(dataligne))
    dataligne = dataligne - np.matrix(np.array([mu[1,]]))
    transpose = dataligne.transpose()
    den2 = den2*math.exp(-((dataligne)*np.linalg.inv(sigma[1])*transpose))*phi[0,1]

    den=den1+den2

    #print(den)

    return num/den




def Maximisation_phi (w, phi):

    sommekluster1 = 0
    sommekluster2 = 0

    for i in range (600):

        #print("la valeur en i ",i,"et j 1 est ",w[i,0])
        sommekluster1 = sommekluster1 + w[i,0]

        #print("la valeur en i ",i,"et j 2 est ",w[i,1])
        sommekluster2 = sommekluster2 + w[i,1]

    #print("Cluster 1 ",sommekluster1/600," Cluster 2 ",sommekluster2/600)

    phi[0,0]=sommekluster1/600
    phi[0,1]=sommekluster2/600

    return phi



def Maximisation_mu (w,data,mu):

    sommekluster1 = 0
    sommekluster2 = 0

    sommeproba1 = 0
    sommeproba2 = 0

    for i in range (600):

        #print(dataligne)
        sommekluster1=sommekluster1+w[i,0]*data[i]
        sommekluster2=sommekluster2+w[i,1]*data[i]

        sommeproba1=sommeproba1+w[i,0]
        sommeproba2=sommeproba2+w[i,1]

    mu[0]=sommekluster1/sommeproba1
    mu[1]=sommekluster2/sommeproba2

    print("Somme des valeur en cluster 1", sommekluster1)
    print("Somme des valeur en cluster 2", sommekluster2)
    print("Somme en valeur cluster 1", sommeproba1)
    print("Somme en valeur cluster 2", sommeproba2)

    return mu



def Maximisation_sigma (w,data,mu,sigma):

    sommekluster1 = 0
    sommekluster2 = 0
    temp1 = 0
    temp2 = 0

    sommeproba1 = 0
    sommeproba2 = 0

    for i in range (600):

        dataligne = [data[i,]]
        dataligne = np.matrix(np.array(dataligne))
        dataligne = dataligne - np.matrix(np.array([mu[j,]]))

        transpose = dataligne.transpose()

        #print("prob ", w[i,0])
        #print("ecart au centroid", dataligne)
        #print("transpose", transpose)

        temp1=w[i,0]*transpose*dataligne

        sommekluster1=sommekluster1+temp1
        sommeproba1=sommeproba1+w[i,0]

        sigma[0]=sommekluster1/sommeproba1



        dataligne = [data[i,]]
        dataligne = np.matrix(np.array(dataligne))
        dataligne = dataligne - np.matrix(np.array([mu[j,]]))

        transpose = dataligne.transpose()

        #print("prob ", w[i,0])
        #print("ecart au centroid", dataligne)
        #print("transpose", transpose)

        temp2=w[i,1]*transpose*dataligne

        sommekluster2=sommekluster2+temp2
        sommeproba2=sommeproba2+w[i,1]

        sigma[1]=sommekluster2/sommeproba2


    return sigma



# generate random sample, two components
np.random.seed(1)

# generate spherical data centered on (20, 20)
shifted_gaussian = np.random.randn(n_samples, dimension_de_x) + np.array([20, 20])

# generate zero centered stretched Gaussian data
C = np.array([[0., -0.7], [3.5, .7]])
stretched_gaussian = np.dot(np.random.randn(n_samples, dimension_de_x), C)

#probablities
w = np.random.randn(600, 2)

# concatenate the two datasets into the final training set
X_train = np.vstack([shifted_gaussian, stretched_gaussian])

sigma = np.array([np.random.rand(k,dimension_de_x), np.random.rand(k,dimension_de_x)])
mu = np.array([np.random.rand(1,dimension_de_x), np.random.rand(1,dimension_de_x)])
phi = np.matrix([[ 0.3333333, 0.6666667]])

for iteration in range(10):
    for i in range(600) :
        for j in range(2):

            test=Expectation(X_train,sigma,mu,phi,i,j)
            w[i,j]=test


    print ("Before maximisation phi ", phi)
    phi=Maximisation_phi(w,phi)
    print("After maximisation phi ", phi)

    print ("Before maximisation sigma ", sigma)
    sigma=Maximisation_sigma(w,X_train,mu,sigma)
    print ("After maximisation sigma ", sigma)

    print ("Before maximisation mu ", mu)
    mu=Maximisation_mu(w,X_train,mu)
    print("After maximisation mu ", mu)