{"cells":[{"metadata":{"_uuid":"37b84883ea624627949fd0796fa052111f998590"},"cell_type":"markdown","source":"#  Low Parameter Inception Model for MNIST\n\nThis kernel provides an overview of a low parameter Inception model that is capable of 99.7% test accuracy in Kaggle's [Digit Recognizer competition](https://www.kaggle.com/c/digit-recognizer).\n\nI have a slow computer but I still wanted to work on deep learning competitions. \nTo make things interesting (and feasible with my computer), I tried to create the lowest parameter model I could for MNIST.\nThe end result is an Inception model that uses 76,264 parameters.\nAs far as I can tell, other models that get state of the art results on MNIST all use 200,000 or more parameters.\n\n- [Load and Prepare Data](#Load-and-Prepare-Data)   \n\n- [Inception](#Inception)\n    - [Inception Modules](#Inception-Modules)\n    - [Inception Model](#Inception-Model)  \n\n- [Training](#Training)\n    - [Dropout](#Dropout)\n    - [Label Smoothing](#Label-Smoothing)\n    - [Data Augmentation](#Data-Augmentation)\n    - [Cosine Annealing](#Cosine-Annealing)\n    - [Adamax](#Adamax)   \n    - [Train](#Train)\n\n- [Model Evaluation](#Model-Evaluation) \n\n- [Predict on Test Set](#Predict-on-Test-Set)"},{"metadata":{"trusted":true,"_uuid":"79b9c6fd40494e6d06605a787908762bb2e2b12e"},"cell_type":"code","source":"%matplotlib inline\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport time\n\nimport keras\nfrom keras import backend as K\nfrom keras import layers\nfrom keras import optimizers\nfrom keras.models import Model\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.utils import to_categorical\n\n\nkeras_version = keras.__version__\ntf_version = K.tensorflow_backend.tf.VERSION\n\nprint(\"keras version:\", keras_version)\nprint(K.backend(), \"version:\", tf_version)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"30238db9bb4087ba2ca12c09d89f0de53827ce77"},"cell_type":"markdown","source":"## Load and Prepare Data\n\nThis part is pretty standard but there are a few things to note. The pixel values are scaled to the range $[-1, 1]$. The validation set is selected with a fixed seed so as to be reproducible. I use a validation set that is balanced across classes because it makes it easier to evaluate model behavior on MNIST&mdash;particularly when looking at misclassified examples and class bias."},{"metadata":{"trusted":true,"_uuid":"735d0822a68ebbd6b177567008324113e7e59a3a"},"cell_type":"code","source":"# load data\nrawdata = np.loadtxt('../input/train.csv', dtype=int, delimiter=',', skiprows=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"20fc9ff6ce487027fbb05b1bbc64e197a28ff9d0"},"cell_type":"code","source":"# inspect data\nprint(\"Raw data shape:\", rawdata.shape)\n\n# split labels and pixel values\ny = rawdata[:, 0]\nX = rawdata[:, 1:]\nprint(\"Labels shape:\", y.shape)\nprint(\"Pixels shape:\", X.shape)\n\n# convert pixel values to 2d arrays\nX = np.reshape(X, (42000, 28, 28))\nprint(\"Pixels reshaped shape:\", X.shape)\n\n# display random sample of images\nplt.rcParams['figure.figsize'] = (10.0, 4.0)\nplt.rcParams['image.interpolation'] = 'nearest'\nplt.rcParams['image.cmap'] = 'gray'\n\nnum_classes = 10\nsamples_per_class = 4\nfor cls in range(num_classes):\n    idxs = np.flatnonzero(y == cls)\n    idxs = np.random.choice(idxs, samples_per_class, replace=False)\n    for i, idx in enumerate(idxs):\n        plt_idx = i * num_classes + cls + 1\n        plt.subplot(samples_per_class, num_classes, plt_idx)\n        plt.imshow(X[idx])\n        plt.axis('off')\n        if i == 0:\n            plt.title(cls)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9efb8fe4ba56a2596f3a415e37f38bfd0cc6f804"},"cell_type":"code","source":"# one hot encode labels\ny_oh = to_categorical(y, num_classes)\n\n# scale pixel values to be between -1 and 1\nX_scaled = X / 127.5 - 1\nX_scaled = np.expand_dims(X_scaled, -1) # channels last\n\n# split data into train set and balanced validation set\nnum_val = int(y.shape[0] * 0.1)\nvalidation_mask = np.zeros(y.shape[0], np.bool)\nnp.random.seed(1)\nfor c in range(num_classes):\n    idxs = np.random.choice(np.flatnonzero(y == c), num_val // 10, replace=False)\n    validation_mask[idxs] = 1\nnp.random.seed(None)  \n    \nX_train = X_scaled[~validation_mask]\nX_val = X_scaled[validation_mask]\nprint(\"Train/val pixel shapes:\", X_train.shape, X_val.shape)\n\ny_train = y_oh[~validation_mask]\ny_val = y_oh[validation_mask]\nprint(\"Train/val label shapes:\", y_train.shape, y_val.shape)\n\n# confirm validation set is balanced across classes\nprint(\"Validation Set Class Distribution:\", np.bincount(y[validation_mask]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a755caabaf3d31da3cc965d327055b29948ac0d3"},"cell_type":"markdown","source":"## Inception\n\nThe inception model was introduced in [\"Going deeper with convolutions\"](https://arxiv.org/abs/1409.4842) (Szegedy et al, 2014).\nThe model evolved in [\"Rethinking the Inception Architecture for Computer Vision\"](https://arxiv.org/abs/1512.00567) (Szegedy et al, 2015) and mature models were outlined in [\"Inception-v4, Inception-ResNet and the Impact of Residual Connections on Learning\"](https://arxiv.org/abs/1602.07261) (Szegedy et al, 2016). \nInception modules are the core of an inception model, but the way the modules and other pieces fit together is also important.\nThe current basic inception module consists of:\n- Batch Normalized Convolutions\n- 1 x 1 Convolutions for Dimensionality Reduction\n- Internal Branches with Different Spatial Coverages\n- Pooling Branch\n- Some Modules use Asymmetric Convolutions\n\nThe basic model format is as follows:\n\n`Stem` &rarr; `Inception Module Stacks` &rarr; `Global Average Pooling` &rarr; `Softmax`,\n\nwhere the Stem is convolutional and pooling layers and the Inception Module Stacks include reduction modules/pooling. \nInception models lend themselves well to low parameter models because:\n- The 1 x 1 convolutions do efficient dimensionality reduction which means the next convolution uses fewer parameters for the incoming data.\n- Asymmetric convolutions in modules use fewer parameters for the same spatial coverage as a symmetric convolution.\n- Global average pooling results in fewer parameters than flattening.\n- There are no hidden layers."},{"metadata":{"_uuid":"cc0f83d3ebc46b564a4d0a9119baba2e00f55615"},"cell_type":"markdown","source":"### Inception Modules"},{"metadata":{"trusted":true,"_uuid":"5f1b5d5e2baef58006756c9eb89822ed84fca856"},"cell_type":"code","source":"def conv2D_bn_relu(x, filters, kernel_size, strides, padding='valid', kernel_initializer='glorot_uniform', name=None):\n    \"\"\"2D convolution with batch normalization and ReLU activation.\n    \"\"\"\n    \n    x = layers.Conv2D(filters=filters, \n                      kernel_size=kernel_size, \n                      strides=strides, \n                      padding=padding, \n                      kernel_initializer=kernel_initializer,\n                      name=name,\n                      use_bias=False)(x)\n    x = layers.BatchNormalization(scale=False)(x)\n    return layers.Activation('relu')(x)\n\n\ndef inception_module_A(x, filters=None, kernel_initializer='glorot_uniform'):\n    \"\"\"Inception module A as described in Figure 4 of \"Inception-v4, Inception-ResNet \n    and the Impact of Residual Connections on Learning\" (Szegedy, et al. 2016).\n    \n    # Arguments\n        x: 4D tensor with shape: `(batch, rows, cols, channels)`.\n        filters: Number of output filters for the module.\n        kernel_initializer: Weight initializer for all convolutional layers in module.\n    \"\"\"\n    \n    if filters is None:\n        filters = int(x.shape[-1])\n    branch_filters = filters // 4\n        \n    b1 = conv2D_bn_relu(x, \n                        filters=(branch_filters // 3) * 2, \n                        kernel_size=1, \n                        strides=1, \n                        kernel_initializer=kernel_initializer)\n    b1 = conv2D_bn_relu(b1, \n                        filters=branch_filters, \n                        kernel_size=3, \n                        strides=1, \n                        padding='same', \n                        kernel_initializer=kernel_initializer)\n    \n    b2 = conv2D_bn_relu(x, \n                        filters=(branch_filters // 3) * 2, \n                        kernel_size=1, \n                        strides=1, \n                        kernel_initializer=kernel_initializer)\n    b2 = conv2D_bn_relu(b2, \n                        filters=branch_filters, \n                        kernel_size=3, \n                        strides=1, \n                        padding='same', \n                        kernel_initializer=kernel_initializer)\n    b2 = conv2D_bn_relu(b2, \n                        filters=branch_filters, \n                        kernel_size=3, \n                        strides=1, \n                        padding='same', \n                        kernel_initializer=kernel_initializer)\n        \n    b3 = conv2D_bn_relu(x, \n                        filters=branch_filters, \n                        kernel_size=1, \n                        strides=1, \n                        kernel_initializer=kernel_initializer)\n    \n    pool = layers.AveragePooling2D(pool_size=(3, 3), strides=1, padding='same')(x)\n    pool = conv2D_bn_relu(pool, \n                          filters=branch_filters, \n                          kernel_size=1, \n                          strides=1, \n                          kernel_initializer=kernel_initializer)\n\n    return layers.concatenate([b1, b2, b3, pool])\n\n\ndef inception_module_C(x, filters=None, kernel_initializer='glorot_uniform'):\n    \"\"\"Inception module C as described in Figure 6 of \"Inception-v4, Inception-ResNet \n    and the Impact of Residual Connections on Learning\" (Szegedy, et al. 2016).\n    \n    # Arguments\n        x: 4D tensor with shape: `(batch, rows, cols, channels)`.\n        filters: Number of output filters for the module.\n        kernel_initializer: Weight initializer for all convolutional layers in module.\n    \"\"\"\n        \n    if filters is None:\n        filters = int(x.shape[-1])\n    branch_filters = filters // 6\n        \n    b1 = conv2D_bn_relu(x, \n                        filters=(branch_filters // 2) * 3, \n                        kernel_size=1, \n                        strides=1, \n                        kernel_initializer=kernel_initializer)\n        \n    b1a = conv2D_bn_relu(b1, \n                         filters=branch_filters, \n                         kernel_size=(1, 3), \n                         strides=1, \n                         padding='same', \n                         kernel_initializer=kernel_initializer)\n    \n    b1b = conv2D_bn_relu(b1, \n                         filters=branch_filters, \n                         kernel_size=(3, 1), \n                         strides=1, \n                         padding='same', \n                         kernel_initializer=kernel_initializer)\n    \n    b2 = conv2D_bn_relu(x, \n                        filters=(branch_filters // 2) * 3, \n                        kernel_size=1, \n                        strides=1, \n                        kernel_initializer=kernel_initializer)\n    b2 = conv2D_bn_relu(b2, \n                        filters=(branch_filters // 4) * 7, \n                        kernel_size=(1, 3), \n                        strides=1, \n                        padding='same', \n                        kernel_initializer=kernel_initializer)\n    b2 = conv2D_bn_relu(b2, \n                        filters=branch_filters * 2, \n                        kernel_size=(3, 1), \n                        strides=1, \n                        padding='same', \n                        kernel_initializer=kernel_initializer)\n\n    b2a = conv2D_bn_relu(b2, \n                         filters=branch_filters, \n                         kernel_size=(1, 3), \n                         strides=1, \n                         padding='same', \n                         kernel_initializer=kernel_initializer)\n    \n    b2b = conv2D_bn_relu(b2, \n                         branch_filters, \n                         kernel_size=(3, 1), \n                         strides=1, \n                         padding='same', \n                         kernel_initializer=kernel_initializer)\n        \n    b3 = conv2D_bn_relu(x, \n                        filters=branch_filters, \n                        kernel_size=1, \n                        strides=1, \n                        kernel_initializer=kernel_initializer)\n    \n    pool = layers.AveragePooling2D(pool_size=(3, 3), strides=1, padding='same')(x)\n    pool = conv2D_bn_relu(pool, \n                          filters=branch_filters, \n                          kernel_size=1, \n                          strides=1, \n                          kernel_initializer=kernel_initializer)\n    \n    return layers.concatenate([b1a, b1b, b2a, b2b, b3, pool])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f9386b1ef04e05bb50f18a3e689a7632412053fb"},"cell_type":"markdown","source":"### Inception Model"},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"867495ac5c1e7376b1bac6d2420d24944c07b3e5"},"cell_type":"code","source":"K.clear_session()\n\nstem_width = 64\n\ninputs = layers.Input(shape=X_scaled.shape[1:])\nx = conv2D_bn_relu(inputs,\n                   filters=stem_width,\n                   kernel_size=5,\n                   strides=2,\n                   padding='valid',\n                   name='conv_1')\n\nx = inception_module_A(x, filters=int(1.5*stem_width))\nx = layers.SpatialDropout2D(0.2)(x)\n\nx = inception_module_A(x, filters=int(1.5*stem_width))\nx = layers.MaxPooling2D(pool_size=(3, 3), strides=2, padding='same')(x)\nx = layers.SpatialDropout2D(0.2)(x)\n\nx = inception_module_C(x, filters=int(2.25*stem_width))\nx = layers.GlobalAveragePooling2D()(x)\nx = layers.Dropout(0.3)(x)\n\nx = layers.Dense(num_classes, name='logits')(x)\nx = layers.Activation('softmax', name='softmax')(x)\n\nmodel = Model(inputs=inputs, outputs=x)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"65d1c2cde49eaaf0263042b979901a91f1021787"},"cell_type":"markdown","source":"## Training"},{"metadata":{"_uuid":"351c5067b3fff91c5b10a4e932a8f91d5fec06b4"},"cell_type":"markdown","source":"### Dropout\n\nThe reader is probably familar with using [dropout](http://www.jmlr.org/papers/volume15/srivastava14a/srivastava14a.pdf) (Srivastava et al, 2014) to prevent overfitting. \nFor convolutional neural networks there is also [spatial dropout](https://arxiv.org/abs/1411.4280) (Tompson et al, 2014), which is similar to dropout but drops entire channels instead of individual coordinates. "},{"metadata":{"_uuid":"1825bfeb7ebe4652e59aec14eb55b1fd61f950ed"},"cell_type":"markdown","source":"### Label Smoothing\nThis is a regularization method from section *7. Model Regularization via Label Smoothing* in [\"Rethinking the Inception Architecture for Computer Vision\"](https://arxiv.org/abs/1512.00567) (Szegedy et al, 2015). I found that it resulted in lower validation loss compared to validation loss in models that did not use label smoothing, though I'm not sure if it helped validation accuracy."},{"metadata":{"trusted":true,"_uuid":"6fef148a8a03afeb0826c1e051e2ad610b46d010"},"cell_type":"code","source":"epsilon = 0.001\ny_train_smooth = y_train * (1 - epsilon) + epsilon / 10\nprint(y_train_smooth)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68df92e18b1e5fbc5bec4880fd7e89825a50ba20"},"cell_type":"markdown","source":"### Data Augmentation\n\nMy Kaggle kernel [MNIST Data Augmentation with Elastic Distortion](https://www.kaggle.com/babbler/mnist-data-augmentation-with-elastic-distortion) has an overview with visualizations of the data augmentation methods used here. To get the best results, I had to:\n- Use [elastic distortion](https://www.microsoft.com/en-us/research/wp-content/uploads/2003/08/icdar03.pdf) (Simard et al, 2003).\n- Shift images by pixel values to avoid interpolation blur.\n\n#### Elastic Distortion Function\n\nCredit to the following gists for the basic function:\n\n- https://gist.github.com/fmder/e28813c1e8721830ff9c\n- https://gist.github.com/chsasank/4d8f68caf01f041a6453e67fb30f8f5a\n- https://gist.github.com/erniejunior/601cdf56d2b424757de5"},{"metadata":{"trusted":true,"_uuid":"0c007a8dafa7602845369519edccc64ded241058"},"cell_type":"code","source":"from scipy.ndimage.filters import gaussian_filter\nfrom scipy.ndimage.interpolation import map_coordinates\n\ndef elastic_transform(image, alpha_range, sigma, random_state=None):\n    \"\"\"Elastic deformation of images as described in [Simard2003]_.\n    .. [Simard2003] Simard, Steinkraus and Platt, \"Best Practices for\n       Convolutional Neural Networks applied to Visual Document Analysis\", in\n       Proc. of the International Conference on Document Analysis and\n       Recognition, 2003.\n       \n   # Arguments\n       image: Numpy array with shape (height, width, channels). \n       alpha_range: Float for fixed value or [lower, upper] for random value from uniform distribution.\n           Controls intensity of deformation.\n       sigma: Float, sigma of gaussian filter that smooths the displacement fields.\n       random_state: `numpy.random.RandomState` object for generating displacement fields.\n    \"\"\"\n    \n    if random_state is None:\n        random_state = np.random.RandomState(None)\n        \n    if np.isscalar(alpha_range):\n        alpha = alpha_range\n    else:\n        alpha = np.random.uniform(low=alpha_range[0], high=alpha_range[1])\n\n    shape = image.shape\n    dx = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma) * alpha\n    dy = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma) * alpha\n\n    x, y, z = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]), np.arange(shape[2]), indexing='ij')\n    indices = np.reshape(x+dx, (-1, 1)), np.reshape(y+dy, (-1, 1)), np.reshape(z, (-1, 1))\n\n    return map_coordinates(image, indices, order=1, mode='reflect').reshape(shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e976e38bcc6301efd0780bcab632e26173e250b8"},"cell_type":"markdown","source":"### Cosine Annealing\n\n[Cosine annealing](https://arxiv.org/abs/1608.03983) (Loshchilov & Hutter, 2017) is a relatively new learning rate annealing technique that does a more thorough job of exploring the model's solution space by using warm restarts to break out of local minimums.\nAs the learning rate decreases, the model gets more precise but may also get stuck in a particular state. Warm restarts raise the learning rate to get the model unstuck. I found that, as long as the model doesn't overfit on the training set, continual warm restarts can potentially discover better and better models. "},{"metadata":{"trusted":true,"_uuid":"d1fafc5e97d86b1633cd75009b5aedb2855ce11e"},"cell_type":"code","source":"class CosineAnneal(keras.callbacks.Callback):\n    \"\"\"\"Cosine annealing with warm restarts.\n    \n    As described in section 3 of \"SGDR: Stochastic Gradient Descent with Warm Restarts\" (Loshchilov & Hutter 2017).\n    \n    # Arguments\n        max_lr: Maximum value of learning rate range.\n        min_lr: Minimum value of learning rate range.\n        T: Number of epochs between warm restarts.\n        T_mul: At warm restarts, multiply `T` by this amount.\n        decay_rate: At warm restarts, multiply the max_lr and min_lr by this amount.\n    \"\"\"\n    def __init__(self, max_lr, min_lr, T, T_mul=1, decay_rate=1.0):\n        self.max_lr = max_lr\n        self.min_lr = min_lr\n        self.decay_rate = decay_rate\n        self.T = T\n        self.T_cur = 0\n        self.T_mul = T_mul\n        self.step = 0\n        \n    def on_batch_begin(self, batch, logs=None):\n        if self.T <= self.T_cur:\n            self.max_lr *= self.decay_rate\n            self.min_lr *= self.decay_rate\n            self.T *= self.T_mul\n            self.T_cur = 0\n            self.step = 0\n        lr = self.min_lr + 0.5 * (self.max_lr - self.min_lr) * (1 + np.cos(self.T_cur * np.pi / self.T))        \n        K.set_value(self.model.optimizer.lr, lr)\n        # use self.step to avoid floating point arithmetic errors at warm restarts\n        self.step += 1\n        self.T_cur = self.step / self.params['steps']\n            \n    def on_epoch_end(self, epoch, logs=None):\n        logs = logs or {}\n        logs['lr'] = K.get_value(self.model.optimizer.lr)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"53b93a3e30a0d43c310d0ca12bd2dd5f095a63a1"},"cell_type":"markdown","source":"### Adamax\nThe Adamax optimizer was introduced alongside Adam in [\"Adam: A Method for Stochastic Optimization\"](https://arxiv.org/abs/1412.6980) (Kingma & Ba, 2014). While Adam is more popular, I found that Adamax is less prone to overfitting. Adamax also worked better than SGD with momentum when used in conjunction with cosine annealing."},{"metadata":{"trusted":true,"_uuid":"5c7726dfe7cc6a3c7b091899172fc1af79197c87"},"cell_type":"code","source":"model.compile(loss='categorical_crossentropy', \n              optimizer=optimizers.Adamax(lr=0.01, beta_1=0.49, beta_2=0.999),\n              metrics=['accuracy'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ab87edf633b715e9ffd58d7cad4a68f512cf82be"},"cell_type":"markdown","source":"### Train\n\nI save the best model for final evaluation and use enough epochs so cosine annealing can thoroughly explore the solution space."},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"aa5015ba5e9bb5a40d6840efdfba800dd7773a69"},"cell_type":"code","source":"batch_size = 32\nepochs = 100\ntime_id = time.strftime(\"%Y-%m-%d_%H-%M-%S\")\nbest_model_filename = 'mnist-inception-best-' + time_id + '.hdf5'\n\n# save weights for visualization\npre_train_weights = model.get_layer('conv_1').get_weights()[0]\npre_train_weights = pre_train_weights.transpose(3, 2, 0, 1)\n\n# setup callbacks\nannealer = CosineAnneal(max_lr=0.014, min_lr=0.003, T=5, T_mul=1, decay_rate=0.99)\nchkpt = keras.callbacks.ModelCheckpoint(best_model_filename, monitor='val_acc', \n                                        save_best_only=True, verbose=False)\n\n# define data augmentations\ndatagen = ImageDataGenerator(\n    width_shift_range=2,\n    height_shift_range=2,\n    preprocessing_function=lambda x: elastic_transform(x, alpha_range=[8, 10], sigma=3)\n)\n\n# train model\nhistory = model.fit_generator(\n    datagen.flow(X_train, y_train_smooth, batch_size=batch_size, shuffle=True),\n    epochs=epochs,\n    steps_per_epoch=(len(y_train) - 1) // batch_size + 1,\n    validation_data=(X_val, y_val),\n    callbacks=[annealer, chkpt]\n)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a22f6736a6f9713c62bd2148fde710ccfe2948c9"},"cell_type":"markdown","source":"## Model Evaluation"},{"metadata":{"trusted":true,"_uuid":"8a9498291283f2b5f1b1a6b5dca471c68181ec25"},"cell_type":"code","source":"def plot_training_log(log, acc=True, loss=True, lr=True, figsize=(10.0, 2.5)):\n    \"\"\"Plot training history.\n    \n    # Arguments\n        log: Dictionary of training history, same as ` History.history`  that  \n            is returned when fitting a Keras model. Should have records for 'acc', \n            'val_acc', 'loss', 'val_loss', and optionally 'lr'.\n        acc: if true, plot both 'acc' and 'val_acc' on one plot.\n        loss: if true, plot both 'loss' and 'val_loss' on one plot.\n        lr: if true, and if 'lr' is in log, plot both 'lr' and 'val_acc' on one plot.\n        figsize: size of each plot.\n    \"\"\"\n    \n    plt.rcParams['figure.figsize'] = figsize\n    max_val_acc_epoch = np.argmax(list(log['val_acc'])) + 1\n    epochs = range(1, len(log['acc']) + 1)\n    \n    def plot(ytype, ylabel, max_val_acc_epoch):\n        plt.axvline(x=max_val_acc_epoch, color='0.5', linestyle='--')\n        plt.plot(epochs, log[ytype], label='Train')\n        plt.plot(epochs, log['val_' + ytype], label='Validation')\n        plt.minorticks_on()\n        plt.grid(b=True, axis='x', which='both', color='0.8', linestyle='-')\n        plt.xlabel('Epoch')\n        plt.ylabel(ylabel)\n        plt.xlim(0, epochs[-1] + 1)\n        plt.legend()\n        plt.tight_layout()\n        plt.show()\n        \n    if acc:\n        plot('acc', 'Accuracy', max_val_acc_epoch)\n        \n    if loss:\n        plot('loss', 'Loss', max_val_acc_epoch)\n    \n    if lr and 'lr' in log.keys():\n        fig, ax1 = plt.subplots()\n        plt.axvline(x=max_val_acc_epoch, color='0.5', linestyle='--')\n\n        ln1 = ax1.plot(epochs, log['lr'], 'C4o-', label='Learning Rate')   \n        ax1.set_xticks(epochs, minor=True)\n        ax1.grid(b=True, axis='x', which='both', color='0.8', linestyle='-')\n        ax1.set_xlabel('Epoch')\n        ax1.set_ylabel('Learning Rate')\n        ax1.set_xlim(0, epochs[-1] + 1)\n\n        ax2 = ax1.twinx()\n        ln2 = ax2.plot(epochs, log['val_acc'], 'C1o-', \n                       label='Validation Accuracy [0.98 - 1.0]')\n        ax2.set_yticks([])\n        ax2.set_ylim(0.98, 1.0)        \n\n        lns = ln1 + ln2\n        plt.legend(lns, [l.get_label() for l in lns])\n        plt.tight_layout()\n        plt.show()\n        \n        \ndef evaluate_model(model, X, y, log=None, pre_train_weights=None):\n    \"\"\"Display accuracy, display misclassified digits, and visualize weights in \n    first convolution layer for the given model with X and y as prediction inputs.\n    Optionally, accuracy and loss is plotted with the log argument and additional\n    weights can be visualized with the pre_train_weights argument.\n    \n    # Arguments\n        model: Keras model with first convolution layer named 'conv_1'.\n        X: Numpy array of MNIST digits, in channels last format.\n        y: Labels for X.\n        log: (optional) Dictionary of training history, same as ` History.history`  that  \n            is returned when fitting a Keras model. Should have records for 'acc', \n            'val_acc', 'loss', 'val_loss', and optionally 'lr'.\n        pre_train_weights: (optional) Numpy array of weights in first convolution layer \n            of the model.\n    \"\"\"\n \n    if log is not None:\n        print(\"Max Validation Accuracy:\", max(log['val_acc']))\n        plot_training_log(log)\n\n    scores = model.predict(X)\n    predictions = np.argmax(scores, axis=1)\n    y_digits = np.nonzero(y)[1]       \n    print(\"Accuracy:\",  np.mean(predictions == y_digits))\n    \n    post_train_weights = model.get_layer('conv_1').get_weights()[0]\n    post_train_weights = post_train_weights.transpose(3, 2, 0, 1)\n    num_weights = len(post_train_weights)\n    if pre_train_weights is not None:\n        plt.rcParams['figure.figsize'] = (10.0, 3.0)\n        for i in range(num_weights):\n            plt.subplot(4, num_weights // 4, i + 1)\n            ker = pre_train_weights[i, 0]\n            low, high = np.amin(ker), np.max(ker)\n            plt.imshow(255 * (ker - low) / (high - low))\n            plt.axis('off')\n        plt.suptitle('Pre-Training Weights from First Convolutional Layer')\n        plt.show()\n    plt.rcParams['figure.figsize'] = (10.0, 3.0)\n    for i in range(num_weights):\n        plt.subplot(4, num_weights // 4, i + 1)\n        ker = post_train_weights[i, 0]\n        low, high = np.amin(ker), np.max(ker)\n        plt.imshow(255 * (ker - low) / (high - low))\n        plt.axis('off')\n    plt.suptitle('Post-Training Weights from First Convolutional Layer')\n    plt.show()\n\n    misclassified_mask = (predictions != y_digits)\n    samples_per_class = 7\n    num_classes = 10\n    plt.rcParams['figure.figsize'] = (10.0, 7.0)    \n    counts = [];\n    for cls in range(10):\n        idxs = np.flatnonzero(y_digits[misclassified_mask] == cls)\n        counts.append(len(idxs))\n        if len(idxs) > samples_per_class:\n            idxs = np.random.choice(idxs, samples_per_class, replace=False)       \n        for i, idx in enumerate(idxs):\n            plt_idx = i * num_classes + cls + 1\n            plt.subplot(samples_per_class, num_classes, plt_idx)\n            plt.imshow(X[misclassified_mask][idx,:,:,0])\n            plt.axis('off')\n            plt.text(14, 27, 'Pred:' + str(predictions[misclassified_mask][idx]), \n                     horizontalalignment='center', verticalalignment='top')        \n            if i == 0:\n                plt.title(cls)            \n    plt.suptitle('Misclassified Digits | Norm of Counts = ' + str(np.linalg.norm(counts)))\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d9a01f879c9bdbdf81331327dc72b89928f95f46"},"cell_type":"markdown","source":"#### How to Tell if the Model is Good:\n\nWhile this method is capable of generating a model that can get 99.7% test accuracy, it's not guaranteed.\nHere are some of the ways I evaluate models.\nI should note that the model that got 99.7% test accuracy only got 99.62% validation accuracy. However, it performed very well in all my other evaluation measures.\n\n- Model Behavior:\n\t1. Check that misclassification examples make sense. In most cases a human won't mistake a 3 for an 8. On the other hand, to me at least, 4 and 9 sometimes look similar. The model should have the same behavior.\n\t2. Avoid bias for or against particular classes. The norm of misclassified counts was a quick and dirty way to see this. I knew from experience with the data set that I wanted to see norms less than 8.\n- Training Metrics (loss & accuracy):\n\t1. Review the metric history from previous training runs and look for models with metrics in top X%.\n\t2. Look at metric noise over the training run&mdash;I found that huge jumps up and down from epoch to epoch would indicate the model was poor, even if validation accuracy was high. Note that I'm not refering to the jumps caused by warm restarts.\n\t3. Look for models with similar test and validation accuracy. Validation accuracy should be highr than test accuracy but not by much.\n\t4. The rate of convergence can be useful but I find it too subjective to give guidelines. Pay attention to how quickly the metrics change and eventually you'll get a feel for when a model is good.\n- Does the Model Improve with Additional Tuning Measures:\n\t1. Check for improvement with warm restarts&mdash;if there's still an upward trend in validation accuracy after 10 warm restarts, keep going.\n\t2. Check for improvement after fine tuning with SGD/Adamax and a low learning rate.\n\t3. Compare to an ensemble of models. For example, save multiple top performing models and use them to make ensembles."},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"05ee9db9637ca96addf94fce4ad0617898b49f31"},"cell_type":"code","source":"model = keras.models.load_model(best_model_filename)\nevaluate_model(model, X_val, y_val, log=history.history, pre_train_weights=pre_train_weights)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"292486ad3b0f7241a87ace9c0909294cfd9ede63"},"cell_type":"markdown","source":"## Predict on Test Set"},{"metadata":{"trusted":true,"_uuid":"13972602bfdd6e5a6a56f273242d9979e92fa148"},"cell_type":"code","source":"# load test set data\ntest_data = np.loadtxt('../input/test.csv', dtype=int, delimiter=',', skiprows=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9f2f42a6229af9ab0154777e3f819ab23993c2bf"},"cell_type":"code","source":"# inspect data\nprint(\"Raw data shape:\", test_data.shape)\n\n# convert pixel values to 2d arrays and scale data\ntest_images = np.reshape(test_data, (28000, 28, 28))\ntest_images = test_images / 127.5 - 1\ntest_images = np.expand_dims(test_images, -1)\nprint(\"Pixels reshaped shape:\", test_images.shape)\n\n# display some images\nplt.rcParams['figure.figsize'] = (10.0, 4.0)\nfor i in range(40):\n    plt.subplot(4, 10, i + 1)\n    plt.imshow(test_images[i,:,:,0])\n    plt.axis('off')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2cf86b02532a55a7b13fc73f5b4523d0481f713e"},"cell_type":"code","source":"# get predictions\ntest_scores = model.predict(test_images)\ntest_predictions = np.argmax(test_scores, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a13401d2df098a0182ee5e83c50ce8e86897253"},"cell_type":"code","source":"# save submission csv\nheader = 'ImageId,Label'\nsubmission = np.stack((range(1, 28001), test_predictions), axis=1)\nnp.savetxt('submission-' + time_id + '.csv', submission, fmt='%i', delimiter=',', header=header, comments='')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"widgets":{"application/vnd.jupyter.widget-state+json":{"state":{},"version_major":2,"version_minor":0}}},"nbformat":4,"nbformat_minor":1}