import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

######### Load the data ##########

infile = open('faces.csv','r')
img_data = infile.read().strip().split('\n')
img = [map(int,a.strip().split(',')) for a in img_data]
pixels = []
for p in img:
    pixels += p
faces = np.reshape(pixels,(400,4096))


######### Global Variable ##########

image_count = 0

######### Function that normalizes a vector x (i.e. |x|=1 ) #########

# > numpy.linalg.norm(x, ord=None, axis=None, keepdims=False) 
#   This function is able to return one of eight different matrix norms, 
#   or one of an infinite number of vector norms (described below), 
#   depending on the value of the ord parameter.
#   Note: in the given functionm, U should be a vector, not a array. 
#         You can write your own normalize function for normalizing 
#         the colomns of an array.

def normalize(U):
	return U / LA.norm(U) 

######### Display first face #########

# Useful functions:
# > numpy.reshape(a, newshape, order='C')
#   Gives a new shape to an array without changing its data.
# > matplotlib.pyplot.figure()
# 	Creates a new figure.
# > matplotlib.pyplot.title()
#	Set a title of the current axes.
# > matplotlib.pyplot.imshow()
#	Display an image on the axes.
#	Note: You need a matplotlib.pyplot.show() at the end to display all the figures.

first_face = np.reshape(faces[0],(64,64),order='F')
image_count+=1
plt.figure(image_count)
plt.title('First_face')
plt.imshow(first_face,cmap=plt.cm.gray)
#plt.show()


########## display a random face ###########

# Useful functions:
# > numpy.random.choice(a, size=None, replace=True, p=None)
#   Generates a random sample from a given 1-D array
# > numpy.ndarray.shape()
#   Tuple of array dimensions.
#   Note: There are two ways to order the elements in an array: 
#         column-major order and row-major order. In np.reshape(), 
#         you can switch the order by order='C' for row-major(default), 
#         or by order='F' for column-major. 


#### Your Code Here ####
rand_idx = np.random.choice(400)
rand_face = np.reshape(faces[rand_idx],(64,64),order = 'F')
image_count += 1
plt.figure(image_count)
plt.title('Random_face')
plt.imshow(rand_face,cmap = plt.cm.gray)
#plt.show()


########## compute and display the mean face ###########

# Useful functions:
# > numpy.mean(a, axis='None', ...)
#   Compute the arithmetic mean along the specified axis.
#   Returns the average of the array elements. The average is taken over 
#   the flattened array by default, otherwise over the specified axis. 
#   float64 intermediate and return values are used for integer inputs.

#### Your Code Here ####
mean = np.mean(faces,axis = 0)
mean_face = np.reshape(mean, (64,64), order='F')
image_count+=1
plt.figure(image_count)
plt.title('Mean_face')
plt.imshow(mean_face,cmap=plt.cm.gray)

######### substract the mean from the face images and get the centralized data matrix A ###########

# Useful functions:
# > numpy.repeat(a, repeats, axis=None)
#   Repeat elements of an array.

#### Your Code Here ####
faces_ctr = np.repeat(faces,1,axis = 0)
faces_ctr = faces_ctr - mean

######### calculate the eigenvalues and eigenvectors of the covariance matrix #####################

# Useful functions:
# > numpy.matrix()
#   Returns a matrix from an array-like object, or from a string of data. 
#   A matrix is a specialized 2-D array that retains its 2-D nature through operations. 
#   It has certain special operators, such as * (matrix multiplication) and ** (matrix power).

# > numpy.matrix.transpose(*axes)
#   Returns a view of the array with axes transposed.

# > numpy.linalg.eig(a)[source]
#   Compute the eigenvalues and right eigenvectors of a square array.
#   The eigenvalues, each repeated according to its multiplicity. 
#   The eigenvalues are not necessarily ordered. 

#### Your Code Here ####
faces_matrix = np.matrix(faces_ctr)
temp = np.transpose(faces_matrix)
'''
cov_mat = temp * faces_matrix
cov_mat = cov_mat / 400
'''
L = faces_matrix * temp

lamda, v_temp = LA.eig(L)
eigen_vectors = []
#print(lamda)
#print(v_temp)

for i in range(len(v_temp)):
    temp = (np.transpose(faces_matrix)) * (np.transpose(v_temp[i]))
    eigen_vectors.append(normalize(temp))
#print(eigen_vectors)

########## Display the first 16 principal components ##################

#### Your Code Here ####
eigen_index = np.argsort(lamda)[::-1]
for i in range(16):
    temp_face = np.reshape(eigen_vectors[ eigen_index[i] ],(64,64),order='F')
    image_count+=1
    plt.figure(image_count)
    plt.title(str(i+1)+ " principle compoent")
    plt.imshow(temp_face*255,cmap=plt.cm.gray)


########## Reconstruct the first face using the first two PCs #########

#### Your Code Here ####

z = np.array([eigen_vectors[0].T,eigen_vectors[1].T])
z = np.matrix(z)
q = z * (np.matrix(faces[0]).T)
result = (z.T * q).T
#print(result.shape) 
result += mean
result = np.reshape(result, (64,64),'F')
image_count+=1
plt.figure(image_count)
plt.title("the first reconstruct")
plt.imshow(result, cmap=plt.cm.gray)

########## Reconstruct random face using the first 5, 10, 25, 50, 100, 200, 300, 399  PCs ###########

#### Your Code Here ####
reconstruct_face = np.matrix(faces[99])
#print(reconstruct_face.shape)
#print(eigen_index)

parameters = [5,10,25,50,100,200,300,399]
for i in range(len(parameters)):
	z = []
	for j in range(parameters[i]):
		index = (np.argwhere(eigen_index == j))[0][0]
		z.append(eigen_vectors[index])	
	z = np.array(z)
	z = np.matrix(z)
	#print(z.shape)
	#print(z)
	q = z * reconstruct_face.T
	result = (z.T * q).T
	result += mean
	result = np.reshape(result, (64,64), 'F')
	image_count += 1
	plt.figure(image_count)
	plt.title("reconstruct face of index 99 using " + str(parameters[i]) + " PCs")
	plt.imshow(result, cmap=plt.cm.gray)
	
#plt.show()
			

######### Plot proportion of variance of all the PCs ###############

# Useful functions:
# > matplotlib.pyplot.plot(*args, **kwargs)
#   Plot lines and/or markers to the Axes. 
# > matplotlib.pyplot.show(*args, **kw)
#   Display a figure. 
#   When running in ipython with its pylab mode, 
#   display all figures and return to the ipython prompt.

#### Your Code Here ####
total_cov = sum(lamda)
y = lamda / total_cov
x = [i for i in range(400)]
image_count+=1
plt.figure(image_count)
plt.title("plt")
plt.plot(x,y)

plt.show()


