mi@0: #!/usr/bin/python
mi@0: #
mi@0: # Copyright (C) Christian Thurau, 2010. 
mi@0: # Licensed under the GNU General Public License (GPL). 
mi@0: # http://www.gnu.org/licenses/gpl.txt
mi@0: """
mi@0: PyMF several distance functions
mi@0: 
mi@0:     kl_divergence(): KL Divergence
mi@0:     l1_distance(): L1 distance
mi@0:     l2_distance(): L2 distance
mi@0:     cosine_distance(): Cosine distance 
mi@0:     pdist(): Pairwise distance computation
mi@0:     vq(): Vector quantization
mi@0:     
mi@0: """
mi@0: 
mi@0: 
mi@0: import numpy as np
mi@0: import scipy.sparse
mi@0: 
mi@0: __all__ = ["abs_cosine_distance", "kl_divergence", "l1_distance", "l2_distance", 
mi@0:            "weighted_abs_cosine_distance","cosine_distance","vq", "pdist"]
mi@0: 
mi@0: def kl_divergence(d, vec):    
mi@0:     b = vec*(1/d)    
mi@0:     b = np.where(b>0, np.log(b),0)    
mi@0:     b = vec * b
mi@0:     b = np.sum(b - vec + d, axis=0).reshape((-1))            
mi@0:     return b
mi@0:         
mi@0: def l1_distance(d, vec):            
mi@0:     ret_val = np.sum(np.abs(d - vec), axis=0)                    
mi@0:     ret_val = ret_val.reshape((-1))                        
mi@0:     return ret_val
mi@0:     
mi@0: def sparse_l2_distance(d, vec):
mi@0:     # compute the norm of d
mi@0:     nd = (d.multiply(d)).sum(axis=0)
mi@0:     nv = (vec.multiply(vec)).sum(axis=0)
mi@0:     ret_val = nd + nv -  2.0*(d.T * vec).T
mi@0:     return np.sqrt(ret_val)
mi@0: 
mi@0: def approx_l2_distance(d, vec):
mi@0:     # Use random projections to approximate the conventional l2 distance
mi@0:     k = np.round(np.log(d.shape[0]))
mi@0:     #k = d.shape[0]
mi@0:     R = np.random.randn(k, d.shape[0])
mi@0:     R = R / np.sqrt((R**2).sum(axis=0))  
mi@0:     A = np.dot(R,d)
mi@0:     B = np.dot(R, vec)
mi@0:     ret_val = np.sum( (A - B)**2, axis=0)                    
mi@0:     ret_val = np.sqrt(R.shape[1]/R.shape[0]) * np.sqrt(ret_val)
mi@0:     ret_val = ret_val.reshape((-1))    
mi@0:     return ret_val
mi@0:     
mi@0: def l2_distance(d, vec):    
mi@0:     if scipy.sparse.issparse(d):
mi@0:         ret_val = sparse_l2_distance(d, vec)
mi@0:     else:
mi@0:         ret_val = np.sqrt(((d[:,:] - vec)**2).sum(axis=0))
mi@0:             
mi@0:     return ret_val.reshape((-1))        
mi@0: 
mi@0: def l2_distance_new(d,vec):
mi@0:     # compute the norm of d
mi@0:     nd = (d**2).sum(axis=0)
mi@0:     nv = (vec**2).sum(axis=0)
mi@0:     ret_val = nd + nv - 2.0*np.dot(d.T,vec.reshape((-1,1))).T
mi@0: 
mi@0:     return np.sqrt(ret_val)
mi@0:     
mi@0: def cosine_distance(d, vec):
mi@0:     tmp = np.dot(np.transpose(d), vec)
mi@0:     a = np.sqrt(np.sum(d**2, axis=0))
mi@0:     b = np.sqrt(np.sum(vec**2))
mi@0:     k = (a*b).reshape(-1) + (10**-9)
mi@0:     
mi@0:     # compute distance
mi@0:     ret_val = 1.0 - tmp/k
mi@0:     
mi@0:     return ret_val.reshape((-1))
mi@0: 
mi@0: def abs_cosine_distance(d, vec, weighted=False):
mi@0:     if scipy.sparse.issparse(d):
mi@0:         tmp = np.array((d.T * vec).todense(), dtype=np.float32).reshape(-1)
mi@0:         a = np.sqrt(np.array(d.multiply(d).sum(axis=0), dtype=np.float32).reshape(-1))
mi@0:         b = np.sqrt(np.array(vec.multiply(vec).sum(axis=0), dtype=np.float32).reshape(-1))
mi@0:     else:
mi@0:         tmp = np.dot(np.transpose(d), vec).reshape(-1)
mi@0:         a = np.sqrt(np.sum(d**2, axis=0)).reshape(-1)
mi@0:         b = np.sqrt(np.sum(vec**2)).reshape(-1)
mi@0:         
mi@0:     k = (a*b).reshape(-1) + 10**-9
mi@0:     
mi@0:     # compute distance
mi@0:     ret_val = 1.0 - np.abs(tmp/k)    
mi@0:     
mi@0:     if weighted:    
mi@0:         ret_val = ret_val * a        
mi@0:     return ret_val.reshape((-1))
mi@0: 
mi@0: def weighted_abs_cosine_distance(d, vec):
mi@0:     ret_val = abs_cosine_distance(d, vec, weighted=True)        
mi@0:     return ret_val
mi@0: 
mi@0: def pdist(A, B, metric='l2' ):
mi@0:     # compute pairwise distance between a data matrix A (d x n) and B (d x m).
mi@0:     # Returns a distance matrix d (n x m).
mi@0:     d = np.zeros((A.shape[1], B.shape[1]))
mi@0:     if A.shape[1] <= B.shape[1]:
mi@0:         for aidx in xrange(A.shape[1]):
mi@0:             if metric == 'l2':
mi@0:                 d[aidx:aidx+1,:] = l2_distance(B[:,:], A[:,aidx:aidx+1]).reshape((1,-1))
mi@0:             if metric == 'l1':
mi@0:                 d[aidx:aidx+1,:] = l1_distance(B[:,:], A[:,aidx:aidx+1]).reshape((1,-1))
mi@0:     else:
mi@0:         for bidx in xrange(B.shape[1]):
mi@0:             if metric == 'l2':
mi@0:                 d[:, bidx:bidx+1] = l2_distance(A[:,:], B[:,bidx:bidx+1]).reshape((-1,1))
mi@0:             if metric == 'l1':
mi@0:                 d[:, bidx:bidx+1] = l1_distance(A[:,:], B[:,bidx:bidx+1]).reshape((-1,1))               
mi@0:     
mi@0:     return d
mi@0: 
mi@0: def vq(A, B, metric='l2'):
mi@0:     # assigns data samples in B to cluster centers A and
mi@0:     # returns an index list [assume n column vectors, d x n]
mi@0:     assigned = np.argmin(pdist(A,B, metric=metric), axis=0)
mi@0:     return assigned