view utils/gmmdist.py @ 19:890cfe424f4a tip

added annotations

#!/usr/bin/env python
# encoding: utf-8
"""
gmmdist.py

Created by George Fazekas on 2014-07-03.
Copyright (c) 2014 . All rights reserved.
"""

import sys,os

from numpy import sum,isnan,isinf,vstack
from numpy.random import rand
import numpy as np
from numpy import log, power, pi, exp, transpose, zeros, log, ones, dot
from sklearn.mixture import GMM
from sklearn.metrics.pairwise import pairwise_distances
from scipy.linalg import *
from pylab import plt
from scipy.stats import norm
#from gmmplot import plot_gmm
from scipy.io import savemat,loadmat
from numpy import trace
from numpy.linalg import det, inv, eigvals

FLOAT_MAX = np.finfo(np.float64).max

def is_pos_def(x):
	'''Check if matrix is positive definite.'''
	return np.all(eigvals(x) > 0)
	
def skl_models(g1,g2):
	'''Wrapper function with error checking and adaptation to the GmmDistance/GaussianComponent classes.
	This function compares two Gaussian mixture models with and equal number of components and full covariance matrices.
	Covariance matrices must be positive definite.
	'''
	m1,m2 = g1.means,g2.means
	v1,v2 = g1.covars.swapaxes(0,2), g2.covars.swapaxes(0,2)
	w1,w2 = g1.weights[:,np.newaxis],g2.weights[:,np.newaxis]
	assert m1.shape[1] > 1, "The minimum number of features is 2."
	assert w1.shape == w2.shape, "Models must have the same number of components"
	# print 'v1, v2', v1.shape, v2.shape
	# assert (is_pos_def(v1) and is_pos_def(v2)) == True, "Covariance matrices must be positive definite."
	d = skl_gmm(m1,m2,v1,v2,w1,w2)
	if isnan(d): #or isinf(d):
		return FLOAT_MAX
	return d

# def kldiv_full(m0,m1,s0,s1):
# 	'''Naive (unoptimised) implementation of the KL divergence between two single Gaussians with fully defined covariances (s).'''
# 	return 0.5*(np.trace(s0/s1)+np.trace(s1/s0)+ np.dot( np.dot((m0-m1).T, np.linalg.inv(s0+s1)), (m0-m1)))-m0.shape[0]
# 	
# def skl_full(p0,p1):
# 	'''Symmetrised KL divergence computed from 2 KL divergences using mean( KL(p||q), KL(q||p) )'''
# 	d = (kldiv_full(p0.means,p0.covars,p1.means,p1.covars) + kldiv_full(p1.means,p1.covars,p0.means,p0.covars)) * 0.5
# 	d = sum(d)
# 	if isnan(d) : 
# 		return np.finfo(np.float64).max
# 	return d

def kldiv_full(m1,m2,s1,s2):
	m1,m2 = m1[:,None],m2[:,None]
	logdet1, logdet2 = log(det(s1)), log(det(s2))
	inv1, inv2 = inv(s1), inv(s2)
	m = m1-m2
	d = m.shape[0] # number of dimensions
	return 0.5 * ((logdet1-logdet2) + trace(dot(inv1,s2)) + dot(dot(m.T, inv1), m) - d)[0][0]

def _skl_full(m1,m2,s1,s2):
	m1,m2 = m1[:,None],m2[:,None]
	logdet1, logdet2 = log(det(s1)), log(det(s2))
	inv1, inv2 = inv(s1), inv(s2)
	m12 = m1-m2
	m21 = m2-m1
	d = m12.shape[0] # number of dimensions
	kl12 = 0.5 * ((logdet1-logdet2) + trace(dot(inv1,s2)) + dot(dot(m12.T, inv1), m12) - d)[0][0]
	kl21 = 0.5 * ((logdet2-logdet1) + trace(dot(inv2,s1)) + dot(dot(m21.T, inv2), m21) - d)[0][0]
	return 0.5 * (kl12+kl21)

def skl_full(p1,p2):
	m1,m2 = p1.means,p2.means
	s1,s2 = p1.covars,p2.covars
	return _skl_full(m1,m2,s1,s2)	

	
		
def skl_gmm(m1,m2,v1,v2,w1,w2):
	'''Take the mean of	 KL(g1mm||gmm2) & KL(gmm2||gmm1) to symmetrise the divergence.'''
	return (abs(kldiv_gmm(m1,m2,v1,v2,w1,w2)) + abs(kldiv_gmm(m2,m1,v2,v1,w2,w1))) * 0.5
		
def kldiv_gmm(m1,m2,v1,v2,w1,w2):
	'''Low level implementation of variational approximation of KL divergence between Gaussian Mixture Models.
	See first: J. R. Hershey and P. A. Olsen. "Approximating the Kullback-Leibler Divergence Between Gaussian Mixture Models."
	IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Volume 4, pp. 317–320, April, 2007.
	Further theory and refinement: J. L. Durrieu, J. P. Thiran, F. Kelly. "Lower and Upper Bounds for Approximation 
	of the Kullback-Leibler Divergence Between Gaussian Mixture Models", ICASSP, 2012.
	This implementation is by George Fazekas, Centre for Digital Music, QMUL, London, UK.
	
	Inputs: 
	m(x) : mean vector of gmm(x)
	v(x) : covariance matrices of gmm(x)
	w(x) : weight vector of gmm(x)
	
	Ouptut:
	kl_full_12 : Kullback-Leibler divergence of the PDFs approximated by two Gaussian Mixture Models.
	
	This implementation is using a variational approximation rather than the conventional (and expensive) 
	Monte Carlo simulation based approach. See cited paper for details. The divergence is not symmetrised 
	and may be negative (unlike the closed form) or inf. In case the output is complex, the coavriance matrices
	do not fulfill required criteria, i.e. somehow badly formed, sigular, not positive definite etc...	  
	'''
	# TODO: consider dieling better with inf/-inf outcomes in the final distance computation 
	# - note: the max of the rows of kl12 are not always zero like that of kl11
	# - eliminate the need for swapaxes of the covariances
	
	n = m1.shape[0] # number of components
	d = m1.shape[1] # number of dimensions

	ixm = ones((n,1),dtype=int).T # selector of mean matrix components
	ixd = range(0,d*d,d+1) # indices of diagonal elements of DxD matrix
	t1 = v1.swapaxes(1,2).reshape(d,n*d) # concatenate gmm1 covariance matrices
	loopn = xrange(n)

	# step 1) precompute log(determinant()) of covariance matrices of gmm1
	logdet1 = zeros((n,1))
	for i in loopn :
		logdet1[i] = log(det(v1[:,:,i]))

	# step 2) compute reference kldiv between individual components of gmm1
	kl11 = zeros((n,n))
	for i in loopn :
		inv1 = inv(v1[:,:,i])
		mm1 = m1 - m1[i*ixm,:][0]
		b1 = dot(inv1,t1).swapaxes(0,1).reshape(n,power(d,2)).T
		kl11[:,i] = 0.5 * ( (logdet1[i]-d-logdet1)[:,0] + sum(b1[ixd,:],0).T + sum(dot(mm1,inv1) * mm1,1))
	# print kl11

	# step 3) precompute log(determinant()) of covariance matrices of gmm2
	logdet2 = zeros((n,1))
	for i in loopn :
		logdet2[i] = log(det(v2[:,:,i]))

	# step 4) compute pair-wise kldiv between components of gmm1 and gmm2
	kl12 = zeros((n,n))
	for i in loopn :
		inv2 = inv(v2[:,:,i])
		m12 = m1 - m2[i*ixm,:][0]
		b2 = dot(inv2,t1).swapaxes(0,1).reshape(n,power(d,2)).T
		kl12[:,i] = 0.5 * ( (logdet2[i]-d-logdet1)[:,0] + sum(b2[ixd,:],0).T + sum(dot(m12,inv2) * m12,1))
	# print kl12

	# step 5) compute the final variational distance between gmm1 and gmm2
	kl_full_12 = dot(w1.T, (log(sum(exp(-kl11)*w1,1))) - log(sum(exp(-kl12)*w2,1)))[0]
	# print "KL divergence between gmm1 || gmm2:", kl_full_12
	return kl_full_12
	

# models = loadmat("gmms.mat")
# # print models.keys()
# 
# X = models['X']
# Y = models['Y']
# 
# print "Data shape:"
# print X.shape
# print Y.shape
# 
# # # plot the fitted model
# # gmm1 = GMM(n_components = 3, covariance_type='full')
# # model1 = gmm1.fit(X)
# # plot_gmm(gmm1,X)
# # 
# # # plot the fitted model
# # gmm2 = GMM(n_components = 3, covariance_type='full')
# # model2 = gmm2.fit(Y)
# # plot_gmm(gmm2,Y)
# 
# # print "KL=",kldiv_full(gmm1.means_[0],gmm1.means_[1],gmm1.covars_[0],gmm1.covars_[1])
# # 
# # print "gmm1_covars:\n", gmm1.covars_, gmm1.covars_.shape
# 
# 
# # m1 = gmm1.means_
# # v1 = gmm1.covars_.swapaxes(0,2)
# # w1 = gmm1.weights_
# # 
# # m2 = gmm2.means_
# # v2 = gmm2.covars_.swapaxes(0,2)
# # w2 = gmm2.weights_
# 
# m1 = models['gmm1_means']
# v1 = models['gmm1_covars']
# w1 = models['gmm1_weights']
# 
# m2 = models['gmm2_means']
# v2 = models['gmm2_covars']
# w2 = models['gmm2_weights']
# 
# print "KL divergence between gmm1 || gmm2:", kldiv_gmm(m1,m2,v1,v2,w1,w2)
# print "KL divergence between gmm2 || gmm1:", kldiv_gmm(m2,m1,v2,v1,w2,w1)
# print "Symmetrised KL distance between gmm1 || gmm2:", skl_gmm(m1,m2,v1,v2,w1,w2)

def main():
	pass


if __name__ == '__main__':
	main()