from theano.sandbox import cuda
cuda.use('gpu2')

%matplotlib inline
import utils; reload(utils)
from utils import *
from __future__ import division, print_function

"""
Explanation: Digit Recognizer (MNIST Data)
This notebook is heavily based on Jeremy Howard's notebooks from Practical Deep Learning For Coders, Part 1 (course.fast.ai). Download the competitions files into data/mnist folder.
Kaggle Competition: https://www.kaggle.com/c/digit-recognizer
Score: 0.994
Imports
End of explanation
"""


path="data/mnist/"
model_path = path + 'models/'
results_path = path + 'results/'
submissions_path = path + 'submissions/'
if not os.path.exists(model_path): os.mkdir(model_path)
if not os.path.exists(results_path): os.mkdir(results_path)
if not os.path.exists(submissions_path): os.mkdir(submissions_path)
batch_size=64

"""
Explanation: Setup
End of explanation
"""


training = pd.read_csv(path+'train.csv')
training.head()

test = pd.read_csv(path+"test.csv")
test.head()

"""
Explanation: Reading Data
End of explanation
"""


Y_train = training['label'].values

Y_train

X_train = training.iloc[:,1:].values
X_train = X_train.reshape(X_train.shape[0], 28, 28, 1) #reshape to rectangular
X_train = X_train/255 #pixel values are 0 - 255 - this makes puts them in the range 0 - 1

X_val = X_train[32000:]

X_train = X_train[0:32000]

Y_val = Y_train[32000:]
Y_train = Y_train[0:32000]

X_test = test.iloc[:,0:].values
X_test = X_test.reshape(X_test.shape[0], 28, 28, 1) #reshape to rectangular
X_test = X_test/255 #pixel values are 0 - 255 - this makes puts them in the range 0 - 1

(X_train.shape, Y_train.shape, X_val.shape, Y_val.shape)

"""
Explanation: Normalize
End of explanation
"""


X_train = np.expand_dims(X_train,1)
X_val = np.expand_dims(X_val,1)

X_train = np.squeeze(X_train,4)
X_val = np.squeeze(X_val,4)

(X_train.shape, Y_train.shape, X_val.shape, Y_val.shape)

"""
Explanation: To match the axis that theano expects (channel on axis 1)
End of explanation
"""


Y_train = onehot(Y_train)
Y_val = onehot(Y_val)

Y_train[:5]

mean_px = X_train.mean().astype(np.float32)
std_px = X_train.std().astype(np.float32)

def norm_input(x): return (x-mean_px)/std_px

"""
Explanation: As expected from theano
End of explanation
"""


def plots(ims, interp=False, titles=None):
    ims=np.array(ims)
    mn,mx=ims.min(),ims.max()
    f = plt.figure(figsize=(12,24))
    for i in range(len(ims)):
        sp=f.add_subplot(1, len(ims), i+1)
        if not titles is None: sp.set_title(titles[i], fontsize=18)
        plt.imshow(ims[i], interpolation=None if interp else 'none', vmin=mn,vmax=mx)

def plot(im, interp=False):
    f = plt.figure(figsize=(3,6), frameon=True)
    plt.imshow(im, interpolation=None if interp else 'none')
    
plt.gray()
plt.close()

digits = np.squeeze(X_train, 1)

digits.shape

plots(digits[0:8], titles=Y_train.argmax(axis=1)[0:8])

"""
Explanation: A look at the data
End of explanation
"""


def get_model_bn_do():
    model = Sequential([
        Lambda(norm_input, input_shape=(1,28,28)),
        Convolution2D(32,3,3, activation='relu'),
        BatchNormalization(axis=1),
        Convolution2D(32,3,3, activation='relu'),
        MaxPooling2D(),
        BatchNormalization(axis=1),
        Convolution2D(64,3,3, activation='relu'),
        BatchNormalization(axis=1),
        Convolution2D(64,3,3, activation='relu'),
        MaxPooling2D(),
        Flatten(),
        BatchNormalization(),
        Dense(512, activation='relu'),
        BatchNormalization(),
        Dropout(0.5),
        Dense(10, activation='softmax')
        ])
    model.compile(Adam(), loss='categorical_crossentropy', metrics=['accuracy'])
    return model

"""
Explanation: Batchnorm + dropout + data augmentation
End of explanation
"""


def fit_model():
    model = get_model_bn_do()
    model.fit_generator(batches, batches.N, nb_epoch=1,
                    validation_data=val_batches, nb_val_samples=val_batches.N)
    model.optimizer.lr=0.1
    model.fit_generator(batches, batches.N, nb_epoch=4,
                    validation_data=val_batches, nb_val_samples=val_batches.N)
    model.optimizer.lr=0.01
    model.fit_generator(batches, batches.N, nb_epoch=12,
                    validation_data=val_batches, nb_val_samples=val_batches.N)
    model.optimizer.lr=0.001
    model.fit_generator(batches, batches.N, nb_epoch=18,
                    validation_data=val_batches, nb_val_samples=val_batches.N)
    return model

models = [fit_model() for i in range(6)]

for i,m in enumerate(models):
    m.save_weights(model_path+'ensmb-'+str(i)+'.h5')

evals = np.array([m.evaluate(X_val, Y_val, batch_size=256) for m in models])

evals.mean(axis=0)

all_preds = np.stack([m.predict(X_val, batch_size=256) for m in models])

all_preds.shape

avg_preds = all_preds.mean(axis=0)

keras.metrics.categorical_accuracy(Y_val, avg_preds).eval()

"""
Explanation: Ensembling
End of explanation
"""


X_test.shape

X_test = np.expand_dims(X_test,1)
X_test = np.squeeze(X_test,4)

all_preds = np.stack([m.predict(X_test, batch_size=256) for m in models])

predictions = all_preds.mean(axis=0)

"""
Explanation: So we should have around 99.5% accuracy. Let's see:
Prediction Time
End of explanation
"""


