# Code from Chapter 6 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008, 2014

# The Kernel PCA algorithm

import numpy as np
import pylab as pl

def kernelmatrix(data,kernel,param=np.array([3,2])):
    
    if kernel=='linear':
        return np.dot(data,transpose(data))
    elif kernel=='gaussian':
        K = np.zeros((np.shape(data)[0],np.shape(data)[0]))
        for i in range(np.shape(data)[0]):
            for j in range(i+1,np.shape(data)[0]):
                K[i,j] = np.sum((data[i,:]-data[j,:])**2)
                K[j,i] = K[i,j]
        return np.exp(-K**2/(2*param[0]**2))
    elif kernel=='polynomial':
        return (np.dot(data,np.transpose(data))+param[0])**param[1]
    
def kernelpca(data,kernel,redDim):
    
    nData = np.shape(data)[0]
    nDim = np.shape(data)[1]
    
    K = kernelmatrix(data,kernel)
    
    # Compute the transformed data
    D = np.sum(K,axis=0)/nData
    E = np.sum(D)/nData
    J = np.ones((nData,1))*D
    K = K - J - np.transpose(J) + E*np.ones((nData,nData))
    
    # Perform the dimensionality reduction
    evals,evecs = np.linalg.eig(K) 
    indices = np.argsort(evals)
    indices = indices[::-1]
    evecs = evecs[:,indices[:redDim]]
    evals = evals[indices[:redDim]]
    
    sqrtE = np.zeros((len(evals),len(evals)))
    for i in range(len(evals)):
        sqrtE[i,i] = np.sqrt(evals[i])
       
    #print shape(sqrtE), shape(data)
    newData = np.transpose(np.dot(sqrtE,np.transpose(evecs)))
    
    return newData

#data = array([[0.1,0.1],[0.2,0.2],[0.3,0.3],[0.35,0.3],[0.4,0.4],[0.6,0.4],[0.7,0.45],[0.75,0.4],[0.8,0.35]])
#newData = kernelpca(data,'gaussian',2)
#plot(data[:,0],data[:,1],'o',newData[:,0],newData[:,0],'.')
#show()