
# Principal Component Analysis Code : 

from numpy import mean,cov,double,cumsum,dot,linalg,array,rank,size,flipud
from pylab import *
import numpy as np
import matplotlib.pyplot as pp
#from enthought.mayavi import mlab

import scipy.ndimage as ni

import roslib; roslib.load_manifest('sandbox_tapo_darpa_m3')
import rospy
#import hrl_lib.mayavi2_util as mu
import hrl_lib.viz as hv
import hrl_lib.util as ut
import hrl_lib.matplotlib_util as mpu
import pickle

from mvpa.clfs.knn import kNN
from mvpa.datasets import Dataset
from mvpa.clfs.transerror import TransferError
from mvpa.misc.data_generators import normalFeatureDataset
from mvpa.algorithms.cvtranserror import CrossValidatedTransferError
from mvpa.datasets.splitters import NFoldSplitter

import sys
sys.path.insert(0, '/home/tapo/svn/robot1_data/usr/tapo/data_code/BMED_8813_HAP/Data')

from data_method_I import Fmat_original

def pca(X):

    #get dimensions
    num_data,dim = X.shape

    #center data
    mean_X = X.mean(axis=1)
    M = (X-mean_X) # subtract the mean (along columns)
    Mcov = cov(M)

    ###### Sanity Check ######
    i=0
    n=0
    while i < 123:
        j=0
        while j < 90:
            if X[i,j] != X[i,j]:
                print X[i,j]
                print i,j
                n=n+1
            j = j+1
        i=i+1

    print n
    ##########################

    print 'PCA - COV-Method used'
    val,vec = linalg.eig(Mcov)
  
    #return the projection matrix, the variance and the mean
    return vec,val,mean_X, M, Mcov

if __name__ == '__main__':

    Fmat = Fmat_original

# Checking the Data-Matrix

    m_tot, n_tot = np.shape(Fmat)
    
    print 'Total_Matrix_Shape:',m_tot,n_tot
    
    eigvec_total, eigval_total, mean_data_total, B, C = pca(Fmat)
    
    #print eigvec_total
    #print eigval_total
    #print mean_data_total
    m_eigval_total, n_eigval_total = np.shape(np.matrix(eigval_total))
    m_eigvec_total, n_eigvec_total = np.shape(eigvec_total)
    m_mean_data_total, n_mean_data_total = np.shape(np.matrix(mean_data_total))
    print 'Eigenvalue Shape:',m_eigval_total, n_eigval_total
    print 'Eigenvector Shape:',m_eigvec_total, n_eigvec_total
    print 'Mean-Data Shape:',m_mean_data_total, n_mean_data_total

#Recall that the cumulative sum of the eigenvalues shows the level of variance accounted by each of the corresponding eigenvectors. On the x axis there is the number of eigenvalues used.     

    perc_total = cumsum(eigval_total)/sum(eigval_total)
    
# Reduced Eigen-Vector Matrix according to highest Eigenvalues..(Considering First 20 based on above figure)

    W = eigvec_total[:,0:2]
    m_W, n_W = np.shape(W)
    print 'Reduced Dimension Eigenvector Shape:',m_W, n_W
    
#Projected Data:
    Y = (W.T)*B 
    m_Y, n_Y = np.shape(Y.T)
    print 'Transposed Projected Data Shape:', m_Y, n_Y

#Using PYMVPA
    PCA_data = np.array(Y.T)
    PCA_label_1 = ['Edge-1']*30 + ['Surface']*30 + ['Edge-2']*30
    PCA_chunk_1 = ['Can-Edge-1']*5 + ['Book-Edge-1']*5 + ['Brown-Cardboard-Box-Edge-1']*5 + ['Cinder-Block-Edge-1']*5 + ['Tin-Box-Edge-1']*5 + ['White-Cardboard-Box-Edge-1']*5 + ['Can-Surface']*5 + ['Book-Surface']*5 + ['Brown-Cardboard-Box-Surface']*5 + ['Cinder-Block-Surface']*5 + ['Tin-Box-Surface']*5 + ['White-Cardboard-Box-Surface']*5 + ['Can-Edge-2']*5 + ['Book-Edge-2']*5 + ['Brown-Cardboard-Box-Edge-2']*5 + ['Cinder-Block-Edge-2']*5 + ['Tin-Box-Edge-2']*5 + ['White-Cardboard-Box-Edge-2']*5
    clf = kNN(k=10)
    terr = TransferError(clf)
    ds1 = Dataset(samples=PCA_data,labels=PCA_label_1,chunks=PCA_chunk_1)
    print ds1.samples.shape
    cvterr = CrossValidatedTransferError(terr,NFoldSplitter(cvtype=1),enable_states=['confusion'])
    error = cvterr(ds1)
    print error
    print cvterr.confusion.asstring(description=False)
    figure(1)
    cvterr.confusion.plot(numbers='True')
    
# Variances
    figure(2)
    title('Variances of PCs')
    stem(range(len(perc_total)),perc_total,'--b')
    axis([-0.3,30.3,0,1.2])
    grid('True')
    show()


