
__version__ = "$Revision: 0.71 - 28 Oct 2010$"
__name__ = "pca"
__all__ = ["pca"]
__doc__ = """This module implements the Principal Component Analysis (PCA) algorithm.
           In the following functions input data M is shaped as in the following:
           M=[
                [var1, var2, ...var-mth] for sample 1, 
                [var1, var2, ...var-mth] for sample 2, 
                [...],
                [var1, var2, ...var-mth] for sample n, 
            ].
            Explanations of the resulting COEFF, SCORE, latent follow
            what reported in Matkab documentation:
            - COEFF = princomp(X) performs principal components analysis on the n-by-p data matrix X,
            and returns the principal component coefficients, also known as loadings. Rows of X
            correspond to observations, columns to variables. COEFF is a p-by-p matrix, each column
            containing coefficients for one principal component. The columns are in order of
            decreasing component variance;
            - SCORE = the principal component scores; that is, the representation of X in the principal
            component space. Rows of SCORE correspond to observations, columns to components;
            - latent = a vector containing the eigenvalues of the covariance matrix of X.
            """


try:
    import numpy
    from numpy.linalg import svd, eig
    __all__.append('numpy')
except ImportError:
    raise ImportError, "Impossible to import numpy/scipy module"


class pca:

    def __init__(self, M, num_type = numpy.float32): # cov_method = False):

        M0 = numpy.asarray(M, dtype = num_type)
        X = M0 - M0.mean(axis = 0)

        #################################################################################
        # The following notation is mutuated from that of Matlab:
        # - SCORE are the principal component scores, that is, the representation of the
        #   original matrix in the principal component space;
        # - COEFF are the principal component coefficients, also known as loadings, in a
        #   p-by-p matrix, where p is the number of variables equal to len(A[0]), and 
        #   containing the coefficients for each principal component;
        # - latent is a vector of eigenvalues of the covariance matrix of original data.
        #################################################################################
        
        #################################################################################
        # Here applying the Singular Value Decomposition (SVD) method to perform PCA
        # The Karhunen-Loeve transform is equivalent to finding the SVD of X.
        #################################################################################
        try:    self._svd_method(X)
        except: self._cov_method(X) ; print "PCA switched to covariance matrix based method.."
    
    def _svd_method(self, X):
        
        U, sigma, V = svd(X)
        self.COEFF  = V.T
        self.SCORE  = numpy.dot(X, self.COEFF)
        self.latent = self.SCORE.var(axis = 0, ddof = 1)
    
    def _cov_method(self, X):

        #X = numpy.copy(numpy.asarray(M, dtype = num_type)) ; X -= X.mean(axis = 0)
        #################################################################################
        # Covariance Method.
        # Here the eigenvalues, each repeated according to its multiplicity. 
        # The eigenvalues are not necessarily ordered, nor are they necessarily real for
        # real matrices
        # (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig)
        #################################################################################
        Cov = numpy.cov(X, rowvar = 0)
        latent, COEFF = numpy.linalg.eig(Cov)
        ixLatentSorted = list(zip(*sorted([(l, i) for i, l in enumerate(latent)])[::-1])[1])
        self.latent = latent[ixLatentSorted]
        self.COEFF = COEFF[:,ixLatentSorted] 
        self.SCORE = numpy.dot(X, self.COEFF)

    def __getitem__(self, key):
        
        if hasattr(self, 'SCORE'): return self.SCORE[key]
        return None


"""
if __name__ == 'pca':
   
    import random as rnd
    from pylab import *
    
    ##############################################################
    # Some example where PCA is not seleceting as first component 
    # the most useful to a subsequent clustering
    ##############################################################
    D = [[rnd.gauss(0, 1), rnd.gauss(0, 1)] for i in range(1000)]

    #X = [rnd.gauss(0, 1) for i in range(1000)]
    #Y = [0.5 * t + rnd.gauss(0,1) for t in X]
    #D = zip(X, Y)
    
    #D = [[rnd.gauss(3, 0.5), rnd.gauss(2, 0.5)] for j in range(100)] + \
    #    [[rnd.gauss(6, 0.5), rnd.gauss(2, 0.5)] for j in range(100)] + \
    #    [[rnd.gauss(3, 0.5), rnd.gauss(6, 1)] for j in range(100)] + \
    #    [[rnd.gauss(6, 0.5), rnd.gauss(6, 1)] for j in range(100)] + \
    #    [[rnd.gauss(4, 1), rnd.gauss(4, 1)] for j in range(100)]

    #figure(1)
    figure()
    axis('equal')
    N = array(D, dtype = float)
    plot(N[:, 0], N[:, 1], 'r.', markeredgecolor = 'r')

    P = pca(D)
"""
"""
    print 'PCA coeff.:', P.COEFF
    print 'PCA score:', P.SCORE
    print 'PCA % latent:',
    for p in P.latent: print '%f' % (p * 100 / sum(P.latent)),
    print '\n'

    #A = [list(i) for i in P.SCORE]
    #print len(A), len(A[0]), A[0]

    #figure(2)
    score = P.SCORE
    axis('equal')
    plot(score[:, 0], score[:, 1], 'b.', markeredgecolor = 'b')
    show()
    
    #close('all')
    
    ######################
    # PCA with 4-dim data
    ######################
    D = [[rnd.gauss(3, 0.5), rnd.gauss(2, 0.5), rnd.gauss(2, 0.5), rnd.gauss(2, 0.5)] for j in range(100)]
    D.extend([[rnd.gauss(6, 0.5), rnd.gauss(2, 0.5), rnd.gauss(1, 0.5), rnd.gauss(1, 0.5)] for j in range(100)])
    D.extend([[rnd.gauss(3, 0.5), rnd.gauss(6, 1), rnd.gauss(0, 0.5), rnd.gauss(6, 1)] for j in range(100)])
    D.extend([[rnd.gauss(6, 0.5), rnd.gauss(6, 1), rnd.gauss(6, 0.5), rnd.gauss(6, 0.5)] for j in range(100)])
    P = pca(D)
    print P.SCORE
    print P.COEFF
    print P.latent
"""
