import numpy as np
from scipy.io import loadmat
from scipy.optimize import minimize
from sklearn.svm import SVC
import pickle


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
    """

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

    n_feature = mat.get("train1").shape[1]
    n_sample = 0
    for i in range(10):
        n_sample = n_sample + mat.get("train" + str(i)).shape[0]
    n_validation = 1000
    n_train = n_sample - 10 * n_validation

    # Construct validation data
    validation_data = np.zeros((10 * n_validation, n_feature))
    for i in range(10):
        validation_data[i * n_validation:(i + 1) * n_validation, :] = mat.get("train" + str(i))[0:n_validation, :]

    # Construct validation label
    validation_label = np.ones((10 * n_validation, 1))
    for i in range(10):
        validation_label[i * n_validation:(i + 1) * n_validation, :] = i * np.ones((n_validation, 1))

    # Construct training data and label
    train_data = np.zeros((n_train, n_feature))
    train_label = np.zeros((n_train, 1))
    temp = 0
    for i in range(10):
        size_i = mat.get("train" + str(i)).shape[0]
        train_data[temp:temp + size_i - n_validation, :] = mat.get("train" + str(i))[n_validation:size_i, :]
        train_label[temp:temp + size_i - n_validation, :] = i * np.ones((size_i - n_validation, 1))
        temp = temp + size_i - n_validation

    # Construct test data and label
    n_test = 0
    for i in range(10):
        n_test = n_test + mat.get("test" + str(i)).shape[0]
    test_data = np.zeros((n_test, n_feature))
    test_label = np.zeros((n_test, 1))
    temp = 0
    for i in range(10):
        size_i = mat.get("test" + str(i)).shape[0]
        test_data[temp:temp + size_i, :] = mat.get("test" + str(i))
        test_label[temp:temp + size_i, :] = i * np.ones((size_i, 1))
        temp = temp + size_i

    # Delete features which don't provide any useful information for classifiers
    sigma = np.std(train_data, axis=0)
    index = np.array([])
    for i in range(n_feature):
        if (sigma[i] > 0.001):
            index = np.append(index, [i])
    train_data = train_data[:, index.astype(int)]
    validation_data = validation_data[:, index.astype(int)]
    test_data = test_data[:, index.astype(int)]

    # Scale data to 0 and 1
    train_data /= 255.0
    validation_data /= 255.0
    test_data /= 255.0

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


def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

#def safe_ln(x, minval=0.0000000001): 
#    return np.log(x.clip(min=minval))

def blrObjFunction(initialWeights, *args):
    """
    blrObjFunction computes 2-class Logistic Regression error function and
    its gradient.

    Input:
        initialWeights: the weight vector (w_k) of size (D + 1) x 10 
        train_data: the data matrix of size N x D
        labeli: the label vector (y_k) of size N x 1 where each entry can be either 0 or 1 representing the label of corresponding feature vector

    Output: 
        error: the scalar value of error function of 2-class logistic regression
        error_grad: the vector of size (D+1) x 1 representing the gradient of
                    error function
    """
    train_data, labeli = args

    n_data = train_data.shape[0]
    n_features = train_data.shape[1]
    error = 0
    error_grad = np.zeros((n_features + 1, 1))
    minval = 0.0000000001
  
    initialWeights = initialWeights.reshape((n_features+1),1)
    
    column_of_ones = np.ones((n_data,1))
    train_data_bias_added = np.append(column_of_ones, train_data,axis=1)
    y = sigmoid(np.dot(train_data_bias_added,initialWeights))
    error_class_1 = np.dot(np.transpose(np.log(y.clip(min=minval))),labeli)
    error_class_2 = np.dot(np.transpose(np.log((1-y).clip(min=minval))),1-labeli)

    error = -1.0*(np.sum(error_class_1) + np.sum(error_class_2))/n_data
    
    y = y.reshape(y.shape[0],1)
    error_grad = np.dot(np.transpose(train_data_bias_added), y-labeli)/n_data
    error_grad = error_grad.flatten()
    return error, error_grad


def blrPredict(W, data):
    """
     blrObjFunction predicts the label of data given the data and parameter W 
     of Logistic Regression
     
     Input:
         W: the matrix of weight of size (D + 1) x 10. Each column is the weight 
         vector of a Logistic Regression classifier.
         X: the data matrix of size N x D
         
     Output: 
         label: vector of size N x 1 representing the predicted label of 
         corresponding feature vector given in data matrix

    """
    label = np.zeros((data.shape[0], 1))
    n_data = data.shape[0]
    
    column_of_ones = np.ones((n_data,1))
    data_bias_added = np.append(column_of_ones, data,axis=1)
    y = sigmoid(np.dot(data_bias_added,W))
    label = y.argmax(axis=1)
    label = label.reshape((data.shape[0],1))
    return label


def mlrObjFunction(params, *args):
    """
    mlrObjFunction computes multi-class Logistic Regression error function and
    its gradient.

    Input:
        initialWeights: the weight vector of size (D + 1) x 1
        train_data: the data matrix of size N x D
        labeli: the label vector of size N x 1 where each entry can be either 0 or 1
                representing the label of corresponding feature vector

    Output:
        error: the scalar value of error function of multi-class logistic regression
        error_grad: the vector of size (D+1) x 10 representing the gradient of
                    error function
    """
    train_data, labeli = args
    n_data = train_data.shape[0]
    n_feature = train_data.shape[1]
    error = 0
    minval = 0.0000000001
    error_grad = np.zeros((n_feature + 1, n_class))
    params = params.reshape((n_feature+1,10))
    column_of_ones = np.ones((n_data,1))
    train_data_bias_added = np.append(column_of_ones, train_data,axis=1)
    
    theta_k = np.exp(np.dot(train_data_bias_added,params))
    theta_sum_for_all_k = (np.sum(theta_k,axis=1).astype(float)).reshape((n_data,1))
    theta_nk = np.divide(theta_k,theta_sum_for_all_k)
    
    temp = np.multiply(labeli, (np.log(theta_nk.clip(min=minval))))
    error = -1.0* (np.sum(np.sum(temp))).astype(float)/n_data
    error_grad = (np.dot(np.transpose(theta_nk - labeli),train_data_bias_added))/n_data
  
    return error, np.transpose(error_grad).flatten()


def mlrPredict(W, data):
    """
     mlrObjFunction predicts the label of data given the data and parameter W
     of Logistic Regression

     Input:
         W: the matrix of weight of size (D + 1) x 10. Each column is the weight
         vector of a Logistic Regression classifier.
         X: the data matrix of size N x D

     Output:
         label: vector of size N x 1 representing the predicted label of
         corresponding feature vector given in data matrix

    """
    label = np.zeros((data.shape[0], 1))
    W = W.reshape((n_feature+1,10))
    n_data = data.shape[0]    
    column_of_ones = np.ones((n_data,1))
    data_bias_added = np.append(column_of_ones, data,axis=1)
    
    theta_k = np.exp(np.dot(data_bias_added, W))
    theta_sum_for_all_k = (np.sum(theta_k,axis=1).astype(float)).reshape((n_data,1))
    theta_nk = np.divide(theta_k,theta_sum_for_all_k)
    label = np.argmax(theta_nk,axis=1)
    
    label = label.reshape((n_data,1))
    return label


"""
Script for Logistic Regression
"""
train_data, train_label, validation_data, validation_label, test_data, test_label = preprocess()

# number of classes
n_class = 10

# number of training samples
n_train = train_data.shape[0]

# number of features
n_feature = train_data.shape[1]

Y = np.zeros((n_train, n_class))
for i in range(n_class):
    Y[:, i] = (train_label == i).astype(int).ravel()

# Logistic Regression with Gradient Descent
W = np.zeros((n_feature + 1, n_class))
initialWeights = np.zeros((n_feature + 1, 1))
opts = {'maxiter': 100}
for i in range(n_class):
    labeli = Y[:, i].reshape(n_train, 1)
    args = (train_data, labeli)
    nn_params = minimize(blrObjFunction, initialWeights, jac=True, args=args, method='CG', options=opts)
    W[:, i] = nn_params.x.reshape((n_feature + 1,))

# Find the accuracy on Training Dataset
predicted_label = blrPredict(W, train_data)
print('\n Training set Accuracy:' + str(100 * np.mean((predicted_label == train_label).astype(float))) + '%')

# Find the accuracy on Validation Dataset
predicted_label = blrPredict(W, validation_data)
print('\n Validation set Accuracy:' + str(100 * np.mean((predicted_label == validation_label).astype(float))) + '%')

# Find the accuracy on Testing Dataset
predicted_label = blrPredict(W, test_data)
print('\n Testing set Accuracy:' + str(100 * np.mean((predicted_label == test_label).astype(float))) + '%')

#with open('params.pickle', 'wb') as f1: 
#    pickle.dump(W, f1)
#    
"""
Script for Support Vector Machine
"""

print('\n\n--------------SVM-------------------\n\n')
##################
# YOUR CODE HERE #
##################
train_label = train_label.ravel()
validation_label = validation_label.ravel()
test_label = test_label.ravel()


print ('1. Linear kernel (all other parameters are kept default)\n')
clf = SVC(kernel='linear')
clf.fit(train_data, train_label)
predicted_label = clf.predict(train_data)
print('\n Training set Accuracy:' + str(100 * np.mean((predicted_label == train_label).astype(float))) + '%')
predicted_label = clf.predict(validation_data)
print('\n Validation set Accuracy:' + str(100*np.mean((predicted_label == validation_label).astype(float))) + '%')
predicted_label = clf.predict(test_data)
print('\n Testing set Accuracy:' + str(100*np.mean((predicted_label == test_label).astype(float))) + '%')

print ('2. Radial basis function with value of gamma setting to 1. \n')
clf2 = SVC(gamma=1)
clf2.fit(train_data, train_label)
predicted_label2 = clf2.predict(train_data)
print('\n Training set Accuracy:' + str(100 * np.mean((predicted_label2 == train_label).astype(float))) + '%')
predicted_label2 = clf2.predict(validation_data)
print('\n Validation set Accuracy:' + str(100*np.mean((predicted_label2 == validation_label).astype(float))) + '%')
predicted_label2 = clf2.predict(test_data)
print('\n Testing set Accuracy:' + str(100*np.mean((predicted_label2 == test_label).astype(float))) + '%')


print ('3. Radial basis function with value of gamma as default. \n')
clf3 = SVC()
clf3.fit(train_data, train_label)
predicted_label3 = clf3.predict(train_data)
print('\n Training set Accuracy:' + str(100 * np.mean((predicted_label3 == train_label).astype(float))) + '%')
predicted_label3 = clf3.predict(validation_data)
print('\n Validation set Accuracy:' + str(100*np.mean((predicted_label3 == validation_label).astype(float))) + '%')
predicted_label3 = clf3.predict(test_data)
print('\n Testing set Accuracy:' + str(100*np.mean((predicted_label3 == test_label).astype(float))) + '%')


print ('4. RBF. Gamma=auto and C varies from 10 to 100\n')
for c in range(10, 100, 10):
    print ('\nValue of C is: '+str(c))
    clf4 = SVC(C=c)
    clf4.fit(train_data, train_label)
    print('\nDone predicting. Calculating score now')
    predicted_label4 = clf4.predict(train_data)
    print('\n Training set Accuracy:' + str(100 * np.mean((predicted_label4 == train_label).astype(float))) + '%')
    predicted_label4 = clf4.predict(validation_data)
    print('\n Validation set Accuracy:' + str(100*np.mean((predicted_label4 == validation_label).astype(float))) + '%')
    predicted_label4 = clf4.predict(test_data)
    print('\n Testing set Accuracy:' + str(100*np.mean((predicted_label4 == test_label).astype(float))) + '%')

print('\n\n--------------END of SVM-------------------\n\n')


"""
Script for Extra Credit Part
"""
# FOR EXTRA CREDIT ONLY
W_b = np.zeros((n_feature + 1, n_class))
initialWeights_b = np.zeros((n_feature + 1, n_class))
opts_b = {'maxiter': 100}

args_b = (train_data, Y)
nn_params = minimize(mlrObjFunction, initialWeights_b, jac=True, args=args_b, method='CG', options=opts_b)
W_b = nn_params.x.reshape((n_feature + 1, n_class))

# Find the accuracy on Training Dataset
predicted_label_b = mlrPredict(W_b, train_data)
print('\n Training set Accuracy:' + str(100 * np.mean((predicted_label_b == train_label).astype(float))) + '%')

# Find the accuracy on Validation Dataset
predicted_label_b = mlrPredict(W_b, validation_data)
print('\n Validation set Accuracy:' + str(100 * np.mean((predicted_label_b == validation_label).astype(float))) + '%')

# Find the accuracy on Testing Dataset
predicted_label_b = mlrPredict(W_b, test_data)
print('\n Testing set Accuracy:' + str(100 * np.mean((predicted_label_b == test_label).astype(float))) + '%')
##with open('params_bonus.pickle', 'wb') as f2:
##    pickle.dump(W_b, f2)
