
# 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/Classification/Data/Single_Contact_kNN/Scaled')

from data_method_V 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 < 140:
            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 = ['Styrofoam-Fixed']*5 + ['Books-Fixed']*5 + ['Bucket-Fixed']*5 + ['Bowl-Fixed']*5 + ['Can-Fixed']*5 + ['Box-Fixed']*5 + ['Pipe-Fixed']*5 + ['Styrofoam-Movable']*5 + ['Container-Movable']*5 + ['Books-Movable']*5 + ['Cloth-Roll-Movable']*5 + ['Black-Rubber-Movable']*5 + ['Can-Movable']*5 + ['Box-Movable']*5 + ['Rug-Fixed']*5 + ['Bubble-Wrap-1-Fixed']*5 + ['Pillow-1-Fixed']*5 + ['Bubble-Wrap-2-Fixed']*5 + ['Sponge-Fixed']*5 + ['Foliage-Fixed']*5 + ['Pillow-2-Fixed']*5 + ['Rug-Movable']*5 + ['Bubble-Wrap-1-Movable']*5 + ['Pillow-1-Movable']*5 + ['Bubble-Wrap-2-Movable']*5 + ['Pillow-2-Movable']*5 + ['Cushion-Movable']*5 + ['Sponge-Movable']*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 
(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)
        
# Normalizes the data set with respect to its variance (Not an Integral part of PCA, but useful)
    length = len(eigval_total)
    s = np.matrix(np.zeros(length)).T
    i = 0
    while i < length:
        s[i] = sqrt(C[i,i])
        i = i+1
    Z = np.divide(B,s)
    m_Z, n_Z = np.shape(Z)
                
#Projected Data:
    Y = (W.T)*B # 'B' for my Laptop: otherwise 'Z' instead of 'B'
    m_Y, n_Y = np.shape(Y.T)
    return Y.T
  

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)
    num_PC=1
    while num_PC <=20:
         Proj = np.zeros((140,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()




result