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 LAESA
mi@0: """
mi@0: 
mi@0: 
mi@0: import scipy.sparse
mi@0: import numpy as np
mi@0: 
mi@0: from dist import *
mi@0: from sivm import SIVM
mi@0: 
mi@0: __all__ = ["LAESA"]
mi@0: 
mi@0: class LAESA(SIVM):
mi@0:     """      
mi@0:     LAESA(data, num_bases=4)
mi@0:     
mi@0:     
mi@0:     Simplex Volume Maximization. Factorize a data matrix into two matrices s.t.
mi@0:     F = | data - W*H | is minimal. H is restricted to convexity. W is iteratively
mi@0:     found by maximizing the volume of the resulting simplex (see [1]).
mi@0:     
mi@0:     Parameters
mi@0:     ----------
mi@0:     data : array_like, shape (_data_dimension, _num_samples)
mi@0:         the input data
mi@0:     num_bases: int, optional
mi@0:         Number of bases to compute (column rank of W and row rank of H).
mi@0:         4 (default)    
mi@0:     
mi@0:     Attributes
mi@0:     ----------
mi@0:     W : "data_dimension x num_bases" matrix of basis vectors
mi@0:     H : "num bases x num_samples" matrix of coefficients
mi@0:     ferr : frobenius norm (after calling .factorize())        
mi@0:     
mi@0:     Example
mi@0:     -------
mi@0:     Applying LAESA to some rather stupid data set:
mi@0:     
mi@0:     >>> import numpy as np
mi@0:     >>> data = np.array([[1.0, 0.0, 2.0], [0.0, 1.0, 1.0]])
mi@0:     >>> laesa_mdl = LAESA(data, num_bases=2)
mi@0:     >>> laesa_mdl.factorize()
mi@0:     
mi@0:     The basis vectors are now stored in laesa_mdl.W, the coefficients in laesa_mdl.H. 
mi@0:     To compute coefficients for an existing set of basis vectors simply    copy W 
mi@0:     to laesa_mdl.W, and set compute_w to False:
mi@0:     
mi@0:     >>> data = np.array([[1.5, 1.3], [1.2, 0.3]])
mi@0:     >>> W = np.array([[1.0, 0.0], [0.0, 1.0]])
mi@0:     >>> laesa_mdl = LAESA(data, num_bases=2)
mi@0:     >>> laesa_mdl.W = W
mi@0:     >>> laesa_mdl.factorize(niter=1, compute_w=False)
mi@0:     
mi@0:     The result is a set of coefficients laesa_mdl.H, s.t. data = W * laesa_mdl.H.
mi@0:     """
mi@0:     def update_w(self):    
mi@0:         # initialize some of the recursively updated distance measures     
mi@0:         self.init_sivm()
mi@0:         distiter = self._distance(self.select[-1])                
mi@0:         
mi@0:         for l in range(self._num_bases-1):                                        
mi@0:             d = self._distance(self.select[-1])                                
mi@0:         
mi@0:             # replace distances in distiter
mi@0:             distiter = np.where(d<distiter, d, distiter)
mi@0:             
mi@0:             # detect the next best data point                        
mi@0:             self.select.append(np.argmax(distiter))
mi@0:             self._logger.info('cur_nodes: ' + str(self.select))
mi@0: 
mi@0:         # sort indices, otherwise h5py won't work
mi@0:         self.W = self.data[:, np.sort(self.select)]
mi@0:         
mi@0:         # but "unsort" it again to keep the correct order
mi@0:         self.W = self.W[:, np.argsort(np.argsort(self.select))]    
mi@0:                            
mi@0:                    
mi@0: if __name__ == "__main__":
mi@0:     import doctest  
mi@0:     doctest.testmod()