ex4

import os
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from sklearn.preprocessing import OneHotEncoder

# Read the dataset
path1 = os.path.dirname(os.path.realpath(__file__)) + '\..\data\ex4data1'
data = sio.loadmat(path1)
# Load weights
path2 = os.path.dirname(os.path.realpath(__file__)) + '\..\data\ex4weights'
weights = sio.loadmat(path2)
theta1_pt = weights['Theta1']
theta2_pt = weights['Theta2']

# Visualizing the data
idx = np.random.randint(0, data['X'].shape[0], size=100)

def displayData(data, m, n):
    # Create 10x10 sub plots
    f, ax = plt.subplots(m,n)
    for j in range(m):
        for k in range(n):
            ax[j,k].imshow(data['X'][idx[j*m+k]].reshape(20,20))

#displayData(data,10,10)

X = data['X']
y = data['y']

learning_rate = 1

#from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder(sparse=False)
y_onehot = encoder.fit_transform(y)
y_onehot.shape

# Sigmoid function
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def forward_propagate(X, theta1, theta2):
    m = X.shape[0]

    a1 = np.insert(X, 0, values=np.ones(m), axis=1)
    z2 = a1 * theta1.T
    a2 = np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1)
    z3 = a2 * theta2.T
    h = sigmoid(z3)

    return a1, z2, a2, z3, h

# cost function with pre-trained weights and no regularization
def costPreTrained(X, y, theta1, theta2):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)

    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)

    first = np.sum(np.multiply(-y, np.log(h)))
    second = np.sum(np.multiply(1-y, np.log(1-h)))

    J = 1/m * (first - second)

    return J

print(costPreTrained(X, y_onehot, theta1_pt, theta2_pt))

# Regularized cost function with pre-trained weights
def costRegPreTrained(X, y, theta1, theta2, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)

    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    first = np.sum(np.multiply(-y, np.log(h)))
    second = np.sum(np.multiply(1-y, np.log(1-h)))

    tmp1 = theta1[:,1:]
    tmp2 = theta2[:,1:]
    reg = learning_rate/(2*m) * (np.sum(np.square(tmp1)) + np.sum(np.square(tmp2)))

    J = 1/m * (first - second) + reg

    return J

print(costRegPreTrained(X, y_onehot, theta1_pt, theta2_pt, learning_rate))

# Sigmoid Gradient
def sigmoidGradient(z):
    return np.multiply(sigmoid(z), (1 - sigmoid(z)))

print(sigmoidGradient(0))

# Normally initialize
epsilon_init = 0.25
input_size = 400
hidden_size = 25
num_labels = 10
learning_rate = 1

params = np.random.random(size=((input_size+1)*hidden_size + (hidden_size+1)*num_labels)) * 2 * epsilon_init - epsilon_init

# Backpropagation
def backpropagate(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)

    theta1 = np.reshape(params[:(input_size+1)*hidden_size], (hidden_size, input_size+1))
    theta2 = np.reshape(params[(input_size+1)*hidden_size:], (num_labels, hidden_size+1))

    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)

    # cost
    J = costRegPreTrained(X, y, theta1, theta2, learning_rate)

    # initializations
    delta_1 = np.zeros(theta1.shape)  # shape (25, 401)
    delta_2 = np.zeros(theta2.shape)  # shape (10, 26)

    # perform backpropagation
    for i in range(m):
        a1i = a1[i,:]  # (1, 401)
        z2i = z2[i,:]  # (1, 25)
        a2i = a2[i,:]  # (1, 26)
        hi = h[i,:]  # (1, 10)
        yi = y[i,:]  # (1, 10)

        d3i = hi - yi  # (1, 10)

        # add bias unit
        z2i = np.insert(z2i, 0, values=np.ones(1))  # (1, 26)
        d2i = np.multiply((d3i * theta2), sigmoidGradient(z2i))  # (1, 26)

        delta_1 = delta_1 + d2i[:,1:].T * a1i # (25, 401)
        delta_2 = delta_2 + d3i.T * a2i # (10, 26)

    # add the gradient regularization term
    delta_1[:,1:] = delta_1[:,1:] + (theta1[:,1:] * learning_rate) / m
    delta_2[:,1:] = delta_2[:,1:] + (theta2[:,1:] * learning_rate) / m

    # unravel the gradient matrices into a single array
    grad = np.concatenate((np.ravel(delta_1), np.ravel(delta_2)))

    return J, grad

J, grad = backpropagate(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)

print(J)
print(grad.shape)

from scipy.optimize import minimize

# minimize the objective function
fmin = minimize(fun=backpropagate, x0=params, args=(input_size, hidden_size, num_labels, X, y_onehot, learning_rate), method='TNC', jac=True, options={'maxiter': 250})

#print(fmin)

def predict(h):
    h_argmax = np.argmax(h, axis=1)
    return h_argmax + 1

X = np.matrix(X)
theta1 = np.matrix(np.reshape(fmin.x[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
theta2 = np.matrix(np.reshape(fmin.x[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))

a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
y_pred = predict(h)

correct = [1 if a == b else 0 for (a, b) in zip(y_pred, y)]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print('accuracy = {0}%'.format(accuracy * 100))