{"cells":[{"metadata":{},"cell_type":"markdown","source":"# MNIST digit classification with advanced CNN\n\nThe purpose of this kernel is to explore how CNN behave with MNIST?  \nShort answer : very well with the right layer setup\n\nThe best model achieve a 0.9962 accuracy on test set with 30 epochs / 5 folds and 0.997 with 200 epochs / 10 folds using this layer setup:\n-  [C32-C32-BN-MP-Dr0.25]-[C64-C64-BN-MP-Dr0.25]-[C128-C128-BN-MP-Dr0.25]-D300-D10\n\nWhere the layers are abbreviated using *Chris Deotte* convention with  \n- C32 = Conv2D(32)\n- BN = Batch Normalization\n- MP = MaxPoolin\n- Dr0.25  = Dropout(0.25)\n- D300 = Dense(300)\n\nThanks to the many public kernels and online ressources where inspiration has been taken from.\n\nAditionally a very good description of accuracy reachable for various models can be found on [How to score 97%, 98%, 99%, and 100%](https://www.kaggle.com/c/digit-recognizer/discussion/61480)\n\n# Dataset load"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"\nimport math\nimport time\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set()\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import confusion_matrix, accuracy_score\n\nimport os\nprint(os.listdir(\"../input\"))\n","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/train.csv\")\ntrain.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv(\"../input/test.csv\")\ntest.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_width  = 28\nimg_height = 28\nnclasses   = 10\nclasses    = np.arange(10)\nlabels     = [ str(c) for c in classes ]\ntrain_X = train.iloc[:,1:].values.reshape(train.shape[0], img_width, img_height, 1)\ntrain_y = train.iloc[:,0].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_X = test.iloc[:,0:].values.reshape(test.shape[0], img_width, img_height, 1)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Digit distribution among examples"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 1, figsize=(10,4))\nsns.countplot(train_y, ax=ax)\ntxt = ax.set_title(\"Train Digits distribution, std(count)={0:0.2f}\".format(np.std(np.histogram(train_y, np.arange(11))[0])))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Values are well distributed although with some drift between digits, notably the 5 is underrepresented"},{"metadata":{},"cell_type":"markdown","source":"# Input visualization\n\n- green number is the digit train label"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_digits(X, y, sel, predict=None, suptitle=\"\"):\n    nrow = math.ceil(len(sel)/20)\n    fig, axes = plt.subplots(nrow, 20, figsize=(14, nrow*1.2))\n    for i,ax in enumerate(axes.flatten()):\n        ax.imshow(X[sel[i],:,:].squeeze(), cmap='gray')\n        ax.set_xticks([])\n        ax.set_yticks([])\n        if y is not None:\n            ax.text(0, -2, str(y[sel[i]]), color=\"white\", bbox=dict(facecolor='green', alpha=0.5))\n        if predict is not None:\n            ax.text(img_width/2, -2, str(predict[sel[i]]), color=\"white\", bbox=dict(facecolor='purple', alpha=0.5))\n    fig.subplots_adjust(hspace=0, wspace=0)\n    fig.suptitle(suptitle)\n\ndef get_shuffle(n, seed=12):\n    rstate = np.random.RandomState(seed)\n    sel = np.arange(n)\n    rstate.shuffle(sel)\n    return sel\n    \n    \n# TRAIN\ntrain_sel = get_shuffle(len(train_X))\nplot_digits(train_X, train_y, train_sel[0:60], suptitle=\"Train digits\")\n\n# TEST\ntest_sel = get_shuffle(len(test_X))\nplot_digits(test_X, None, test_sel[0:60], suptitle=\"Test digits\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Normalization step1: check intensity distribution\n\nFirst let check pixels intensity distribution along  and Y axis"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,4))\nfor digit in range(10):\n    ax1.plot( (train_X[train_y==digit]**2).mean(axis=(0,1,3)), label=str(digit) )\n    ax2.plot( (train_X[train_y==digit]**2).mean(axis=(0,2,3)), label=str(digit) )\nax1.set_title(\"TRAIN: mean digit Power distribution over X axis\")\nax1.set_xlabel(\"X\")\nax2.set_title(\"TRAIN: mean digit Power distribution over Y axis\")\nax2.set_xlabel(\"Y\")\nleg = ax1.legend()\nleg = ax2.legend()\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,4))\nax1.plot( (test_X**2).mean(axis=(0,1,3)), label=\"TEST\" )\nax1.plot( (train_X**2).mean(axis=(0,1,3)), label=\"TRAIN\" )\nax2.plot( (test_X**2).mean(axis=(0,2,3)), label=\"TEST\" )\nax2.plot( (train_X**2).mean(axis=(0,2,3)), label=\"TRAIN\" )\nax1.set_title(\"TEST: mean Power distribution over X axis\")\ntext = ax1.set_xlabel(\"X\")\nax2.set_title(\"TEST: mean Power distribution over Y axis\")\ntext = ax2.set_xlabel(\"Y\")\nleg = ax1.legend()\nleg = ax2.legend()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The distribution looks very different among digits of the training set.\nThis is very encouraging for classification. It also explains why strategies like linear regressions and random forrest are very successfull.\n\nBoth sets have a surprisingly similar envelope distribution: this is also a good hint for generalization!\n\n\nLet have a look at pixel values distribution:"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(14,3))\ntrain_count, train_bins = np.histogram(train_X, range(256))\nplt.plot(train_bins[1:-1], train_count[1:], label=\"train\")\n\ntest_count, test_bins = np.histogram(test_X, range(256))\nplt.plot(test_bins[1:-1], test_count[1:], label=\"test\")\nplt.title(\"pixel values distribution\")\n\nleg = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The zero value has been omitted as most of the pixels are black. "},{"metadata":{},"cell_type":"markdown","source":"# Max normalization\n\nThe data has a clean distribution in pixel space and a spiked distribution in pixel values, but in the first try I noticed a very low accuracy rate of 0.1 which was really surprising.  \nAfter checking, the dtype input was `int64` and setting it to `float64` by using max normalization got accuracy over 0.97 on the first epoch with the Model 1 (see bellow) simple CNN :-)"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X = train_X / train_X.max()\ntest_X = test_X / test_X.max()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Keras CNN Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Dense, Flatten, Conv2D, Dropout, BatchNormalization, MaxPool2D\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom keras.utils import to_categorical\nfrom tensorflow import set_random_seed\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class MNIST_CNN:\n    \"\"\"All models have a final Dense(300) layer with Dropout\n    Model 1: [Conv2d(16/32/64) + Dropout(0.5)] x 3 \n        => 0.991 : wow simple is beautiful!\n    Model2: data aug(rot/shift/zoom) +  [Conv2d(16/32/64) + Dropout(0.5)] x 3  \n        => 0.994\n    Model3: data aug(rot/shift/zoom) +  [Conv2d(16/32/64) + BatchNorm + MaxPool + Dropout(0.25)] x 3\n        => 0.995\n    Model3: data aug(rot/shift/zoom) +  [Conv2d(32/64/128) + BatchNorm + MaxPool + Dropout(0.25)] x 3\n        => 0.99657\n    \"\"\"\n    def build(self, img_rows, img_cols, nclasses=10):\n        input_shape=(img_rows, img_cols, 1)\n        model = Sequential()\n        \n        # Suite of 3x Conv2D + BatchNorm + Dropout\n        model.add(Conv2D(32, (3,3), activation=\"relu\", input_shape=(img_rows, img_cols, 1),\n                        kernel_initializer=\"glorot_normal\", padding=\"same\"))\n        model.add(Conv2D(32, (3,3), activation=\"relu\",\n                        kernel_initializer=\"glorot_normal\", padding=\"same\"))\n        model.add(BatchNormalization())\n        model.add(MaxPool2D(pool_size=(2,2), strides=(2,2)))\n        model.add(Dropout(0.25))\n\n        model.add(Conv2D(64, (3,3), activation=\"relu\",\n                        kernel_initializer=\"glorot_normal\", padding=\"same\"))\n        model.add(Conv2D(64, (3,3), activation=\"relu\",\n                        kernel_initializer=\"glorot_normal\", padding=\"same\"))\n        model.add(BatchNormalization())\n        model.add(MaxPool2D(pool_size=(2,2), strides=(2,2)))\n        model.add(Dropout(0.25))\n        \n        model.add(Conv2D(128, (3,3), activation=\"relu\",\n                        kernel_initializer=\"glorot_normal\", padding=\"same\"))\n        model.add(Conv2D(128, (3,3), activation=\"relu\",\n                        kernel_initializer=\"glorot_normal\", padding=\"same\"))\n        model.add(BatchNormalization())\n        model.add(MaxPool2D(pool_size=(2,2), strides=(2,2)))\n        model.add(Dropout(0.25))\n      \n        # Final layers\n        model.add(Flatten())\n        model.add(Dense(300, activation=\"relu\"))\n        model.add(Dropout(0.25))\n\n        model.add(Dense(nclasses, activation=\"softmax\"))\n        self.model = model\n    \n    def compile(self, **kw):\n        kw1 = dict(optimizer=\"adam\", loss=\"categorical_crossentropy\", metrics=[\"categorical_accuracy\"])\n        kw1.update(kw)\n        self.model.compile(**kw1)\n        \n    def fit(self, *arg, **kw):\n        return self.model.fit(*arg, **kw)\n        \nmdl = MNIST_CNN()\nmdl.build(img_width, img_height)\nmdl.compile()\nmdl.model.summary()\n        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Use data augmentation\n\nImplementing data augmentation with ImageDataGenerator adds a continuous flow of slightly modified images from the original inputs. \n\nWhen enabling data augmentation it has the following effect that \n- the cross validation accuracy is actually better than the trainig set    \n\nThis may be explained by the fact that the network does never see twice the same training input, however generalization becomes better while averaging the multiple inputs.\nOn the contrary the cross validation set is never augmented, hence it is easier to classify.\n\nOne open question is how much similar should be both metrics in order to be otimal?"},{"metadata":{"trusted":true},"cell_type":"code","source":"datagen = ImageDataGenerator(\n    rotation_range=30,\n    #shear_range=30,\n    zoom_range=[0.9, 1.1],\n    width_shift_range = 2,\n    height_shift_range = 2)\ndatagen.fit(train_X)\n\naug_X, aug_y = datagen.flow(train_X, train_y, batch_size=60).next()\nplot_digits(aug_X, aug_y, np.arange(60), suptitle=\"Data Augmented digits\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Use a stratified K-Fold validation strategy\n\nK-Fold has the advantage that it provides a reliable estimator of model generalization on test data."},{"metadata":{"trusted":true,"_kg_hide-output":true,"collapsed":true,"_kg_hide-input":true},"cell_type":"code","source":"DEBUG = 0\nK = 3 if DEBUG else 10\nset_random_seed(777)\nskf = StratifiedKFold(n_splits=K, shuffle=True)\ntrain_y_cat = to_categorical(train_y)\nhistories = []\ntest_pred  = np.zeros((len(test_X), nclasses), dtype=np.float64)\ntrain_pred = np.zeros((len(train_X), nclasses), dtype=np.float64)\nval_pred   = np.zeros((len(train_X), nclasses), dtype=np.float64)\ntrain_accs = []\nval_accs   = []\n\nfor i, (train_idx, val_idx) in enumerate(skf.split(train_X, train_y)):\n    t0 = time.time()\n    batch_size=256 if DEBUG else 64\n    mdl = MNIST_CNN()\n    mdl.build(img_width, img_height)\n    mdl.compile()\n    histories.append(\n        mdl.model.fit_generator(\n            datagen.flow(train_X[train_idx], train_y_cat[train_idx], batch_size=batch_size),\n            epochs=5 if DEBUG else 200,\n            verbose=1 if DEBUG else 0,\n            steps_per_epoch=len(train_idx) / batch_size,\n            validation_data=(train_X[val_idx], train_y_cat[val_idx]), \n        ),\n    )\n    test_pred  += mdl.model.predict(test_X)\n    train_pred += mdl.model.predict(train_X)\n    val_pred[val_idx] += mdl.model.predict(train_X[val_idx])\n    train_accs.append( accuracy_score(train_y, mdl.model.predict(train_X).argmax(axis=1)) )\n    val_accs.append( accuracy_score(train_y[val_idx], val_pred[val_idx].argmax(axis=1)) )\n    print(\"******* [dt={0}] accuracy_scores at fold {1}: train={2} val={3}\".format(\n        time.time() - t0, i + 1, train_accs[-1], val_accs[-1] ) )\n    \ntest_y = np.argmax(test_pred, axis=1)\nval_y  = np.argmax(val_pred, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"global fit accuracy_scores: train={0} val={1}\".format(\n    accuracy_score(train_y, train_pred.argmax(axis=1)),\n    accuracy_score(train_y, val_pred.argmax(axis=1)) ))\n\nprint(f\"{K}-Fold train accuracy = {np.mean(train_accs):0.4f} +- {np.std(train_accs):0.4f}\")\nprint(f\"{K}-Fold cross validation accuracy = {np.mean(val_accs):0.4f} +- {np.std(val_accs):0.4f}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Learning metrics\n\nBellow are shown the metrics for train and validation. The variance of each metric is estimated with the k-fold runs. The enveloppe shows the $25$ and $75$ percentiles."},{"metadata":{"trusted":true,"collapsed":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_metrics(histories, metrics, ax, \n                 agg=lambda x: np.median(x, axis=0), \n                 env0=lambda x: np.percentile(x, 25, axis=0), \n                 env1=lambda x: np.percentile(x, 75, axis=0), \n                 process=lambda x: x, colors=None):\n    x = histories[0].epoch\n    for i, metric in enumerate(metrics):\n        # convert metric histories for each k-fold into an (nfold, epochs) array \n        data = np.asarray([h.history[metric] for h in histories])\n        data = process(data)\n        ax.fill_between(x, env0(data), env1(data), \n                        color=colors[i] if colors else None, alpha=0.5)\n        ax.plot(x, agg(data), color=colors[i] if colors else None, label=metric)\n        \nfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14,8))\n\nplot_metrics(histories, ['loss', 'val_loss'], ax1, colors=[\"green\", \"orange\" ])\nplot_metrics(histories, ['categorical_accuracy', 'val_categorical_accuracy'], \n             ax2, colors=[\"green\", \"orange\" ], process=lambda x: 1-x)\n\nax1.set_yscale('log'); ax1.set_ylabel('Loss'); ax1.legend(); ax1.set_ylim(1e-3, 1);\nax2.set_yscale('log'); ax2.set_ylabel(\"Error = 1-acc\");  ax2.legend(); ax2.set_ylim(1e-3, 1e-1);\nax2.set_xlabel('Epoch')\nfig.subplots_adjust(hspace=0.02)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A noticeable fact is that the error levels are higher than the value obtained above corresponding to the *global* accuracy which is obtained by averaging the predictions over all folds.\n\nWhen epochs is 30 the plateau is clearly not reached, meaning that an increase of about 0.001 in accuracy can be expected for values of epoch over 100. \n\nFor recall, in this setup, an epoch correspond to a full training set. Hence  to reach the full potential of data augmentation we need many epochs."},{"metadata":{},"cell_type":"markdown","source":"# Confusion matrix"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_confusion(ytrue, ypred, label=\"\"):\n    cnf_matrix = confusion_matrix(ytrue, ypred)\n    cnf_matrix_norm = cnf_matrix.astype('float') / cnf_matrix.sum(axis=1)[:, np.newaxis]\n    df_cm = pd.DataFrame(cnf_matrix_norm, index=classes, columns=classes)\n\n    plt.figure(figsize=(14, 9))\n    ax = plt.axes()\n    sns.heatmap(df_cm, annot=True, fmt='.4f', cmap=\"flag_r\", ax=ax)\n    ax.set_title(f\"Confusion matrix for {label}\")\n    ax.set_ylabel(\"True digit\")\n    ax.set_xlabel(\"Predicted digit\")\n    plt.show()\n\nplot_confusion(train_y, train_pred.argmax(axis=1), label=\"TRAIN set (averaged over k-folds)\")\nplot_confusion(train_y, val_y, label=\"VALIDATION set (built over k-folds)\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most digit are very well predicted on the **validation** set.  \nSome expected confusions are visible:\n- `0` => `6`  \n- `1` => `7`\n- `2` => `7`\n- `4` => `9`\n- `6` => `0`\n- `7` => `9`, `1`\n- `8` => `2`, `4`, `6`\n- `9` => `4`, `5`, `7`\n\nThe `9` and the `1` have the highest error rate, which is rather expected, \nbut the `1` digit is only confused with the `7`. \nHence the `9` is problematic because of the proximity to several other classes.\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Examples of predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_digits(test_X, None, test_sel[0:60], predict=test_y, suptitle=\"Test digits\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Nice!\n\nThe predictions all match the visual digit for 60 random cases."},{"metadata":{},"cell_type":"markdown","source":"# Misclassified digits in the validation set\nBellow are the first 60 missed predictions\n- green = true value\n- purple = predicted"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"miss = np.where(train_y != val_y)[0]\nplot_digits(train_X[miss], train_y[miss], np.arange(60), predict=val_y[miss], \n            suptitle=\"Misclassified validation digits\")\nprint(f\"There are {len(miss)} missclassified digits\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Curiously many of the misclassified digits shown here can be correctly classified by a human. \nWhile the global accuracy is good, this means that the classifier has still room for improvement. It is a bit surprising as the accurcy is already so close to the expected maximum reachable accuracy (0.998) by deep networks."},{"metadata":{},"cell_type":"markdown","source":"# Check predicted digits distribution"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def countplot_ratio(df, x, y, hue, ax=None):\n    # cook book : https://github.com/mwaskom/seaborn/issues/1027\n    df = df[x].groupby(df[hue]) \\\n       .value_counts(normalize=True) \\\n       .rename(y) \\\n       .reset_index()\n    sns.barplot()\n    sns.barplot(x=x, y=y, hue=hue, data=df, ax=ax)\n    \n        \ndf_count= pd.DataFrame({'digit': np.concatenate([train_y, test_y]), \n                        'data set': np.concatenate([ [\"train\"] * len(train_X), [\"test\"] * len(test_X)])})\nfig, ax = plt.subplots(1, 1, figsize=(10,4))\ncountplot_ratio(df_count, 'digit', 'ratio', 'data set' )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Distribution of digits is very close between test and train data sets.  \nThis is a very clean result, indicating that the test split has been built so the digit ratios are the same in both sets.\n"},{"metadata":{},"cell_type":"markdown","source":"# Submit predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test = pd.DataFrame({'ImageId': 1 + np.arange(len(test_X)), \n                        'Label':   test_y })\ndf_test.to_csv('submit.csv', index=False)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusions\n\nThe kaggle MNIST dataset is really easy to go along with due to its nice preprocessing as seen in the exploratory section.  \nOn overall this is a fun dataset to play with :-)\n\nIt behaves very well with convolutional neural networks:\n- 0.991 is reaches easily on test set with no particular adjustments other than using 5-folds of 30 epochs\n- 0.994 is attained with data augmentation by including rotation, shear, zoom and shift\n- 0.995 Next milestone is with a sandwich of Conv2D + BatchNorm + MaxPool  using 16/32/64 convs\n    - I noticed that BatchNorm behaved correctly if coupled with a MaxPool Layer, otherwise accuracy was as low as 0.1\n- 0.9965 A significant improvement comes by doubling the size of the conv layers to 32/64/128\n- 0.9970 is finally reached by pushing epochs to 200 meaning that the nb of epochs is not very significant above 30\n"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}