import sys

import PIL.Image
import scipy.misc, scipy.optimize, scipy.io, scipy.special
from numpy import *

import pylab
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.mlab as mlab


def displayData( X, theta = None ):
	width = 20
	rows, cols = 10, 10
	out = zeros(( width * rows, width*cols ))
	rand_indices = random.permutation( 5000 )[0:rows * cols]

	counter = 0
	for y in range(0, rows):
		for x in range(0, cols):
			start_x = x * width
			start_y = y * width
			out[start_x:start_x+width, start_y:start_y+width] = X[rand_indices[counter]].reshape(width, width).T
			counter += 1

img = scipy.misc.toimage( out )
figure = pyplot.figure()
axes = figure.add_subplot(111)
axes.imshow( img )


if theta is not None:
result_matrix = []
X_biased = c_[ ones( shape(X)[0] ), X ]

for idx in rand_indices:
result = (argmax( theta.T.dot(X_biased[idx]) ) + 1) % 10
result_matrix.append( result )

result_matrix = array( result_matrix ).reshape( rows, cols ).transpose()
print result_matrix

pyplot.show( )


def sigmoid( z ):
return scipy.special.expit(z)
# return 1.0 / (1.0 + exp( -z ))

def computeCost( theta, X, y, lamda ):
m = shape( X )[0]
hypo = sigmoid( X.dot( theta ) )
term1 = log( hypo ).dot( -y )
term2 = log( 1.0 - hypo ).dot( 1 - y )
left_hand = (term1 - term2) / m
right_hand = theta.T.dot( theta ) * lamda / (2*m)
return left_hand + right_hand

def gradientCost( theta, X, y, lamda ):
m = shape( X )[0]
grad = X.T.dot( sigmoid( X.dot( theta ) ) - y ) / m
grad[1:] = grad[1:] + ( (theta[1:] * lamda ) / m )
return grad

def oneVsAll( X, y, num_classes, lamda ):
m,n = shape( X )
X = c_[ones((m, 1)), X]
all_theta = zeros((n+1, num_classes))

for k in range(0, num_classes):
theta = zeros(( n+1, 1 )).reshape(-1)
temp_y = ((y == (k+1)) + 0).reshape(-1)
result = scipy.optimize.fmin_cg( computeCost, fprime=gradientCost, x0=theta, \
args=(X, temp_y, lamda), maxiter=50, disp=False, full_output=True )
all_theta[:, k] = result[0]
print "%d Cost: %.5f" % (k+1, result[1])

# save( "all_theta.txt", all_theta )
return all_theta

def predictOneVsAll( theta, X, y ):
m,n = shape( X )
X = c_[ones((m, 1)), X]

correct = 0
for i in range(0, m ):
prediction = argmax(theta.T.dot( X[i] )) + 1
actual = y[i]
# print "prediction = %d actual = %d" % (prediction, actual)
if actual == prediction:
correct += 1
print "Accuracy: %.2f%%" % (correct * 100.0 / m )


def part1_1():
mat = scipy.io.loadmat( "ex3data1.mat" )
X, y = mat['X'], mat['y']

def part1_2():
mat = scipy.io.loadmat( "ex3data1.mat" )
X, y = mat['X'], mat['y']
displayData( X )

def part1_3():
mat = scipy.io.loadmat( "ex3data1.mat" )
X, y = mat['X'], mat['y']
m, n = shape( X )
input_layer_size = 400
num_labels = 10
lamda = 0.1

theta = oneVsAll( X, y, num_labels, lamda )

def part1_4():
mat = scipy.io.loadmat( "ex3data1.mat" )
X, y = mat['X'], mat['y']
m, n = shape( X )
input_layer_size = 400
num_labels = 10
lamda = 0.1

theta = oneVsAll( X, y, num_labels, lamda )
predictOneVsAll( theta, X, y )

displayData( X, theta )


def main():
set_printoptions(precision=6, linewidth=200)
part1_1()
part1_2()
part1_3()
part1_4()

if __name__ == '__main__':
main(