# -*- coding: utf-8 -*-
"""
Created on Fri Nov 11 16:58:06 2016

@author: ibackus
"""
import numpy as np

    
def reconstruct(x, v, nUse):
    """
    Reconstruct x using nUse principle components
    Parameters
    ----------
    x : array
    """
    vshort = v[:, 0:nUse]
    vvT = np.dot(vshort, vshort.T)
    return np.dot(x, vvT)

def PCA(x):
    """
    Perform PCA on x.  Also see reconstruct() to use the output to reconstruct
    x after projecting onto nUseprinciple components
    
    Parameters
    ----------
    x : array
        Feature array shape (nDataPts x nFeatures)
    
    Returns
    -------
    w : array
        Eigen-values for the principle components
    v : array
        Principle component array.  Components are stored along columns
    """
    # Get shape
    n, d = x.shape
    # Get principle components (columns of v)
    sigma = np.dot(x.T, x)/n
    # eigenvalues, eigenvectors
    w, v = np.linalg.eig(sigma)
    v = np.real(v)
    w = np.real(w)
    # make sure to sort!
    ind = (-w).argsort()
    w = w[ind]
    v = v[:, ind]
    print 'Sum of eigenvalues (direct): {0:.5g}'.format(w.sum())
    print 'Sum of eigenvalues (from xTx): {0:.5g}'.format(np.trace(sigma))
    return w, v


def proj(x, v, nUse=None):
    """
    Project x onto nUse principle components
    """
    if nUse is None:
        
        nUse = v.shape[1]
        
    return np.dot(x, v[:, 0:nUse])
    
def reconstructionError(w):
    """
    From the eigenvalues generated by PCA(), calculate the reconstruction error
    """
    a = np.zeros(w.shape)
    a[1:] = np.cumsum(w)[0:-1]
    a /= a[-1]
    error = 1 - a
    
    return error