#!/usr/bin/python

import numpy
numpy.random.seed(1337)

from keras.datasets import mnist
from keras.models import Sequential
from keras.layers.core import Dense, Dropout, Activation
from keras.optimizers import Adam, RMSprop, SGD
from keras.utils import np_utils

batch_size = 128
nb_epoch = 6

trn_X = numpy.loadtxt("mnist_trn_X.csv", delimiter = ",")
tst_X = numpy.loadtxt("mnist_tst_X.csv", delimiter = ",")
trn_y = numpy.loadtxt("mnist_trn_y.txt", delimiter = ",")
tst_y = numpy.loadtxt("mnist_tst_y.txt", delimiter = ",")
nb_classes = len(numpy.unique(trn_y))
trn_X = trn_X.astype("float32")
tst_X = tst_X.astype("float32")

# normalizing to convert 0..255 to 0..1
trn_X /= 255
tst_X /= 255
print(trn_X.shape[0], " trn samples")
print(tst_X.shape[0], " tst samples")

# one hot encoding
trn_Y = np_utils.to_categorical(trn_y, nb_classes)
tst_Y = np_utils.to_categorical(tst_y, nb_classes)

# 2 layers of rectified linear units and a softmax output layer
model = Sequential()
model.add(Dense(512, input_shape = (trn_X.shape[1],)))
model.add(Activation("relu"))
model.add(Dropout(0.2))
model.add(Dense(512))
model.add(Activation("relu"))
model.add(Dropout(0.2))
model.add(Dense(10))
model.add(Activation("softmax"))

print(model.summary())

model.compile(loss = "categorical_crossentropy",
              optimizer = RMSprop(),
              metrics = [ "accuracy" ])
training = model.fit(trn_X, trn_Y,
                     batch_size=batch_size,
                     nb_epoch=nb_epoch,
                     verbose=1,
                     validation_data=(tst_X, tst_Y))
print(training.history["val_acc"])

score = model.evaluate(tst_X, tst_Y, verbose = 0)
print("Test loss:", score[0])
print("Test accuracy:", score[1])
import math
print(math.exp(- score[0]))

predictions = model.predict(tst_X)
output = open("predictions.csv", "w")
output.write("index," + ",".join([ "probability" + str(i) for i in range(predictions.shape[1]) ]) + "\n")
for i in range(predictions.shape[0]):
    output.write(str(i) + "," + ",".join([ str(prediction) for prediction in predictions[i,:] ]) + "\n")
output.close()

total = 0.0
for i in range(predictions.shape[0]):
    total = total - numpy.log(predictions[i,int(tst_y[i])])
print(total)