#CS428
#Assignment 2
#Student ID 15315901

"""
Import the images (as vectors) into Python and perform a principal component analysis.
Let P(n) denote the vector space generated by those eigenvectors corresponding to the n
largest eigenvalues. For n = 10,50,100 and 300 determine how much of the variability of
the database is captured by projecting onto P(n) ?

2. Take an image of yourself and store it in the same format as the AT&T images. Display,
as an image (rather than a vector), the projection of your original image onto P(n) for
n = 10,50,100 and 300.

3. Take an image of a friend and determine the distance between the projections of your own
image and your friend’s image onto P(300). Specify which metric you are using to compute
this distance



"""
# Relevant package import
import os, sys
import numpy as np
import PIL.Image as Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm



sys.path.append("..")


class AbstractDistance(object):
    def __init__(self,name):
        self._name = name

    def __call__(self, p, q):
        raise NotImplementedError("Every AbstractDistance must implement the __call__ method")

    @property
    def name(self):
        return self._name

    def __repr__(self):
        return self._name

class EuclideanDistance(AbstractDistance):

    def __init__(self):
        AbstractDistance.__init__(self, "EuclideanDistance")

    def __call__(self, p, q):
        p = np.asarray(p).flatten()
        q = np.asarray(q).flatten()
        return np.sqrt(np.sum(np.power((p-q),2)))

class BaseModel(object):
    def __init__(self, X=None, y=None, dist_metric=EuclideanDistance(), num_components=0):
        self.dist_metric = dist_metric
        self.num_components = 0
        self.projections=[]
        self.W = []
        self.mu = []
        if (X is not None) and (y is not None):
            self.compute(X, y)

    def compute(self, X, y):
        raise NotImplementedError("Every BaseModel must implement the compute method.")

    def predict(self, X):
        minDist = np.finfo('float').max
        minClass = -1
        Q = project(self.W, X.reshape(1,-1), self.mu)
        for i in range(len(self.projections)):
            dist = self.dist_metric(self.projections[i], Q)
            if dist < minDist:
                minDist = dist
                minClass = self.y[i]
        return minClass

class EigenfacesModel (BaseModel):
    def __init__(self, X=None, y=None, dist_metric=EuclideanDistance(), num_components=0):
        super(EigenfacesModel, self).__init__(X=X, y=y, dist_metric=dist_metric, num_components=num_components)

    def compute(self, X, y):
        [D, self.W, self.mu] = pca(asRowMatrix(X),y, self.num_components)
        self.y=y
        for xi in X:
            self.projections.append(project(self.W, xi.reshape(1,-1), self.mu))




# Function for reading images
def read_images(path, sz=None):
    c = 0
    X,y = [], []
    for dirname , dirnames , filenames in os.walk(path):
        for subdirname in dirnames:
            subject_path = os.path.join(dirname , subdirname)
            for filename in os.listdir(subject_path):
                try:
                    im = Image.open(os.path.join(subject_path , filename))
                    im = im.convert("L")
                    # resize to given size (if given)
                    if (sz is not None):
                        im = im.resize(sz, Image.ANTIALIAS)
                    X.append(np.asarray(im, dtype=np.uint8))
                    y.append(c)
                except IOError:
                    print("I/O error({0}): {1}".format(errno , strerror))
                except:
                    print("Unexpected error:", sys.exc_info()[0])
                    raise
            c = c+1

    return [X,y]

# Function to transform a given matrix into a row or column fom
def asRowMatrix(X):
    if len(X) == 0:
        return np.array([])
    mat = np.empty((0, X[0].size), dtype=X[0].dtype)
    for row in X:
        mat = np.vstack((mat, np.asarray(row).reshape(1,-1)))
    return mat

def asColumnMatrix(X):
    if len(X) == 0:
        return np.array([])
    mat = np.empty((0, X[0].size), dtype=X[0].dtype)
    for col in X:
        mat = np.hstack((mat, np.asarray(col).reshape(-1,1)))
    return mat

# Creating PCA funciton
def pca(X, y, num_components = 0):
    [n,d] = X.shape
    if (num_components<=0) or (num_components>n):
        num_components = n
    mu = X.mean(axis=0)
    X = X - mu
    if n>d:
        C = np.dot(X.T, X)
        [eigenvalues, eigenvectors] = np.linalg.eigh(C)
    else:
        C = np.dot(X,X.T)
        [eigenvalues, eigenvectors] = np.linalg.eigh(C)
        eigenvectors = np.dot(X.T, eigenvectors)
        for i in range(n):
            eigenvectors[:,i] = eigenvectors[:,i]/np.linalg.norm(eigenvectors[:,i])
    variance = np.linalg.svd(X.T, full_matrices=False)
    idx = np.argsort(-eigenvalues)
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:,idx]
    eigenvalues = eigenvalues[0:num_components].copy()
    eigenvectors = eigenvectors[:,0:num_components].copy()
    #print(num_components)
    #print(variance)
    return [eigenvalues, eigenvectors, mu]

# For rearranging projections
def project(W, X, mu=None):
    if mu is None:
        return np.dot(W, X)
    return np.dot(X - mu, W)


def reconstruct(W, Y, mu=None):
    if mu is None:
        return np.dot(Y, W.T)
    return np.dot(Y, W.T) + mu

def normalize(X, low, high, dtype=None):
    X = np.asarray(X)
    minX, maxX = np.min(X), np.max(X)
    X = X - float(minX)
    X = X / float((maxX - minX))
    X = X * (high-low)
    X = X + low
    if dtype is None:
        return np.asarray(X)
    return np.asarray(X, dtype=dtype)



def create_font(fontname='Tahoma', fontsize=10):
    return { 'fontname': fontname, 'fontsize':fontsize }

def subplot(title, images, rows, cols, sptitle="subplot", sptitles=[], colormap=cm.gray, ticks_visible=True, filename=None):
    fig = plt.figure()
    fig.text(.5, .95, title, horizontalalignment='center')
    for i in range(len(images)):
        ax0 = fig.add_subplot(rows,cols,(i+1))
        plt.setp(ax0.get_xticklabels(), visible=False)
        plt.setp(ax0.get_yticklabels(), visible=False)
        if len(sptitles) == len(images):
            plt.title("%s #%s" % (sptitle, str(sptitles[i])), create_font('Tahoma',10))
        else:
            plt.title("%s #%d" % (sptitle, (i+1)), create_font('Tahoma',10))
        plt.imshow(np.asarray(images[i]), cmap=colormap)

    if filename is None:
        plt.show()
    else:
        fig.savefig(filename)

def imsave(image, title="", filename=None):
    plt.figure()
    plt.imshow(np.asarray(image))
    plt.title(title, create_font('Tahoma',10))
    if filename is None:
        plt.show()
    else:
        fig.savefig(filename)





###
### Function end
###


[X,y] = read_images('orl_faces')
[D, W, mu] = pca(asRowMatrix(X), y)



model = EigenfacesModel(X[1:], y[1:])
#model = EigenfacesModel(X, y)


print("Expected = ", y[40], "/", "Predicted =", model.predict(X[40]))

"""
E = []
for i in range(min(len(X),16)):
    e = W[:,i].reshape(X[0].shape)
    E.append(normalize(e,0,255))

subplot(title="Eigenfaces AT&T Facedatabase", images=E, rows=4, cols=4, sptitle="Eigenface", colormap=cm.jet, filename="python_pca_eignefaces_color_with_personal_and_frank.pdf")

"""
steps = [i for i in range(10, min(len(X), 320), 20)]
E = []
for i in range(min(len(steps),16)):
    numEvs = steps[i]
    P = project(W[:, 0:numEvs], X[0].reshape(1,-1), mu)
    R = reconstruct(W[:,0:numEvs], P, mu)

    R = R.reshape(X[0].shape)
    E.append(normalize(R,0,255))

subplot(title="Reconstruction AT&T Facedatabase", images=E, rows=8, cols=8, sptitle="Eigenvectors", sptitles=steps, colormap=cm.gray, filename="python_pca_eignefaces_with_personal_and_frank.pdf")
