import numpy as np
from scipy.optimize import minimize
from scipy.io import loadmat
from math import sqrt
import time
import pickle

t = time.time()

def initializeWeights(n_in,n_out):
    
    """
    # initializeWeights return the random weights for Neural Network given the
    # number of node in the input layer and output layer

    # Input:
    # n_in: number of nodes of the input layer
    # n_out: number of nodes of the output layer

    # Output: 
    # W: matrix of random initial weights with size (n_out x (n_in + 1))"""

    epsilon = sqrt(6) / sqrt(n_in + n_out + 1);
    W = (np.random.rand(n_out, n_in + 1)*2* epsilon) - epsilon;

    return W
 

def sigmoid(z):

    
    """# Notice that z can be a scalar, a vector or a matrix
    # return the sigmoid of input z"""

    return  1.0/(1.0+np.exp(-z))


def preprocess():

    """ Input:
     Although this function doesn't have any input, you are required to load
     the MNIST data set from file 'mnist_all.mat'.

     Output:
     train_data: matrix of training set. Each row of train_data contains 
       feature vector of a image
     train_label: vector of label corresponding to each image in the training
       set
     validation_data: matrix of training set. Each row of validation_data 
       contains feature vector of a image
     validation_label: vector of label corresponding to each image in the 
       training set
     test_data: matrix of training set. Each row of test_data contains 
       feature vector of a image
     test_label: vector of label corresponding to each image in the testing
       set

     Some suggestions for preprocessing step:
     - divide the original data set to training, validation and testing set
           with corresponding labels
     - convert original data set from integer to double by using double()
           function
     - normalize the data to [0, 1]
     - feature selection"""

    mat = loadmat('./mnist_all.mat') #loads the MAT object as a Dictionary

    #Pick a reasonable size for validation data
    #Your code here

    train_data = np.array([])
    train_label = np.array([])
    validation_data = np.array([])
    validation_label = np.array([])
    test_data = np.array([])
    test_label = np.array([])
    sizen = mat.get("train0").shape[1]

    totalsize = 0
    for i in range(10):
	    m = mat.get('train'+str(i))
	    totalsize = totalsize + m.shape[0]
    #print totalsize

    testsize = 0
    for i in range(10):
	    m = mat.get('test'+str(i))
	    testsize = testsize + m.shape[0]
    #print testsize

    valn = 1000
    trainsize = totalsize - 10*valn

    #print trainsize

    valsize = 10*valn

    #print valsize

    #initialization part

    test_data = np.zeros((0, sizen))
    train_data = np.zeros((0, sizen))
    validation_data = np.zeros((0, sizen))

    test_label = np.zeros((testsize, 1))
    train_label = np.zeros((trainsize, 1))
    validation_label = np.zeros((valsize, 1))

    #label initialization
    validation_label = np.ones((valsize,1))
    for i in range(10):
	    validation_label[(i*valn):((i+1)*valn),:] = i*np.ones((valn,1))

    temp1 = 0
    temp2 = 0
    for i in range(10):
	    sizei = mat.get("train"+str(i)).shape[0]
	    sizej = mat.get("test"+str(i)).shape[0]
	
	    train_label[temp1:temp1+sizei-valn,:] = i*np.ones((sizei-valn,1))
	    temp1 = temp1+sizei-valn
	
	    test_label[temp2:temp2+sizej,:] = i*np.ones((sizej,1))
	    temp2 = temp2 + sizej
	
    #test data part
    for i in range(10):
            A = mat.get('test'+str(i))
            test_data = np.concatenate([test_data, A], 0) 

    # train and val data part
    for i in range(10):
	    B = mat.get('train'+str(i))
	    b = range(B.shape[0])
	    bperm = np.random.permutation(b)
	    B1 = B[bperm[0:1000],:]
	    B2 = B[bperm[1000:],:]
	    validation_data = np.concatenate([validation_data, B1], 0)
	    train_data = np.concatenate([train_data, B2], 0)

    #print str(train_data.shape) + '-' + str(validation_data.shape) + '-' + str(test_data.shape)
    #print str(train_label.shape) + '-' + str(validation_label.shape) + '-' + str(test_label.shape)

    featarr=np.all(train_data==train_data[0,:],axis=0)
    indiarr=[]
    for i in range(featarr.size):
    	    if featarr[i]==True:
		    indiarr.append(i)

    train_data=np.delete(train_data,indiarr,1)
    #print (train_data.shape)
    validation_data=np.delete(validation_data,indiarr,1)
    #print (test_data.shape)
    test_data=np.delete(test_data,indiarr,1)
    #print (test_data.shape)

    train_data = np.double(train_data)/255.0
    validation_data =np.double(validation_data)/255.0
    test_data = np.double(test_data)/255.0
    #print train_label.shape

    return train_data, train_label, validation_data, validation_label, test_data, test_label


def nnObjFunction(params, *args):

    """% nnObjFunction computes the value of objective function (negative log 
    %   likelihood error function with regularization) given the parameters 
    %   of Neural Networks, thetraining data, their corresponding training 
    %   labels and lambda - regularization hyper-parameter.

    % Input:
    % params: vector of weights of 2 matrices w1 (weights of connections from
    %     input layer to hidden layer) and w2 (weights of connections from
    %     hidden layer to output layer) where all of the weights are contained
    %     in a single vector.
    % n_input: number of node in input layer (not include the bias node)
    % n_hidden: number of node in hidden layer (not include the bias node)
    % n_class: number of node in output layer (number of classes in
    %     classification problem
    % training_data: matrix of training data. Each row of this matrix
    %     represents the feature vector of a particular image
    % training_label: the vector of truth label of training images. Each entry
    %     in the vector represents the truth label of its corresponding image.
    % lambda: regularization hyper-parameter. This value is used for fixing the
    %     overfitting problem.

    % Output: 
    % obj_val: a scalar value representing value of error function
    % obj_grad: a SINGLE vector of gradient value of error function
    % NOTE: how to compute obj_grad
    % Use backpropagation algorithm to compute the gradient of error function
    % for each weights in weight matrices.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % reshape 'params' vector into 2 matrices of weight w1 and w2
    % w1: matrix of weights of connections from input layer to hidden layers.
    %     w1(i, j) represents the weight of connection from unit j in input 
    %     layer to unit i in hidden layer.
    % w2: matrix of weights of connections from hidden layer to output layers.
    %     w2(i, j) represents the weight of connection from unit j in hidden 
    %     layer to unit i in output layer."""

    n_input, n_hidden, n_class, training_data, training_label, lambdaval = args

    w1 = params[0:n_hidden * (n_input + 1)].reshape( (n_hidden, (n_input + 1)))
    w2 = params[(n_hidden * (n_input + 1)):].reshape((n_class, (n_hidden + 1)))
    obj_val = 0  

    #print w1.shape
    #print w2.shape
    #Your code here
    #print training_data.shape   
    ncount = training_data.shape[0] 
    #Converting training label to yl
    yl = np.zeros((ncount,10))
    for x in range (0,ncount):
        l = int(training_label[x])
        #print int(l)
        yl[x][l]=1

    #Input layer
    training_data_bias = np.ones((ncount,1))
    training_dataxy = np.hstack((training_data,training_data_bias))
    #print training_data1.shape
    aj = np.dot(training_dataxy, np.transpose(w1))
    #print aj.shape
    zj = sigmoid(aj)
    #print zj.shape

    #Hidden layer 
    hidden_bias = np.ones((ncount,1))
    trainingzj = np.hstack((zj, hidden_bias))
    bj = np.dot(trainingzj, np.transpose(w2))
    oj = sigmoid(bj)    

    #print 'ncount' + str(ncount) + 'trainbias' + str(training_data_bias.shape) + 'traindataxy' + str(training_dataxy.shape) + 'aj' + str(aj.shape) + 'zj' + str(zj.shape)
    #print 'hid' + str(hidden_bias.shape) + 'trainzj' + str(trainingzj.shape) + 'bj' + str(bj.shape) + 'oj' + str(oj.shape) 
    
    #delta l calculation Eq 9
    deltal = ((yl - oj) * (1 - oj) * oj)
    #print deltal.shape
	
	#J calculation Eq 8
    Jp_w2 = (-1)*(np.dot(np.transpose(deltal) , trainingzj))
    #print Jp_w2.shape

	#Error calculation - Eq 5
    J1 = np.square(yl - oj)
    #print yl.shape
    #print oj.shape
    #print J1.shape

    J1 = (np.sum(J1))/2.0
    J1 = np.sum(J1)/ncount
    #print J1

	#Regularization: Eq 15
    temp = ((np.sum(np.square(w1))) + (np.sum(np.square(w2))))
    #temp = ((np.sum(w1 * w1)) + (np.sum(w2 * w2)))

    r = (lambdaval * temp)/(2.0 * ncount)

    #print r
    obj_val = J1 + r 

	#Eq 12
    Z1 = ((-1) * (1-trainingzj) * (trainingzj))
    Z2 = (np.dot(deltal,w2))
  
    Jp_w1 = np.dot(np.transpose(Z1*Z2),training_dataxy) 


	#Eq 16
    deltaJ2 = ((Jp_w2) + (lambdaval * w2))/ ncount

	#Eq 17    
    deltaJ1 = ((np.delete(Jp_w1, (-1), axis = 0)) + (lambdaval * w1))/ ncount

    grad_w1 = deltaJ1
    grad_w2 = deltaJ2
    #print grad_w1.shape
    #print grad_w2.shape

    obj_grad = np.concatenate((grad_w1.flatten(), grad_w2.flatten()),0)

    print 'obj_val ' , obj_val
    #print 'obj_grad' , obj_grad.shape

    #Make sure you reshape the gradient matrices to a 1D array. for instance if your gradient matrices are grad_w1 and grad_w2
    #you would use code similar to the one below to create a flat array    #obj_grad = np.array([])

    return (obj_val,obj_grad)


def nnPredict(w1,w2,data):


    """% nnPredict predicts the label of data given the parameter w1, w2 of Neural
    % Network.

    % Input:
    % w1: matrix of weights of connections from input layer to hidden layers.
    %     w1(i, j) represents the weight of connection from unit i in input 
    %     layer to unit j in hidden layer.
    % w2: matrix of weights of connections from hidden layer to output layers.
    %     w2(i, j) represents the weight of connection from unit i in input 
    %     layer to unit j in hidden layer.
    % data: matrix of data. Each row of this matrix represents the feature 
    %       vector of a particular image

    % Output: 
    % label: a column vector of predicted labels""" 

    labels = np.array([])
    #Your code here
    #Input layer 
    labels = np.zeros((data.shape[0],1))
    ncount1 = data.shape[0]
    training_data_bias = np.ones((ncount1,1))
    data = np.hstack((data,training_data_bias))
    aj = np.dot(data, np.transpose(w1))
    zj = sigmoid(aj)

    #Hidden layer 
    hidden_bias = np.ones((ncount1,1))
    trainingzj = np.hstack((zj, hidden_bias))
    bj = np.dot(trainingzj, np.transpose(w2))
    oj = sigmoid(bj) 

    for i in range(oj.shape[0]):
       m = np.argmax(oj[i])
       labels[i] = m

    return labels

"""**************Neural Network Script Starts here********************************"""

train_data, train_label, validation_data,validation_label, test_data, test_label = preprocess();

#  Train Neural Network

# set the number of nodes in input unit (not including bias unit)
n_input = train_data.shape[1]; 

# set the number of nodes in hidden unit (not including bias unit)
n_hidden = 50;
			
# set the number of nodes in output unit
n_class = 10;				   

# initialize the weights into some random matrices
initial_w1 = initializeWeights(n_input, n_hidden);
initial_w2 = initializeWeights(n_hidden, n_class);

# unroll 2 weight matrices into single column vector
initialWeights = np.concatenate((initial_w1.flatten(), initial_w2.flatten()),0)

# set the regularization hyper-parameter
lambdaval = 0.1


args = (n_input, n_hidden, n_class, train_data, train_label, lambdaval)


#Train Neural Network using fmin_cg or minimize from scipy,optimize module. Check documentation for a working example

opts = {'maxiter' : 50}    # Preferred value.

nn_params = minimize(nnObjFunction, initialWeights, jac=True, args=args,method='CG', options=opts)


#In Case you want to use fmin_cg, you may have to split the nnObjectFunction to two functions nnObjFunctionVal
#and nnObjGradient. Check documentation for this function before you proceed.
#nn_params, cost = fmin_cg(nnObjFunctionVal, initialWeights, nnObjGradient,args = args, maxiter = 50)


#Reshape nnParams from 1D vector into w1 and w2 matrices

w1 = nn_params.x[0:n_hidden * (n_input + 1)].reshape( (n_hidden, (n_input + 1)))
w2 = nn_params.x[(n_hidden * (n_input + 1)):].reshape((n_class, (n_hidden + 1)))

#Test the computed parameters
predicted_label = nnPredict(w1,w2,train_data)

#find the accuracy on Training Dataset

print('\n Training set Accuracy:' + str(100*np.mean((predicted_label == train_label).astype(float))) + '%')

predicted_label = nnPredict(w1,w2,validation_data)

#find the accuracy on Validation Dataset

print('\n Validation set Accuracy:' + str(100*np.mean((predicted_label == validation_label).astype(float))) + '%')

predicted_label = nnPredict(w1,w2,test_data)

#find the accuracy on Validation Dataset

print('\n Test set Accuracy:' + str(100*np.mean((predicted_label == test_label).astype(float))) + '%')

t1 = time.time()
print('\n Time in seconds: ' + str(t1-t) + 's')

pickle_params = open('params.pickle','wb')
pickle.dump([n_hidden,w1,w2,lambdaval],pickle_params)
pickle_params.close()

