# Compute cost function for feedforward net

import scipy.io as sio
import numpy as np
import math
import matplotlib.pyplot as plt

def sigmoid_func(z):
    return float(1) / (1 + math.e**(-z))


def nnCostFunction(Theta1, Theta2,  input_layer_size, hidden_layer_size, numLabels, X, y, lval):

    J = 0

    nhiddennodes,nfeat  = Theta1.shape
    nclasses, nhidden2 = Theta2.shape
    

    m, nxfeat = X.shape
    Aappend = np.ones((m, nxfeat + 1))
    Aappend[:, 1:] = X
    print 'A', Aappend.shape
    z2 = Theta1.dot(Aappend.T)
    a2 = sigmoid_func(z2)

    A2append = np.ones((m, nhidden2))
    A2append[:, 1:] = a2.T
    z3 = Theta2.dot(A2append.T)

    a3 = sigmoid_func(z3)

    # Reshape y to vectors of incidator variables
    # Classes 1-10: zero=10
    yshape = np.eye(nclasses, nclasses)
    ymat = np.zeros((m,nclasses))
    for i in range(0,m):

        curry = yshape[:,y[i]-1].T
        ymat[i,:] = curry





    jtest = 0.0
    a3y = np.zeros((m,nclasses))
    # Test in a loop
    for i in range(0,m):


        yvec = ymat[i,:]
        currhvec = a3[:,i].T
        oneminusy = 1-yvec
        oneminush = 1-currhvec
        a3y[i,:] = (a3[:,i]-yvec).T
        logh = np.log(currhvec)
        log1minus = np.log(oneminush)
        term1 = np.multiply(-yvec,logh)
        term2 = np.multiply(-oneminusy,log1minus)
        currvalvec = np.sum(term1)+np.sum(term2)
        jtest = jtest +currvalvec



    jtest = jtest/m

    print 'Cost function without regularization', jtest


    J = jtest




    # Then include regularization
    regterm = 0
    th1J, th1K = Theta1.shape
    th2J, th2K = Theta2.shape

    t1sum = 0
    t2sum = 0
    #Straightforward implementation
    for j in range(0,th1J):
        for k in range(1,th1K):
            t1sum = t1sum+np.square(Theta1[j,k])

    for j in range(0,th2J):
        for k in range(1,th2K):
            t2sum = t2sum + np.square(Theta2[j,k])

    regterm = (lval*t1sum+t2sum)/(2*m)


    J = jtest+regterm
    print 'Cost function with regularization ', J




    return J



# Set up the number of nodes
input_layer_size = 400
hidden_layer_size = 25
num_labels = 10

# Load and visualize the data
#THis contains the MNIST images of handwritten digits
datastruct = sio.loadmat('ex4data1.mat')

X = datastruct['X']
y=datastruct['y']
nsamp,nfeat = X.shape
#Change label10 to 0 for '0's


#y[y==10]=0

print 'yy', y.min(), y.max()
tt=1
#print np.mean(X[14,:])
#print X.shape
#y = datastruct.y

#Visualize 100 random images
#sel = np.random.permutation(nsamp)
#sel = sel[:100]
#displayData(X[sel,:])


# Load pre-trained net weights
wstruct = sio.loadmat('week6weights.mat')
#wstruct = sio.loadmat('ex4weights.mat')
Theta1 = wstruct['Theta1']
Theta2 = wstruct['Theta2']
#print Theta1.shape, Theta2.shape
lrat = 1

J =  nnCostFunction(Theta1, Theta2, input_layer_size, hidden_layer_size, num_labels, X, y,lrat)

print 'J', J