
# 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 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 < 82:
        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

def my_mvpa(Y,num2):

    #Using PYMVPA
    PCA_data = np.array(Y)
    PCA_label_2 = ['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=num2)
    terr = TransferError(clf)
    ds1 = Dataset(samples=PCA_data,labels=PCA_label_2)
    cvterr = CrossValidatedTransferError(terr,NFoldSplitter(cvtype=1),enable_states=['confusion'])
    error = cvterr(ds1)
    return (1-error)*100

def result(eigvec_total,eigval_total,mean_data_total,B,C,num_PC):
    
# Reduced Eigen-Vector Matrix according to highest Eigenvalues..(Considering First 20 based on above figure)

    W = eigvec_total[:,0:num_PC]
    m_W, n_W = np.shape(W)
        
#Projected Data:
    Y = (W.T)*B
    m_Y, n_Y = np.shape(Y.T)
    return Y.T
  

if __name__ == '__main__':

   Fmat = np.row_stack([Fmat_original[0:41,:], Fmat_original[82:123,:]])

# 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)
    num_PC=1
    while num_PC <=20:
         Proj = np.zeros((90,num_PC))
         Proj = result(eigvec_total,eigval_total,mean_data_total,B,C,num_PC)
         # PYMVPA:
         num=0
         cv_acc = np.zeros(21)
         while num <=20:
             cv_acc[num] = my_mvpa(Proj,num)
             num = num+1
         plot(np.arange(21),cv_acc,'-s')
         grid('True')
         hold('True')
         num_PC = num_PC+1
    legend(('1-PC', '2-PCs', '3-PCs', '4-PCs', '5-PCs', '6-PCs', '7-PCs', '8-PCs', '9-PCs', '10-PCs', '11-PC', '12-PCs', '13-PCs', '14-PCs', '15-PCs', '16-PCs', '17-PCs', '18-PCs', '19-PCs', '20-PCs'))
    ylabel('Cross-Validation Accuracy')
    xlabel('k in k-NN Classifier')
    show()



