#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import itertools
import sys

def basicModel():
	'''Returns toy model for debugging

	Returns:
		(
			np.ndarray: 1D array with clustering coefficients
			np.ndarray: 3D array with Dirichlet distribution for each answer regarding
			    a (cluster, question)
		)

	Returns
	-------
	(mixtureCoefs, probs) où
	- mixtureCoefs is a 1D vector holding the mixture coefficients
	- probs is a 3D array (nClusters, nQuestions, nAnswers)
		holding Dirichlet distribution associated to answers
		for each point (cluster, question)
	'''
	
	# Mixture Cluster Coefficients
	mixtureCoefs = np.array([0.4,0.6])
	
	# Probabilities of 3 answers to four possible questions given their clusters
	probs = np.array([
	# P(answer = j | question = i, cluster = 1)
	[
		[0.2 , 0.5, 0.3],
		[0.1 , 0.2, 0.7 ],
		[0.7 , 0.1, 0.2 ],
		[0.3 , 0.3, 0.4 ]
	],
	# P(answer = j | question = i, cluster = 2)
	[
		[0.7, 0.1, 0.2 ],
		[0.2, 0.6, 0.2 ],
		[0.2, 0.4, 0.4 ],
		[0.2, 0.6, 0.2 ]
	]])
	
	# Model
	return (mixtureCoefs, probs)

def generateModel(nQuestions, nAnswers, nClusters):
	''' Generate a mixture model with Dirichlet distribution
	
	Args:
		nQuestion(int): Number of questions
		nAnswers(int): Number of answers
		nClusters(int): Number of clusters
	
	Returns:
	    (
			np.ndarray: Mixture coefficient array
			np.ndarray: 3D array with dirichlet distribution of answer given their question and cluster
		)
	'''
		
	mixtureCoefs = np.random.dirichlet(np.ones(nClusters))
	probs = np.random.dirichlet(np.ones(nAnswers), (nClusters, nQuestions))
	return (mixtureCoefs, probs)
	
def generateData(model, nPeople):
	''' Generate answer to a survey provided with a model

	Args:
		model(tuple): model generated by 'generateModel'
		nPeople: Number of people who answered the survey
		
	Returns:
		(
			np.ndarray: nPeople*nQuestions array with answers to each question
			np.ndarray: nPeople array with the cluster associated to each person
		)
	'''
	
	mixtureCoefs, probs = model
	nClusters, nQuestions, nAnswers = probs.shape
	clusters = np.random.choice(nClusters, size = nPeople, p = mixtureCoefs)
	data = np.zeros((nPeople, nQuestions))
	answers = range(1, nAnswers+1)
	for i in range(nPeople):
		for j in range(nQuestions):
			data[i,j] = np.random.choice(answers, p=probs[clusters[i],j,:])	
	return (np.asarray(data, 'uint'), clusters + 1)

def drawModel(model):
	''' Plot distribution for each clusters

	Args:
		model(tuple): Model output by 'generateModel'
	'''
	
	plt.ion()
	mixtureCoefs, probs = model
	nClusters, nQuestions, nAnswers = probs.shape	
	questions = np.arange(nQuestions)
	fig = plt.figure()
	for cluster in range(nClusters):
		plt.subplot(nClusters, 1, cluster+1)	
		colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k'])
		total = np.zeros(nQuestions)		
		for answer in range(nAnswers):
			c = next(colors)
			d = probs[cluster][:,answer]
			plt.bar(questions, d, width=1, bottom = total, color=c)
			total += d
		plt.yticks(np.arange(0, 1.0001, 0.5))
		plt.ylim((0,1))
		plt.title('Groupe ' + str(cluster+1))

	plt.xlabel('Questions', fontsize=14)
	fig.tight_layout()		
	plt.show()

def drawCluster(H):
	''' Plot distribution of latent variables
	H can be either a n*1 or n*nc matrix.

	The number n of rows of H is the number of persons
	* If H is a n*nc matrix, 
		the number of columns is the number of clusters
		and the components H[i,j] gives the likelihood 
		P(Person i belongs to cluster j)
	* If H is n*1 matrix, 
		the component H[i] gives the most probable cluster to 
		which person i belongs
	
	Args:
	    H(np.ndarray): (soft or hard) Clusters assignements
	'''
	
	plt.ion()
	if(len(H.shape) >= 2):
		nPeople, nClusters = H.shape
	else:
		nPeople, = H.shape
		nClusters = np.max(H)
		Hnew = np.zeros((nPeople, nClusters))
		for (i,j) in zip(range(nPeople), H-1):
			Hnew[i,j] = 1.
		H = Hnew
	fig = plt.figure()
	people = np.arange(nPeople)
	colors = itertools.cycle(['b', 'r', 'g', 'c', 'm', 'y', 'k'])
	total = np.zeros(nPeople)		
	for cluster in range(nClusters):
		c = next(colors)
		d = H[:,cluster]
		plt.bar(people, d, width=1, bottom = total, color=c)
		total += d
	plt.yticks(np.arange(0, 1.0001, 0.25))
	plt.ylim((0,1))
	plt.title('Distribution of latent variables')
	plt.xlabel('person', fontsize=14)
	fig.tight_layout()		
	plt.show()

def compareClusters(C1,C2):
	''' Renvoie le degré de superposition de deux clusters C1 et C2 
		à une permutation près des clusters
		
		C1 et C2 donnent les clusters associés à chaque exemple.
		La ligne i de C1 et C2 correspond à l'exemple i.
		C1 et C2 peuvent être soit un vecteur soit une matrice.
		* Si C est un vecteur, C[i] contient le cluster auquel i appartient.
		* Si C est une matrice, le nombre de colonnes de C correspond au nombre de clusters
		et C[i,j] contient la probabilité que l'exemple i appartient au cluster j
		Dans ce cas le cluster le plus probable est associé à i
		
		La fonction renvoie (bestOverlap, bestPerm) où
		* bestOverlap est le degré de superposition de C1 et C2
		* bestPerm donne la correspondance des clusters de C2 avec ceux de C1
		bestPerm[n° du cluster de C2] = n° du cluster de C1 correspondant '''

	if(C1.shape[0] != C2.shape[0]):
		raise "C1 and C2 must have the same number of lines"
			
	if(len(C1.shape) >= 2 and C1.shape[1] > 1):
		C1 = np.argmax(C1,1)+1
	if(len(C2.shape) >= 2 and C2.shape[1] > 1):
		C2 = np.argmax(C2,1)+1
			
	nClusters = max(np.max(C1), np.max(C2))
	perms = itertools.permutations(range(1,nClusters+1))
	bestOverlap = 0.
	bestPerm = None
	for perm in perms:
		p = np.array(perm)
		permC2 = p[C2-1]
		overlap = np.sum(C1 == permC2) * 1. / np.size(C1)
		if(overlap > bestOverlap):
			bestOverlap = overlap
			bestPerm = p
	return (bestOverlap, bestPerm)
	
def permuteClusters(M, H, perm):
	'''Permute clusters on a model M
	and latent cluster distributions H using
	cluster mapping perm

	Args:
		M: Model
		H: Latent clusters
		perm: Cluster mapping 
	
	Returns:
		tuple: (permuted model, permuted H)
	'''
	invperm = np.argsort(perm-1)
	mixtureCoefs, probs = M
	mixtureCoefs = mixtureCoefs[invperm]
	probs = probs[invperm, :, :]
	H = H[:, invperm]
	return ((mixtureCoefs, probs), H)

###############
# A compléter #
###############

def em(data, nClusters, nIterations):
	''' EM Algorithm
	the data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		nIterations(int): Number of iteration to run the EM updates
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	'''

	nPeople, nQuestions = data.shape
	nAnswers = int(np.max(data))

	# q[i,j] = P(people i belong to Cluster j)
	# P(Hi|Vi = vi,Θ = θ) ~ qi
	# We store those distributions in an array
	q = np.zeros((nPeople, nClusters))

    # Mixture coefficients ph
	# Parameters of the distribution of the hidden variables
	# P(H) = ph(i) = P(Cluster = i) 
	ph = np.ones(nClusters)/nClusters

	# Probability of the choice of a visible variable
	# knowing that a given latent variable was chosen
	# P(V|H,θ) = pvh(i,j,k) = P(Answer = A(i,j) | Cluster = k)
	# Answer i,j is the the answer given by person i to question j
	# We could roughly model A by couple of variable (P,Q) i.e P answered Q
	pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range(nClusters)])
	# Parameters of different clusters distribution are all equal
	# pvh = np.array([np.ones((nQuestions, nAnswers)) for i in range(nClusters)])
	pvh = normalize(pvh)
    
	for t in range(nIterations):
	    # Step E 
	    # First We focus only on the distribution q to maximize the Likelihood 
		for i in range(nPeople):
			answers = data[i,:]
			for c in range(nClusters):
				# We use log since the probabilities might have small values
				# We use Bayes rule i.e H|V,θ = V|(H,θ) * H|θ / (V,θ) 
				# First we add the log of our prior H
				q[i,c] = np.log(ph[c])
                
				# Then we add the logs of Vi|H
				for j in range(nQuestions):
					q[i,c] += np.log(pvh[c,j,int(answers[j]-1)])
				
			q[i,:] -= np.max(q[i,:])
			q[i,:] = np.exp(q[i,:])
			# normalize
			q[i,:] /= np.sum(q[i,:])
        
		# Step M
		# The distribution of our cluster is the average of our the likelihoods distribution qi
		weights = np.sum(q,0)
		ph = weights/np.sum(weights)

		for cluster in range(nClusters):
			for j in range(nQuestions):
				for k in range(1, nAnswers + 1):
					s = 0
					#total = 0.
					for i in range(nPeople):
						if (data[i,j] == k):
							s+= q[i,cluster]
						#total += q[i,cluster]
					#pvh[cluster, j, k-1] = s/total
					pvh[cluster, j, k-1] = s/weights[cluster]
                
				# Normalization 
				pvh[cluster, j, :] /= sum(pvh[cluster, j, :])

		model = (ph, pvh)
		return (model, q)


def normalize(v):
	"""Normalize a 3D array

	Args:
		v(np.ndarray): 3D array
	
	Returns:
        np.ndarray: Normalized matrix
	"""
	dim = v.shape 
	for i in range(0, dim[0]-1):
		v[i,:,:] = (v[i,:,:].T/np.sum(v[i,:,:],1)).T

	return v

def loglikelihood(model, data, q):
	"""Computes and return the loglikelihood of a model M 
	given the data D and hidden variables distributions H (as returned by EM)

	Args:
		model: Mixture coefficient and visible variable distributions
		data: Matrix holding answers to each questions
		q: Latent cluster assignements
	
	Returns:
		float: Loglikelihood 
	"""
	ph, pvh = model
	nPeople, nQuestions = data.shape
	logL = 0
	for i in range(nPeople):
		answers = data[i,:]
		for k in range(nQuestions):
			logL += np.log(sum(pvh[:, k, int(answers[k] - 1)] * q[i,:].T))
	return logL

def em2(data, nClusters, epsilon):
	''' EM Algorithm with a stopping criterion. 
	The data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		epsilon(float): Threshold to detecte if the logLikelihood is still increasing
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	'''

	nPeople, nQuestions = data.shape
	nAnswers = int(np.max(data))

	# Each latent variable Hi follows a probability of distribution qi
	# q[i,j] = P(personne i belongs to Cluster j)
	# We store those distribution in an array
	q = np.zeros((nPeople, nClusters))

    # Mixture coefficients ph
	# Parameters of the distribution of the hidden variables
	# pc = ph(i) = P(Cluster = i) 
	ph = np.ones(nClusters)/nClusters

	# Probability of the choice of a visible variable
	# knowing that a given latent variable was chosen
	# pvh(i,j,k) = P(Reponse = A(i,j) | Cluster = k)
	pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range(nClusters)])
	# Parameters of different clusters distribution are all equal
	#pvh = np.array([np.ones((nQuestions, nAnswers)) for i in range(nClusters)])
	pvh = normalize(pvh)

	# Threshold epsilon
	oldLogL = 0;
	deltalogL = 100;
    
	while deltalogL > epsilon:

		# Step E 
		# First We focus only on the distribution q to maximize the MLE
		for i in range(nPeople):
			answers = data[i,:]
			for c in range(nClusters):
				# We use log since the probabilities might have small values 
				# First we add the log of our prior 
				q[i,c] = np.log(ph[c])
				
				# Then we add the log-likelihood of H|V
				for j in range(nQuestions):
					q[i,c] += np.log(pvh[c,j,int(answers[j]-1)])
				
			q[i,:] -= np.max(q[i,:])
			q[i,:] = np.exp(q[i,:])
			# normalize
			q[i,:] /= np.sum(q[i,:])
		
		# Step M
		weights = np.sum(q,0)
		ph = weights/np.sum(weights)

		for cluster in range(nClusters):
			for j in range(nQuestions):
				for k in range(1, nAnswers + 1):
					s = 0
					#total = 0.
					for i in range(nPeople):
						if (data[i,j] == k):
							s+= q[i,cluster]
						#total += q[i,cluster]
					#pvh[cluster, j, k-1] = s/total
					pvh[cluster, j, k-1] = s/weights[cluster]
				
				# Normalization 
				pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
		
		model = (ph, pvh)
		optLogL = loglikelihood(model, data, q)
		print(optLogL)
		deltalogL = abs(oldLogL - optLogL)
		oldLogL = optLogL
		return (model, q, optLogL)

	def BIC(model, logL, nPeople):
		"""
		Compute the Bayesian Information Criterion
		In order to  find the ideal number of clusters that produces accurate but not over-specific models, 
		we use a criterion that penalizes models with a large number of parameters. 
		The Bayesian Information Criterion (BIC) is one of these criteria. 
		While BIC has some theorical justifications for the exponential family of parameter distributions 
		whose categorical distributions are member of. We use it here as an empirical criterion.
		"""
		ph, pvh = model
		nClusters, nQuestions, nAnswers = pvh.shape

		nParams = (nClusters - 1) + nClusters * nQuestions *(nAnswers - 1)
		bic = -2. * logL + nParams * np.log(nPeople)
		return bic
	
	def findBestModel(D):
		nPeople, nQuestions = D.shape
		bestBic = 1E-1
		bestNClusters = 1
		for nClusters in range(1,10):
			M, H, logL = em2(D, nClusters, 1E-1)
			bic = BIC(M, logL, nPeople)
			if(bic < bestBic):
				bestBic = bic
				bestModel = M
				bestNClusters = nClusters
		
		return bestModel, bestNClusters





bff4a5f2b5b20b61ac6509f4d70603a6b0c91abf
<mask token>


def basicModel():
    """Returns toy model for debugging

	Returns:
		(
			np.ndarray: 1D array with clustering coefficients
			np.ndarray: 3D array with Dirichlet distribution for each answer regarding
			    a (cluster, question)
		)

	Returns
	-------
	(mixtureCoefs, probs) où
	- mixtureCoefs is a 1D vector holding the mixture coefficients
	- probs is a 3D array (nClusters, nQuestions, nAnswers)
		holding Dirichlet distribution associated to answers
		for each point (cluster, question)
	"""
    mixtureCoefs = np.array([0.4, 0.6])
    probs = np.array([[[0.2, 0.5, 0.3], [0.1, 0.2, 0.7], [0.7, 0.1, 0.2], [
        0.3, 0.3, 0.4]], [[0.7, 0.1, 0.2], [0.2, 0.6, 0.2], [0.2, 0.4, 0.4],
        [0.2, 0.6, 0.2]]])
    return mixtureCoefs, probs


def generateModel(nQuestions, nAnswers, nClusters):
    """ Generate a mixture model with Dirichlet distribution
	
	Args:
		nQuestion(int): Number of questions
		nAnswers(int): Number of answers
		nClusters(int): Number of clusters
	
	Returns:
	    (
			np.ndarray: Mixture coefficient array
			np.ndarray: 3D array with dirichlet distribution of answer given their question and cluster
		)
	"""
    mixtureCoefs = np.random.dirichlet(np.ones(nClusters))
    probs = np.random.dirichlet(np.ones(nAnswers), (nClusters, nQuestions))
    return mixtureCoefs, probs


def generateData(model, nPeople):
    """ Generate answer to a survey provided with a model

	Args:
		model(tuple): model generated by 'generateModel'
		nPeople: Number of people who answered the survey
		
	Returns:
		(
			np.ndarray: nPeople*nQuestions array with answers to each question
			np.ndarray: nPeople array with the cluster associated to each person
		)
	"""
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    clusters = np.random.choice(nClusters, size=nPeople, p=mixtureCoefs)
    data = np.zeros((nPeople, nQuestions))
    answers = range(1, nAnswers + 1)
    for i in range(nPeople):
        for j in range(nQuestions):
            data[i, j] = np.random.choice(answers, p=probs[clusters[i], j, :])
    return np.asarray(data, 'uint'), clusters + 1


def drawModel(model):
    """ Plot distribution for each clusters

	Args:
		model(tuple): Model output by 'generateModel'
	"""
    plt.ion()
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    questions = np.arange(nQuestions)
    fig = plt.figure()
    for cluster in range(nClusters):
        plt.subplot(nClusters, 1, cluster + 1)
        colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k'])
        total = np.zeros(nQuestions)
        for answer in range(nAnswers):
            c = next(colors)
            d = probs[cluster][:, answer]
            plt.bar(questions, d, width=1, bottom=total, color=c)
            total += d
        plt.yticks(np.arange(0, 1.0001, 0.5))
        plt.ylim((0, 1))
        plt.title('Groupe ' + str(cluster + 1))
    plt.xlabel('Questions', fontsize=14)
    fig.tight_layout()
    plt.show()


<mask token>


def em(data, nClusters, nIterations):
    """ EM Algorithm
	the data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		nIterations(int): Number of iteration to run the EM updates
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    for t in range(nIterations):
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        return model, q


def normalize(v):
    """Normalize a 3D array

	Args:
		v(np.ndarray): 3D array
	
	Returns:
        np.ndarray: Normalized matrix
	"""
    dim = v.shape
    for i in range(0, dim[0] - 1):
        v[i, :, :] = (v[i, :, :].T / np.sum(v[i, :, :], 1)).T
    return v


def loglikelihood(model, data, q):
    """Computes and return the loglikelihood of a model M 
	given the data D and hidden variables distributions H (as returned by EM)

	Args:
		model: Mixture coefficient and visible variable distributions
		data: Matrix holding answers to each questions
		q: Latent cluster assignements
	
	Returns:
		float: Loglikelihood 
	"""
    ph, pvh = model
    nPeople, nQuestions = data.shape
    logL = 0
    for i in range(nPeople):
        answers = data[i, :]
        for k in range(nQuestions):
            logL += np.log(sum(pvh[:, k, int(answers[k] - 1)] * q[i, :].T))
    return logL


def em2(data, nClusters, epsilon):
    """ EM Algorithm with a stopping criterion. 
	The data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		epsilon(float): Threshold to detecte if the logLikelihood is still increasing
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    oldLogL = 0
    deltalogL = 100
    while deltalogL > epsilon:
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        optLogL = loglikelihood(model, data, q)
        print(optLogL)
        deltalogL = abs(oldLogL - optLogL)
        oldLogL = optLogL
        return model, q, optLogL

    def BIC(model, logL, nPeople):
        """
		Compute the Bayesian Information Criterion
		In order to  find the ideal number of clusters that produces accurate but not over-specific models, 
		we use a criterion that penalizes models with a large number of parameters. 
		The Bayesian Information Criterion (BIC) is one of these criteria. 
		While BIC has some theorical justifications for the exponential family of parameter distributions 
		whose categorical distributions are member of. We use it here as an empirical criterion.
		"""
        ph, pvh = model
        nClusters, nQuestions, nAnswers = pvh.shape
        nParams = nClusters - 1 + nClusters * nQuestions * (nAnswers - 1)
        bic = -2.0 * logL + nParams * np.log(nPeople)
        return bic

    def findBestModel(D):
        nPeople, nQuestions = D.shape
        bestBic = 0.1
        bestNClusters = 1
        for nClusters in range(1, 10):
            M, H, logL = em2(D, nClusters, 0.1)
            bic = BIC(M, logL, nPeople)
            if bic < bestBic:
                bestBic = bic
                bestModel = M
                bestNClusters = nClusters
        return bestModel, bestNClusters

<mask token>


def basicModel():
    """Returns toy model for debugging

	Returns:
		(
			np.ndarray: 1D array with clustering coefficients
			np.ndarray: 3D array with Dirichlet distribution for each answer regarding
			    a (cluster, question)
		)

	Returns
	-------
	(mixtureCoefs, probs) où
	- mixtureCoefs is a 1D vector holding the mixture coefficients
	- probs is a 3D array (nClusters, nQuestions, nAnswers)
		holding Dirichlet distribution associated to answers
		for each point (cluster, question)
	"""
    mixtureCoefs = np.array([0.4, 0.6])
    probs = np.array([[[0.2, 0.5, 0.3], [0.1, 0.2, 0.7], [0.7, 0.1, 0.2], [
        0.3, 0.3, 0.4]], [[0.7, 0.1, 0.2], [0.2, 0.6, 0.2], [0.2, 0.4, 0.4],
        [0.2, 0.6, 0.2]]])
    return mixtureCoefs, probs


def generateModel(nQuestions, nAnswers, nClusters):
    """ Generate a mixture model with Dirichlet distribution
	
	Args:
		nQuestion(int): Number of questions
		nAnswers(int): Number of answers
		nClusters(int): Number of clusters
	
	Returns:
	    (
			np.ndarray: Mixture coefficient array
			np.ndarray: 3D array with dirichlet distribution of answer given their question and cluster
		)
	"""
    mixtureCoefs = np.random.dirichlet(np.ones(nClusters))
    probs = np.random.dirichlet(np.ones(nAnswers), (nClusters, nQuestions))
    return mixtureCoefs, probs


def generateData(model, nPeople):
    """ Generate answer to a survey provided with a model

	Args:
		model(tuple): model generated by 'generateModel'
		nPeople: Number of people who answered the survey
		
	Returns:
		(
			np.ndarray: nPeople*nQuestions array with answers to each question
			np.ndarray: nPeople array with the cluster associated to each person
		)
	"""
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    clusters = np.random.choice(nClusters, size=nPeople, p=mixtureCoefs)
    data = np.zeros((nPeople, nQuestions))
    answers = range(1, nAnswers + 1)
    for i in range(nPeople):
        for j in range(nQuestions):
            data[i, j] = np.random.choice(answers, p=probs[clusters[i], j, :])
    return np.asarray(data, 'uint'), clusters + 1


def drawModel(model):
    """ Plot distribution for each clusters

	Args:
		model(tuple): Model output by 'generateModel'
	"""
    plt.ion()
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    questions = np.arange(nQuestions)
    fig = plt.figure()
    for cluster in range(nClusters):
        plt.subplot(nClusters, 1, cluster + 1)
        colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k'])
        total = np.zeros(nQuestions)
        for answer in range(nAnswers):
            c = next(colors)
            d = probs[cluster][:, answer]
            plt.bar(questions, d, width=1, bottom=total, color=c)
            total += d
        plt.yticks(np.arange(0, 1.0001, 0.5))
        plt.ylim((0, 1))
        plt.title('Groupe ' + str(cluster + 1))
    plt.xlabel('Questions', fontsize=14)
    fig.tight_layout()
    plt.show()


def drawCluster(H):
    """ Plot distribution of latent variables
	H can be either a n*1 or n*nc matrix.

	The number n of rows of H is the number of persons
	* If H is a n*nc matrix, 
		the number of columns is the number of clusters
		and the components H[i,j] gives the likelihood 
		P(Person i belongs to cluster j)
	* If H is n*1 matrix, 
		the component H[i] gives the most probable cluster to 
		which person i belongs
	
	Args:
	    H(np.ndarray): (soft or hard) Clusters assignements
	"""
    plt.ion()
    if len(H.shape) >= 2:
        nPeople, nClusters = H.shape
    else:
        nPeople, = H.shape
        nClusters = np.max(H)
        Hnew = np.zeros((nPeople, nClusters))
        for i, j in zip(range(nPeople), H - 1):
            Hnew[i, j] = 1.0
        H = Hnew
    fig = plt.figure()
    people = np.arange(nPeople)
    colors = itertools.cycle(['b', 'r', 'g', 'c', 'm', 'y', 'k'])
    total = np.zeros(nPeople)
    for cluster in range(nClusters):
        c = next(colors)
        d = H[:, cluster]
        plt.bar(people, d, width=1, bottom=total, color=c)
        total += d
    plt.yticks(np.arange(0, 1.0001, 0.25))
    plt.ylim((0, 1))
    plt.title('Distribution of latent variables')
    plt.xlabel('person', fontsize=14)
    fig.tight_layout()
    plt.show()


<mask token>


def em(data, nClusters, nIterations):
    """ EM Algorithm
	the data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		nIterations(int): Number of iteration to run the EM updates
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    for t in range(nIterations):
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        return model, q


def normalize(v):
    """Normalize a 3D array

	Args:
		v(np.ndarray): 3D array
	
	Returns:
        np.ndarray: Normalized matrix
	"""
    dim = v.shape
    for i in range(0, dim[0] - 1):
        v[i, :, :] = (v[i, :, :].T / np.sum(v[i, :, :], 1)).T
    return v


def loglikelihood(model, data, q):
    """Computes and return the loglikelihood of a model M 
	given the data D and hidden variables distributions H (as returned by EM)

	Args:
		model: Mixture coefficient and visible variable distributions
		data: Matrix holding answers to each questions
		q: Latent cluster assignements
	
	Returns:
		float: Loglikelihood 
	"""
    ph, pvh = model
    nPeople, nQuestions = data.shape
    logL = 0
    for i in range(nPeople):
        answers = data[i, :]
        for k in range(nQuestions):
            logL += np.log(sum(pvh[:, k, int(answers[k] - 1)] * q[i, :].T))
    return logL


def em2(data, nClusters, epsilon):
    """ EM Algorithm with a stopping criterion. 
	The data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		epsilon(float): Threshold to detecte if the logLikelihood is still increasing
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    oldLogL = 0
    deltalogL = 100
    while deltalogL > epsilon:
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        optLogL = loglikelihood(model, data, q)
        print(optLogL)
        deltalogL = abs(oldLogL - optLogL)
        oldLogL = optLogL
        return model, q, optLogL

    def BIC(model, logL, nPeople):
        """
		Compute the Bayesian Information Criterion
		In order to  find the ideal number of clusters that produces accurate but not over-specific models, 
		we use a criterion that penalizes models with a large number of parameters. 
		The Bayesian Information Criterion (BIC) is one of these criteria. 
		While BIC has some theorical justifications for the exponential family of parameter distributions 
		whose categorical distributions are member of. We use it here as an empirical criterion.
		"""
        ph, pvh = model
        nClusters, nQuestions, nAnswers = pvh.shape
        nParams = nClusters - 1 + nClusters * nQuestions * (nAnswers - 1)
        bic = -2.0 * logL + nParams * np.log(nPeople)
        return bic

    def findBestModel(D):
        nPeople, nQuestions = D.shape
        bestBic = 0.1
        bestNClusters = 1
        for nClusters in range(1, 10):
            M, H, logL = em2(D, nClusters, 0.1)
            bic = BIC(M, logL, nPeople)
            if bic < bestBic:
                bestBic = bic
                bestModel = M
                bestNClusters = nClusters
        return bestModel, bestNClusters

<mask token>


def basicModel():
    """Returns toy model for debugging

	Returns:
		(
			np.ndarray: 1D array with clustering coefficients
			np.ndarray: 3D array with Dirichlet distribution for each answer regarding
			    a (cluster, question)
		)

	Returns
	-------
	(mixtureCoefs, probs) où
	- mixtureCoefs is a 1D vector holding the mixture coefficients
	- probs is a 3D array (nClusters, nQuestions, nAnswers)
		holding Dirichlet distribution associated to answers
		for each point (cluster, question)
	"""
    mixtureCoefs = np.array([0.4, 0.6])
    probs = np.array([[[0.2, 0.5, 0.3], [0.1, 0.2, 0.7], [0.7, 0.1, 0.2], [
        0.3, 0.3, 0.4]], [[0.7, 0.1, 0.2], [0.2, 0.6, 0.2], [0.2, 0.4, 0.4],
        [0.2, 0.6, 0.2]]])
    return mixtureCoefs, probs


def generateModel(nQuestions, nAnswers, nClusters):
    """ Generate a mixture model with Dirichlet distribution
	
	Args:
		nQuestion(int): Number of questions
		nAnswers(int): Number of answers
		nClusters(int): Number of clusters
	
	Returns:
	    (
			np.ndarray: Mixture coefficient array
			np.ndarray: 3D array with dirichlet distribution of answer given their question and cluster
		)
	"""
    mixtureCoefs = np.random.dirichlet(np.ones(nClusters))
    probs = np.random.dirichlet(np.ones(nAnswers), (nClusters, nQuestions))
    return mixtureCoefs, probs


def generateData(model, nPeople):
    """ Generate answer to a survey provided with a model

	Args:
		model(tuple): model generated by 'generateModel'
		nPeople: Number of people who answered the survey
		
	Returns:
		(
			np.ndarray: nPeople*nQuestions array with answers to each question
			np.ndarray: nPeople array with the cluster associated to each person
		)
	"""
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    clusters = np.random.choice(nClusters, size=nPeople, p=mixtureCoefs)
    data = np.zeros((nPeople, nQuestions))
    answers = range(1, nAnswers + 1)
    for i in range(nPeople):
        for j in range(nQuestions):
            data[i, j] = np.random.choice(answers, p=probs[clusters[i], j, :])
    return np.asarray(data, 'uint'), clusters + 1


def drawModel(model):
    """ Plot distribution for each clusters

	Args:
		model(tuple): Model output by 'generateModel'
	"""
    plt.ion()
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    questions = np.arange(nQuestions)
    fig = plt.figure()
    for cluster in range(nClusters):
        plt.subplot(nClusters, 1, cluster + 1)
        colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k'])
        total = np.zeros(nQuestions)
        for answer in range(nAnswers):
            c = next(colors)
            d = probs[cluster][:, answer]
            plt.bar(questions, d, width=1, bottom=total, color=c)
            total += d
        plt.yticks(np.arange(0, 1.0001, 0.5))
        plt.ylim((0, 1))
        plt.title('Groupe ' + str(cluster + 1))
    plt.xlabel('Questions', fontsize=14)
    fig.tight_layout()
    plt.show()


def drawCluster(H):
    """ Plot distribution of latent variables
	H can be either a n*1 or n*nc matrix.

	The number n of rows of H is the number of persons
	* If H is a n*nc matrix, 
		the number of columns is the number of clusters
		and the components H[i,j] gives the likelihood 
		P(Person i belongs to cluster j)
	* If H is n*1 matrix, 
		the component H[i] gives the most probable cluster to 
		which person i belongs
	
	Args:
	    H(np.ndarray): (soft or hard) Clusters assignements
	"""
    plt.ion()
    if len(H.shape) >= 2:
        nPeople, nClusters = H.shape
    else:
        nPeople, = H.shape
        nClusters = np.max(H)
        Hnew = np.zeros((nPeople, nClusters))
        for i, j in zip(range(nPeople), H - 1):
            Hnew[i, j] = 1.0
        H = Hnew
    fig = plt.figure()
    people = np.arange(nPeople)
    colors = itertools.cycle(['b', 'r', 'g', 'c', 'm', 'y', 'k'])
    total = np.zeros(nPeople)
    for cluster in range(nClusters):
        c = next(colors)
        d = H[:, cluster]
        plt.bar(people, d, width=1, bottom=total, color=c)
        total += d
    plt.yticks(np.arange(0, 1.0001, 0.25))
    plt.ylim((0, 1))
    plt.title('Distribution of latent variables')
    plt.xlabel('person', fontsize=14)
    fig.tight_layout()
    plt.show()


<mask token>


def permuteClusters(M, H, perm):
    """Permute clusters on a model M
	and latent cluster distributions H using
	cluster mapping perm

	Args:
		M: Model
		H: Latent clusters
		perm: Cluster mapping 
	
	Returns:
		tuple: (permuted model, permuted H)
	"""
    invperm = np.argsort(perm - 1)
    mixtureCoefs, probs = M
    mixtureCoefs = mixtureCoefs[invperm]
    probs = probs[invperm, :, :]
    H = H[:, invperm]
    return (mixtureCoefs, probs), H


def em(data, nClusters, nIterations):
    """ EM Algorithm
	the data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		nIterations(int): Number of iteration to run the EM updates
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    for t in range(nIterations):
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        return model, q


def normalize(v):
    """Normalize a 3D array

	Args:
		v(np.ndarray): 3D array
	
	Returns:
        np.ndarray: Normalized matrix
	"""
    dim = v.shape
    for i in range(0, dim[0] - 1):
        v[i, :, :] = (v[i, :, :].T / np.sum(v[i, :, :], 1)).T
    return v


def loglikelihood(model, data, q):
    """Computes and return the loglikelihood of a model M 
	given the data D and hidden variables distributions H (as returned by EM)

	Args:
		model: Mixture coefficient and visible variable distributions
		data: Matrix holding answers to each questions
		q: Latent cluster assignements
	
	Returns:
		float: Loglikelihood 
	"""
    ph, pvh = model
    nPeople, nQuestions = data.shape
    logL = 0
    for i in range(nPeople):
        answers = data[i, :]
        for k in range(nQuestions):
            logL += np.log(sum(pvh[:, k, int(answers[k] - 1)] * q[i, :].T))
    return logL


def em2(data, nClusters, epsilon):
    """ EM Algorithm with a stopping criterion. 
	The data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		epsilon(float): Threshold to detecte if the logLikelihood is still increasing
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    oldLogL = 0
    deltalogL = 100
    while deltalogL > epsilon:
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        optLogL = loglikelihood(model, data, q)
        print(optLogL)
        deltalogL = abs(oldLogL - optLogL)
        oldLogL = optLogL
        return model, q, optLogL

    def BIC(model, logL, nPeople):
        """
		Compute the Bayesian Information Criterion
		In order to  find the ideal number of clusters that produces accurate but not over-specific models, 
		we use a criterion that penalizes models with a large number of parameters. 
		The Bayesian Information Criterion (BIC) is one of these criteria. 
		While BIC has some theorical justifications for the exponential family of parameter distributions 
		whose categorical distributions are member of. We use it here as an empirical criterion.
		"""
        ph, pvh = model
        nClusters, nQuestions, nAnswers = pvh.shape
        nParams = nClusters - 1 + nClusters * nQuestions * (nAnswers - 1)
        bic = -2.0 * logL + nParams * np.log(nPeople)
        return bic

    def findBestModel(D):
        nPeople, nQuestions = D.shape
        bestBic = 0.1
        bestNClusters = 1
        for nClusters in range(1, 10):
            M, H, logL = em2(D, nClusters, 0.1)
            bic = BIC(M, logL, nPeople)
            if bic < bestBic:
                bestBic = bic
                bestModel = M
                bestNClusters = nClusters
        return bestModel, bestNClusters

import numpy as np
import matplotlib.pyplot as plt
import itertools
import sys


def basicModel():
    """Returns toy model for debugging

	Returns:
		(
			np.ndarray: 1D array with clustering coefficients
			np.ndarray: 3D array with Dirichlet distribution for each answer regarding
			    a (cluster, question)
		)

	Returns
	-------
	(mixtureCoefs, probs) où
	- mixtureCoefs is a 1D vector holding the mixture coefficients
	- probs is a 3D array (nClusters, nQuestions, nAnswers)
		holding Dirichlet distribution associated to answers
		for each point (cluster, question)
	"""
    mixtureCoefs = np.array([0.4, 0.6])
    probs = np.array([[[0.2, 0.5, 0.3], [0.1, 0.2, 0.7], [0.7, 0.1, 0.2], [
        0.3, 0.3, 0.4]], [[0.7, 0.1, 0.2], [0.2, 0.6, 0.2], [0.2, 0.4, 0.4],
        [0.2, 0.6, 0.2]]])
    return mixtureCoefs, probs


def generateModel(nQuestions, nAnswers, nClusters):
    """ Generate a mixture model with Dirichlet distribution
	
	Args:
		nQuestion(int): Number of questions
		nAnswers(int): Number of answers
		nClusters(int): Number of clusters
	
	Returns:
	    (
			np.ndarray: Mixture coefficient array
			np.ndarray: 3D array with dirichlet distribution of answer given their question and cluster
		)
	"""
    mixtureCoefs = np.random.dirichlet(np.ones(nClusters))
    probs = np.random.dirichlet(np.ones(nAnswers), (nClusters, nQuestions))
    return mixtureCoefs, probs


def generateData(model, nPeople):
    """ Generate answer to a survey provided with a model

	Args:
		model(tuple): model generated by 'generateModel'
		nPeople: Number of people who answered the survey
		
	Returns:
		(
			np.ndarray: nPeople*nQuestions array with answers to each question
			np.ndarray: nPeople array with the cluster associated to each person
		)
	"""
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    clusters = np.random.choice(nClusters, size=nPeople, p=mixtureCoefs)
    data = np.zeros((nPeople, nQuestions))
    answers = range(1, nAnswers + 1)
    for i in range(nPeople):
        for j in range(nQuestions):
            data[i, j] = np.random.choice(answers, p=probs[clusters[i], j, :])
    return np.asarray(data, 'uint'), clusters + 1


def drawModel(model):
    """ Plot distribution for each clusters

	Args:
		model(tuple): Model output by 'generateModel'
	"""
    plt.ion()
    mixtureCoefs, probs = model
    nClusters, nQuestions, nAnswers = probs.shape
    questions = np.arange(nQuestions)
    fig = plt.figure()
    for cluster in range(nClusters):
        plt.subplot(nClusters, 1, cluster + 1)
        colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k'])
        total = np.zeros(nQuestions)
        for answer in range(nAnswers):
            c = next(colors)
            d = probs[cluster][:, answer]
            plt.bar(questions, d, width=1, bottom=total, color=c)
            total += d
        plt.yticks(np.arange(0, 1.0001, 0.5))
        plt.ylim((0, 1))
        plt.title('Groupe ' + str(cluster + 1))
    plt.xlabel('Questions', fontsize=14)
    fig.tight_layout()
    plt.show()


def drawCluster(H):
    """ Plot distribution of latent variables
	H can be either a n*1 or n*nc matrix.

	The number n of rows of H is the number of persons
	* If H is a n*nc matrix, 
		the number of columns is the number of clusters
		and the components H[i,j] gives the likelihood 
		P(Person i belongs to cluster j)
	* If H is n*1 matrix, 
		the component H[i] gives the most probable cluster to 
		which person i belongs
	
	Args:
	    H(np.ndarray): (soft or hard) Clusters assignements
	"""
    plt.ion()
    if len(H.shape) >= 2:
        nPeople, nClusters = H.shape
    else:
        nPeople, = H.shape
        nClusters = np.max(H)
        Hnew = np.zeros((nPeople, nClusters))
        for i, j in zip(range(nPeople), H - 1):
            Hnew[i, j] = 1.0
        H = Hnew
    fig = plt.figure()
    people = np.arange(nPeople)
    colors = itertools.cycle(['b', 'r', 'g', 'c', 'm', 'y', 'k'])
    total = np.zeros(nPeople)
    for cluster in range(nClusters):
        c = next(colors)
        d = H[:, cluster]
        plt.bar(people, d, width=1, bottom=total, color=c)
        total += d
    plt.yticks(np.arange(0, 1.0001, 0.25))
    plt.ylim((0, 1))
    plt.title('Distribution of latent variables')
    plt.xlabel('person', fontsize=14)
    fig.tight_layout()
    plt.show()


def compareClusters(C1, C2):
    """ Renvoie le degré de superposition de deux clusters C1 et C2 
		à une permutation près des clusters
		
		C1 et C2 donnent les clusters associés à chaque exemple.
		La ligne i de C1 et C2 correspond à l'exemple i.
		C1 et C2 peuvent être soit un vecteur soit une matrice.
		* Si C est un vecteur, C[i] contient le cluster auquel i appartient.
		* Si C est une matrice, le nombre de colonnes de C correspond au nombre de clusters
		et C[i,j] contient la probabilité que l'exemple i appartient au cluster j
		Dans ce cas le cluster le plus probable est associé à i
		
		La fonction renvoie (bestOverlap, bestPerm) où
		* bestOverlap est le degré de superposition de C1 et C2
		* bestPerm donne la correspondance des clusters de C2 avec ceux de C1
		bestPerm[n° du cluster de C2] = n° du cluster de C1 correspondant """
    if C1.shape[0] != C2.shape[0]:
        raise 'C1 and C2 must have the same number of lines'
    if len(C1.shape) >= 2 and C1.shape[1] > 1:
        C1 = np.argmax(C1, 1) + 1
    if len(C2.shape) >= 2 and C2.shape[1] > 1:
        C2 = np.argmax(C2, 1) + 1
    nClusters = max(np.max(C1), np.max(C2))
    perms = itertools.permutations(range(1, nClusters + 1))
    bestOverlap = 0.0
    bestPerm = None
    for perm in perms:
        p = np.array(perm)
        permC2 = p[C2 - 1]
        overlap = np.sum(C1 == permC2) * 1.0 / np.size(C1)
        if overlap > bestOverlap:
            bestOverlap = overlap
            bestPerm = p
    return bestOverlap, bestPerm


def permuteClusters(M, H, perm):
    """Permute clusters on a model M
	and latent cluster distributions H using
	cluster mapping perm

	Args:
		M: Model
		H: Latent clusters
		perm: Cluster mapping 
	
	Returns:
		tuple: (permuted model, permuted H)
	"""
    invperm = np.argsort(perm - 1)
    mixtureCoefs, probs = M
    mixtureCoefs = mixtureCoefs[invperm]
    probs = probs[invperm, :, :]
    H = H[:, invperm]
    return (mixtureCoefs, probs), H


def em(data, nClusters, nIterations):
    """ EM Algorithm
	the data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		nIterations(int): Number of iteration to run the EM updates
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    for t in range(nIterations):
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        return model, q


def normalize(v):
    """Normalize a 3D array

	Args:
		v(np.ndarray): 3D array
	
	Returns:
        np.ndarray: Normalized matrix
	"""
    dim = v.shape
    for i in range(0, dim[0] - 1):
        v[i, :, :] = (v[i, :, :].T / np.sum(v[i, :, :], 1)).T
    return v


def loglikelihood(model, data, q):
    """Computes and return the loglikelihood of a model M 
	given the data D and hidden variables distributions H (as returned by EM)

	Args:
		model: Mixture coefficient and visible variable distributions
		data: Matrix holding answers to each questions
		q: Latent cluster assignements
	
	Returns:
		float: Loglikelihood 
	"""
    ph, pvh = model
    nPeople, nQuestions = data.shape
    logL = 0
    for i in range(nPeople):
        answers = data[i, :]
        for k in range(nQuestions):
            logL += np.log(sum(pvh[:, k, int(answers[k] - 1)] * q[i, :].T))
    return logL


def em2(data, nClusters, epsilon):
    """ EM Algorithm with a stopping criterion. 
	The data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		epsilon(float): Threshold to detecte if the logLikelihood is still increasing
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	"""
    nPeople, nQuestions = data.shape
    nAnswers = int(np.max(data))
    q = np.zeros((nPeople, nClusters))
    ph = np.ones(nClusters) / nClusters
    pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range
        (nClusters)])
    pvh = normalize(pvh)
    oldLogL = 0
    deltalogL = 100
    while deltalogL > epsilon:
        for i in range(nPeople):
            answers = data[i, :]
            for c in range(nClusters):
                q[i, c] = np.log(ph[c])
                for j in range(nQuestions):
                    q[i, c] += np.log(pvh[c, j, int(answers[j] - 1)])
            q[i, :] -= np.max(q[i, :])
            q[i, :] = np.exp(q[i, :])
            q[i, :] /= np.sum(q[i, :])
        weights = np.sum(q, 0)
        ph = weights / np.sum(weights)
        for cluster in range(nClusters):
            for j in range(nQuestions):
                for k in range(1, nAnswers + 1):
                    s = 0
                    for i in range(nPeople):
                        if data[i, j] == k:
                            s += q[i, cluster]
                    pvh[cluster, j, k - 1] = s / weights[cluster]
                pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
        model = ph, pvh
        optLogL = loglikelihood(model, data, q)
        print(optLogL)
        deltalogL = abs(oldLogL - optLogL)
        oldLogL = optLogL
        return model, q, optLogL

    def BIC(model, logL, nPeople):
        """
		Compute the Bayesian Information Criterion
		In order to  find the ideal number of clusters that produces accurate but not over-specific models, 
		we use a criterion that penalizes models with a large number of parameters. 
		The Bayesian Information Criterion (BIC) is one of these criteria. 
		While BIC has some theorical justifications for the exponential family of parameter distributions 
		whose categorical distributions are member of. We use it here as an empirical criterion.
		"""
        ph, pvh = model
        nClusters, nQuestions, nAnswers = pvh.shape
        nParams = nClusters - 1 + nClusters * nQuestions * (nAnswers - 1)
        bic = -2.0 * logL + nParams * np.log(nPeople)
        return bic

    def findBestModel(D):
        nPeople, nQuestions = D.shape
        bestBic = 0.1
        bestNClusters = 1
        for nClusters in range(1, 10):
            M, H, logL = em2(D, nClusters, 0.1)
            bic = BIC(M, logL, nPeople)
            if bic < bestBic:
                bestBic = bic
                bestModel = M
                bestNClusters = nClusters
        return bestModel, bestNClusters

#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import itertools
import sys

def basicModel():
	'''Returns toy model for debugging

	Returns:
		(
			np.ndarray: 1D array with clustering coefficients
			np.ndarray: 3D array with Dirichlet distribution for each answer regarding
			    a (cluster, question)
		)

	Returns
	-------
	(mixtureCoefs, probs) où
	- mixtureCoefs is a 1D vector holding the mixture coefficients
	- probs is a 3D array (nClusters, nQuestions, nAnswers)
		holding Dirichlet distribution associated to answers
		for each point (cluster, question)
	'''
	
	# Mixture Cluster Coefficients
	mixtureCoefs = np.array([0.4,0.6])
	
	# Probabilities of 3 answers to four possible questions given their clusters
	probs = np.array([
	# P(answer = j | question = i, cluster = 1)
	[
		[0.2 , 0.5, 0.3],
		[0.1 , 0.2, 0.7 ],
		[0.7 , 0.1, 0.2 ],
		[0.3 , 0.3, 0.4 ]
	],
	# P(answer = j | question = i, cluster = 2)
	[
		[0.7, 0.1, 0.2 ],
		[0.2, 0.6, 0.2 ],
		[0.2, 0.4, 0.4 ],
		[0.2, 0.6, 0.2 ]
	]])
	
	# Model
	return (mixtureCoefs, probs)

def generateModel(nQuestions, nAnswers, nClusters):
	''' Generate a mixture model with Dirichlet distribution
	
	Args:
		nQuestion(int): Number of questions
		nAnswers(int): Number of answers
		nClusters(int): Number of clusters
	
	Returns:
	    (
			np.ndarray: Mixture coefficient array
			np.ndarray: 3D array with dirichlet distribution of answer given their question and cluster
		)
	'''
		
	mixtureCoefs = np.random.dirichlet(np.ones(nClusters))
	probs = np.random.dirichlet(np.ones(nAnswers), (nClusters, nQuestions))
	return (mixtureCoefs, probs)
	
def generateData(model, nPeople):
	''' Generate answer to a survey provided with a model

	Args:
		model(tuple): model generated by 'generateModel'
		nPeople: Number of people who answered the survey
		
	Returns:
		(
			np.ndarray: nPeople*nQuestions array with answers to each question
			np.ndarray: nPeople array with the cluster associated to each person
		)
	'''
	
	mixtureCoefs, probs = model
	nClusters, nQuestions, nAnswers = probs.shape
	clusters = np.random.choice(nClusters, size = nPeople, p = mixtureCoefs)
	data = np.zeros((nPeople, nQuestions))
	answers = range(1, nAnswers+1)
	for i in range(nPeople):
		for j in range(nQuestions):
			data[i,j] = np.random.choice(answers, p=probs[clusters[i],j,:])	
	return (np.asarray(data, 'uint'), clusters + 1)

def drawModel(model):
	''' Plot distribution for each clusters

	Args:
		model(tuple): Model output by 'generateModel'
	'''
	
	plt.ion()
	mixtureCoefs, probs = model
	nClusters, nQuestions, nAnswers = probs.shape	
	questions = np.arange(nQuestions)
	fig = plt.figure()
	for cluster in range(nClusters):
		plt.subplot(nClusters, 1, cluster+1)	
		colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k'])
		total = np.zeros(nQuestions)		
		for answer in range(nAnswers):
			c = next(colors)
			d = probs[cluster][:,answer]
			plt.bar(questions, d, width=1, bottom = total, color=c)
			total += d
		plt.yticks(np.arange(0, 1.0001, 0.5))
		plt.ylim((0,1))
		plt.title('Groupe ' + str(cluster+1))

	plt.xlabel('Questions', fontsize=14)
	fig.tight_layout()		
	plt.show()

def drawCluster(H):
	''' Plot distribution of latent variables
	H can be either a n*1 or n*nc matrix.

	The number n of rows of H is the number of persons
	* If H is a n*nc matrix, 
		the number of columns is the number of clusters
		and the components H[i,j] gives the likelihood 
		P(Person i belongs to cluster j)
	* If H is n*1 matrix, 
		the component H[i] gives the most probable cluster to 
		which person i belongs
	
	Args:
	    H(np.ndarray): (soft or hard) Clusters assignements
	'''
	
	plt.ion()
	if(len(H.shape) >= 2):
		nPeople, nClusters = H.shape
	else:
		nPeople, = H.shape
		nClusters = np.max(H)
		Hnew = np.zeros((nPeople, nClusters))
		for (i,j) in zip(range(nPeople), H-1):
			Hnew[i,j] = 1.
		H = Hnew
	fig = plt.figure()
	people = np.arange(nPeople)
	colors = itertools.cycle(['b', 'r', 'g', 'c', 'm', 'y', 'k'])
	total = np.zeros(nPeople)		
	for cluster in range(nClusters):
		c = next(colors)
		d = H[:,cluster]
		plt.bar(people, d, width=1, bottom = total, color=c)
		total += d
	plt.yticks(np.arange(0, 1.0001, 0.25))
	plt.ylim((0,1))
	plt.title('Distribution of latent variables')
	plt.xlabel('person', fontsize=14)
	fig.tight_layout()		
	plt.show()

def compareClusters(C1,C2):
	''' Renvoie le degré de superposition de deux clusters C1 et C2 
		à une permutation près des clusters
		
		C1 et C2 donnent les clusters associés à chaque exemple.
		La ligne i de C1 et C2 correspond à l'exemple i.
		C1 et C2 peuvent être soit un vecteur soit une matrice.
		* Si C est un vecteur, C[i] contient le cluster auquel i appartient.
		* Si C est une matrice, le nombre de colonnes de C correspond au nombre de clusters
		et C[i,j] contient la probabilité que l'exemple i appartient au cluster j
		Dans ce cas le cluster le plus probable est associé à i
		
		La fonction renvoie (bestOverlap, bestPerm) où
		* bestOverlap est le degré de superposition de C1 et C2
		* bestPerm donne la correspondance des clusters de C2 avec ceux de C1
		bestPerm[n° du cluster de C2] = n° du cluster de C1 correspondant '''

	if(C1.shape[0] != C2.shape[0]):
		raise "C1 and C2 must have the same number of lines"
			
	if(len(C1.shape) >= 2 and C1.shape[1] > 1):
		C1 = np.argmax(C1,1)+1
	if(len(C2.shape) >= 2 and C2.shape[1] > 1):
		C2 = np.argmax(C2,1)+1
			
	nClusters = max(np.max(C1), np.max(C2))
	perms = itertools.permutations(range(1,nClusters+1))
	bestOverlap = 0.
	bestPerm = None
	for perm in perms:
		p = np.array(perm)
		permC2 = p[C2-1]
		overlap = np.sum(C1 == permC2) * 1. / np.size(C1)
		if(overlap > bestOverlap):
			bestOverlap = overlap
			bestPerm = p
	return (bestOverlap, bestPerm)
	
def permuteClusters(M, H, perm):
	'''Permute clusters on a model M
	and latent cluster distributions H using
	cluster mapping perm

	Args:
		M: Model
		H: Latent clusters
		perm: Cluster mapping 
	
	Returns:
		tuple: (permuted model, permuted H)
	'''
	invperm = np.argsort(perm-1)
	mixtureCoefs, probs = M
	mixtureCoefs = mixtureCoefs[invperm]
	probs = probs[invperm, :, :]
	H = H[:, invperm]
	return ((mixtureCoefs, probs), H)

###############
# A compléter #
###############

def em(data, nClusters, nIterations):
	''' EM Algorithm
	the data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		nIterations(int): Number of iteration to run the EM updates
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	'''

	nPeople, nQuestions = data.shape
	nAnswers = int(np.max(data))

	# q[i,j] = P(people i belong to Cluster j)
	# P(Hi|Vi = vi,Θ = θ) ~ qi
	# We store those distributions in an array
	q = np.zeros((nPeople, nClusters))

    # Mixture coefficients ph
	# Parameters of the distribution of the hidden variables
	# P(H) = ph(i) = P(Cluster = i) 
	ph = np.ones(nClusters)/nClusters

	# Probability of the choice of a visible variable
	# knowing that a given latent variable was chosen
	# P(V|H,θ) = pvh(i,j,k) = P(Answer = A(i,j) | Cluster = k)
	# Answer i,j is the the answer given by person i to question j
	# We could roughly model A by couple of variable (P,Q) i.e P answered Q
	pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range(nClusters)])
	# Parameters of different clusters distribution are all equal
	# pvh = np.array([np.ones((nQuestions, nAnswers)) for i in range(nClusters)])
	pvh = normalize(pvh)
    
	for t in range(nIterations):
	    # Step E 
	    # First We focus only on the distribution q to maximize the Likelihood 
		for i in range(nPeople):
			answers = data[i,:]
			for c in range(nClusters):
				# We use log since the probabilities might have small values
				# We use Bayes rule i.e H|V,θ = V|(H,θ) * H|θ / (V,θ) 
				# First we add the log of our prior H
				q[i,c] = np.log(ph[c])
                
				# Then we add the logs of Vi|H
				for j in range(nQuestions):
					q[i,c] += np.log(pvh[c,j,int(answers[j]-1)])
				
			q[i,:] -= np.max(q[i,:])
			q[i,:] = np.exp(q[i,:])
			# normalize
			q[i,:] /= np.sum(q[i,:])
        
		# Step M
		# The distribution of our cluster is the average of our the likelihoods distribution qi
		weights = np.sum(q,0)
		ph = weights/np.sum(weights)

		for cluster in range(nClusters):
			for j in range(nQuestions):
				for k in range(1, nAnswers + 1):
					s = 0
					#total = 0.
					for i in range(nPeople):
						if (data[i,j] == k):
							s+= q[i,cluster]
						#total += q[i,cluster]
					#pvh[cluster, j, k-1] = s/total
					pvh[cluster, j, k-1] = s/weights[cluster]
                
				# Normalization 
				pvh[cluster, j, :] /= sum(pvh[cluster, j, :])

		model = (ph, pvh)
		return (model, q)


def normalize(v):
	"""Normalize a 3D array

	Args:
		v(np.ndarray): 3D array
	
	Returns:
        np.ndarray: Normalized matrix
	"""
	dim = v.shape 
	for i in range(0, dim[0]-1):
		v[i,:,:] = (v[i,:,:].T/np.sum(v[i,:,:],1)).T

	return v

def loglikelihood(model, data, q):
	"""Computes and return the loglikelihood of a model M 
	given the data D and hidden variables distributions H (as returned by EM)

	Args:
		model: Mixture coefficient and visible variable distributions
		data: Matrix holding answers to each questions
		q: Latent cluster assignements
	
	Returns:
		float: Loglikelihood 
	"""
	ph, pvh = model
	nPeople, nQuestions = data.shape
	logL = 0
	for i in range(nPeople):
		answers = data[i,:]
		for k in range(nQuestions):
			logL += np.log(sum(pvh[:, k, int(answers[k] - 1)] * q[i,:].T))
	return logL

def em2(data, nClusters, epsilon):
	''' EM Algorithm with a stopping criterion. 
	The data used here comes from 'generateData'

	Args:
		data(int): nPeople*nQuestion matrix with at each row the answer to each question for one person
		nClusters(int): Number of Clusters
		epsilon(float): Threshold to detecte if the logLikelihood is still increasing
	
	Returns:
		tuple: (model,H) with:
		   (
			   model: Learnt model i.e the mixture coefficient and the visible variable (probability)
			   H: Cluster assignements i.e the probability for each person to belong to a cluster
		   )
	- model est le modèle appris au même format que les modèles produits pare generateModel
	- H est une matrice (nombre de personnes, nombre de clusters) dont les lignes
	correspondent aux distribution de cluster de chaque personne
	'''

	nPeople, nQuestions = data.shape
	nAnswers = int(np.max(data))

	# Each latent variable Hi follows a probability of distribution qi
	# q[i,j] = P(personne i belongs to Cluster j)
	# We store those distribution in an array
	q = np.zeros((nPeople, nClusters))

    # Mixture coefficients ph
	# Parameters of the distribution of the hidden variables
	# pc = ph(i) = P(Cluster = i) 
	ph = np.ones(nClusters)/nClusters

	# Probability of the choice of a visible variable
	# knowing that a given latent variable was chosen
	# pvh(i,j,k) = P(Reponse = A(i,j) | Cluster = k)
	pvh = np.array([np.random.random((nQuestions, nAnswers)) for i in range(nClusters)])
	# Parameters of different clusters distribution are all equal
	#pvh = np.array([np.ones((nQuestions, nAnswers)) for i in range(nClusters)])
	pvh = normalize(pvh)

	# Threshold epsilon
	oldLogL = 0;
	deltalogL = 100;
    
	while deltalogL > epsilon:

		# Step E 
		# First We focus only on the distribution q to maximize the MLE
		for i in range(nPeople):
			answers = data[i,:]
			for c in range(nClusters):
				# We use log since the probabilities might have small values 
				# First we add the log of our prior 
				q[i,c] = np.log(ph[c])
				
				# Then we add the log-likelihood of H|V
				for j in range(nQuestions):
					q[i,c] += np.log(pvh[c,j,int(answers[j]-1)])
				
			q[i,:] -= np.max(q[i,:])
			q[i,:] = np.exp(q[i,:])
			# normalize
			q[i,:] /= np.sum(q[i,:])
		
		# Step M
		weights = np.sum(q,0)
		ph = weights/np.sum(weights)

		for cluster in range(nClusters):
			for j in range(nQuestions):
				for k in range(1, nAnswers + 1):
					s = 0
					#total = 0.
					for i in range(nPeople):
						if (data[i,j] == k):
							s+= q[i,cluster]
						#total += q[i,cluster]
					#pvh[cluster, j, k-1] = s/total
					pvh[cluster, j, k-1] = s/weights[cluster]
				
				# Normalization 
				pvh[cluster, j, :] /= sum(pvh[cluster, j, :])
		
		model = (ph, pvh)
		optLogL = loglikelihood(model, data, q)
		print(optLogL)
		deltalogL = abs(oldLogL - optLogL)
		oldLogL = optLogL
		return (model, q, optLogL)

	def BIC(model, logL, nPeople):
		"""
		Compute the Bayesian Information Criterion
		In order to  find the ideal number of clusters that produces accurate but not over-specific models, 
		we use a criterion that penalizes models with a large number of parameters. 
		The Bayesian Information Criterion (BIC) is one of these criteria. 
		While BIC has some theorical justifications for the exponential family of parameter distributions 
		whose categorical distributions are member of. We use it here as an empirical criterion.
		"""
		ph, pvh = model
		nClusters, nQuestions, nAnswers = pvh.shape

		nParams = (nClusters - 1) + nClusters * nQuestions *(nAnswers - 1)
		bic = -2. * logL + nParams * np.log(nPeople)
		return bic
	
	def findBestModel(D):
		nPeople, nQuestions = D.shape
		bestBic = 1E-1
		bestNClusters = 1
		for nClusters in range(1,10):
			M, H, logL = em2(D, nClusters, 1E-1)
			bic = BIC(M, logL, nPeople)
			if(bic < bestBic):
				bestBic = bic
				bestModel = M
				bestNClusters = nClusters
		
		return bestModel, bestNClusters




