{"cells":[{"metadata":{},"cell_type":"markdown","source":"# MNIST digit recognizer: CNN, grid search and data augmentation\nI am using the MNIST digit classificaton problem as an exercice to implement some intermediate technics for image processing using Keras, Tensorflow and Scikit Learn. We will: \n* start by **implementing a simple MultiLayer Perceptron** with some default values for hyper parameters \n* improve this model by doing a **grid search over some hyperparameters**\n* move on on a more complex model by implementing a **Convolutional Neural Network**\n* improve the CNN model using again a **grid search**\n* use **data augmentation** to reduce overfit and further improve performance\n\nThis kernel is designed for intermediate users having some knowledge and experience of neural networks and optimization of hyperparameters in Scikit Learn. You can find [here](https://victorzhou.com/blog/keras-neural-network-tutorial/) an excellent introduction to MLP. The same autor has also a very good [introduction](https://victorzhou.com/blog/keras-cnn-tutorial/) to CNN. You can also refer to this Kaggle [kernel](https://www.kaggle.com/anebzt/mnist-with-cnn-in-keras-detailed-explanation) for detailed information on implementing a CNN.\n\nFor hyperparameter optimization, this amazing [kernel](https://www.kaggle.com/cdeotte/how-to-choose-cnn-architecture-mnist) is an excellent demonstration of manual tuning. Thanks also to the authors of this [post](https://machinelearningmastery.com/grid-search-hyperparameters-deep-learning-models-python-keras/) which got me started on using grid search for Keras hyperparameters, and of this [kernel](https://www.kaggle.com/yassineghouzam/introduction-to-cnn-keras-0-997-top-6) for data augmentation.\n\nWe will follow the following steps to run this kernel:\n1. EDA: initial look at the data\n1. Implementing a simple MLP\n1. Optimize hyperparameters and add dropout to reduce overfit\n1. Implementing a Convolutional Neural Network\n1. Optimize the hyperparameters of the CNN thru GridSearch\n1. Use data augmentation to improve performance"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nfrom sklearn.metrics import accuracy_score\nfrom sklearn.model_selection import train_test_split\n\n# keras import\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Conv2D, MaxPooling2D, AveragePooling2D, Flatten\nfrom keras.utils import to_categorical\nfrom keras.callbacks import EarlyStopping\n\n# hyperparameter optimization\nfrom sklearn.model_selection import GridSearchCV\nfrom keras.wrappers.scikit_learn import KerasClassifier\n\n# data augmentation\nfrom keras.preprocessing.image import ImageDataGenerator\n\n# visualisation\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nimport seaborn as sns\n#set figure size\nplt.rcParams['figure.figsize'] = 12, 6\nsns.set_style('white')\n\n# others\nfrom random import randrange\nfrom time import time","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I define here global constants for the number of epochs when training an individual model (n_epochs) or doing a grid search with cross validation (n_epochs_cv). Reduce these 2 values if you want to reduce the running time of the kernel.  \n\nI also define the number of runs for cross validation (when using Scikit Learn GridSearchCV) and the size of the validation set for train_test_split.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_epochs = 30 # 30 \nn_epochs_cv = 10 # 10  # reduce number of epochs for cross validation for performance reason\n\nn_cv = 3\nvalidation_ratio = 0.10\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1. EDA: initial look at the data\nWe load the data and display some sample observations to understand the data structure.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"# load dataset and check dimension\ndata_set = pd.read_csv('/kaggle/input/digit-recognizer/train.csv')\nprint(data_set.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"data_set.sample(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We check the distribution of the 10 classes of digits. They are roughly equivalently represented, therefore we do not need to use stratify when splitting the data set into training and validation sets."},{"metadata":{"trusted":true},"cell_type":"code","source":"# segregate training data set in pixel features and label\ny = data_set['label']\nX = data_set.drop(labels = ['label'], axis=1) \n# free memory\ndel data_set\n\n# check distribution of the handwritten digits\nsns.countplot(y, color='skyblue');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, let's plot a random sample of 60 images to get a *visual feeling* of the classification task. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# show multiple images chosen randomly \nfig, axs = plt.subplots(6, 10, figsize=(10, 6)) # 6 rows of 10 images\n\nfor ax in axs.flat:\n    i = randrange(X.shape[0])\n    ax.imshow(X.loc[i].values.reshape(28, 28), cmap='gray_r')\n    ax.set_axis_off()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Normalize pixel value to range 0 to 1\nX = X / 255.0\n\n# extract train and validation set\nX_train, X_val, y_train, y_val = train_test_split(X, y, test_size = validation_ratio)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. Implementing a simple MLP\nLet's start with a simple multilayer perceptron with only 1 hidden layer:\n* the input layer consists in 128 units, with relu\n* the hidden layer is made of 64 units, with sigmoid activation function\n* the output layer contains one unit per expected class, that is 10 units, and uses a softmax activation function to output probabilities  \n[](http://)\nHow many parameters has this model?\n* the input layer is taking values from 28x28 images: the number of parameters is 28x28 (input size) x 128 (output size) weights + 128 bias values = 100,480\n* for the hidden layer, a similar caculation gives 128 x 64 + 64 = 8,256\n* and for the output layer, we have 64 x 10 + 10 = 650 parameters  \n\nOK, this was easy, but we'll see below that counting parameters of CNN is a bit more tricky. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# define model\nmlp = Sequential()\nmlp.add(Dense(128, activation='relu', input_shape=(784,)))\nmlp.add(Dense(64, activation='sigmoid'))  \nmlp.add(Dense(10, activation='softmax'))\n\nmlp.compile(\n  optimizer='adam',\n  loss='categorical_crossentropy',\n  metrics=['accuracy'],\n)\n\nmlp.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's train the model. I am using a early stop callback to reduce training time.  \nLook at the accuracy on the validation set and compare it with the accuracy on the training set. What can you say about the performance and limits of this simple model?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Train the model\n\n#define callbacks\nearly_stop = EarlyStopping(monitor = 'val_accuracy', mode = 'max', patience=5, restore_best_weights=True)\n\nhistory = mlp.fit(\n    X_train,\n    to_categorical(y_train),\n    epochs = n_epochs,  \n    validation_data = (X_val, to_categorical(y_val)),\n    batch_size = 32,\n    callbacks = [early_stop]\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## What can we say from the results of our training:\n* the performance of our model is given by the best accuracy obtained on the validation set: around 0.975 (depending on the run)\n* the best scores are already obtained after more or less 10 epochs, meaning that there is probably no or little gain to increase number of epochs\n* and there is clear sign of overfitting as the loss for the validation set (around 0.1) is roughly 10 times the loss on the training set (around 0.01).  \n\nLet's plot the accuracy for the training and validation to confirm this last point.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare accuracy accuracy on training and validation data\ndf_history = pd.DataFrame(history.history)\nsns.lineplot(data=df_history[['accuracy','val_accuracy']], palette=\"tab10\", linewidth=2.5);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Optimize hyperparameters and add dropout to reduce overfit\nAs our MLP model shows sign of overfit, we will now add dropout layers to try to fix this. The dropout rate will be determined thru a grid search, along with the batch size.  \nFor this, we will use the KerasClassifier wrapper for Scikit Learn, which gives us a Scikit Learn estimator that we can optimize with GridSearchCV (cf Keras [documentation](https://keras.io/scikit-learn-api/)).  \n\nResults are:\n* no more overfit from dropout rate = 0.2 and above \n* no or minimal improvement of the accuracy on the validation set: **we are reaching the limit of a simple MLP model**\n* noticeable degradation of results (both on training and validation sets) for rate of 0.4 and above "},{"metadata":{"trusted":true},"cell_type":"code","source":"start=time()\n\n# define a function to create model, required for KerasClassifier\n# the function takes drop_out rate as argument so we can optimize it  \ndef create_mlp_model(dropout_rate=0):\n    # create model\n    model = Sequential()\n    model.add(Dense(128, activation='relu', input_shape=(784,))) \n    # add a dropout layer if rate is not null\n    if dropout_rate != 0:\n        model.add(Dropout(rate=dropout_rate))        \n    model.add(Dense(64, activation='sigmoid')) \n    # add a dropout layer if rate is not null    \n    if dropout_rate != 0:\n        model.add(Dropout(rate=dropout_rate))           \n    model.add(Dense(10, activation='softmax'))\n    \n    # Compile model\n    model.compile( \n        optimizer='adam',\n        loss='categorical_crossentropy',\n        metrics=['accuracy'],\n        )    \n    return model\n\n# define function to display the results of the grid search\ndef display_cv_results(search_results):\n    print('Best score = {:.4f} using {}'.format(search_results.best_score_, search_results.best_params_))\n    means = search_results.cv_results_['mean_test_score']\n    stds = search_results.cv_results_['std_test_score']\n    params = search_results.cv_results_['params']\n    for mean, stdev, param in zip(means, stds, params):\n        print('mean test accuracy +/- std = {:.4f} +/- {:.4f} with: {}'.format(mean, stdev, param))    \n    \n# create model\nmodel = KerasClassifier(build_fn=create_mlp_model, verbose=1)\n# define parameters and values for grid search \nparam_grid = {\n    'batch_size': [16, 32, 64],\n    'epochs': [n_epochs_cv],\n    'dropout_rate': [0.0, 0.10, 0.20, 0.30],\n}\ngrid = GridSearchCV(estimator=model, param_grid=param_grid, n_jobs=-1, cv=n_cv)\ngrid_result = grid.fit(X, to_categorical(y))  # fit the full dataset as we are using cross validation \n\n# print out results\nprint('time for grid search = {:.0f} sec'.format(time()-start))\ndisplay_cv_results(grid_result)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The best score (around 0.973) is obtained for a drop out rate of 0.10 and a batch size of 16 or 32 (depending on the run).  \nLooking at the mean accuracy for each set of parameters:\n* batch sizes of 16 or 32 give similar results, but we see a degradation for batch size = 64\n* droprout rate of 0.10 consistently gives better results, regardless of the batch size\n\nOne point that puzled me is that the best score is lower than the accuracy that we obtained before without hyperparameter optimization (!?). My guess is that this best score value is an average of the score of the best estimator over all epochs (I'd appreciate a lot if someone can point me to the exact reason). For this reason, I reload the best estimator and train it on the full dataset. Compared to the first network that we trained, there is a slight improvement of the accuracy at 0.978."},{"metadata":{"trusted":true},"cell_type":"code","source":"# reload best model\nmlp = grid_result.best_estimator_ \n\n# retrain best model on the full training set \nhistory = mlp.fit(\n    X_train,\n    to_categorical(y_train),\n    validation_data = (X_val, to_categorical(y_val)),\n    epochs = n_epochs,\n    callbacks = [early_stop]    \n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get prediction on validation dataset \ny_pred = mlp.predict(X_val)\nprint('Accuracy on validation data = {:.4f}'.format(accuracy_score(y_val, y_pred)))\n\n# plot accuracy on training and validation data\ndf_history = pd.DataFrame(history.history)\nsns.lineplot(data=df_history[['accuracy','val_accuracy']], palette=\"tab10\", linewidth=2.5);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load test data and make prediction\ntest = pd.read_csv('/kaggle/input/digit-recognizer/test.csv')\ny_test = mlp.predict(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# convert prediction to df\nsubmission = pd.DataFrame(data=y_test)\n\n# set label as the 0-9 class with highest value \nsubmission['Label'] = submission.idxmax(axis=1)\nsubmission['ImageId'] = np.asarray([i+1 for i in range(submission.shape[0])])\n\nsubmission.to_csv('submission-mlp_dropout.csv', \n                  columns=['ImageId','Label'],\n                  header=True,\n                  index=False)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4. Implementing a Convolutional Neural Network\nWe implement now a CNN since MLP model is limited at 0.975 accuracy. I selected the architecture thru trials based on several examples. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Reshape the images\nimg_size = 28\nX_cnn = X.values.reshape(-1, img_size, img_size, 1)\n# check \nprint(X_cnn.shape)\n\nX_train, X_val, y_train, y_val = train_test_split(X_cnn, y, test_size = validation_ratio)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# function to create the model for Keras wrapper to scikit learn\n# we will optimize the type of pooling layer (max or average) and the activation function of the 2nd and 3rd convolution layers \ndef create_cnn_model(pool_type='max', conv_activation='sigmoid', dropout_rate=0.10):\n    # create model\n    model = Sequential()\n    \n    # first layer: convolution\n    model.add(Conv2D(16, kernel_size=(5, 5), activation='relu', input_shape=(28, 28, 1))) \n        \n    # second series of layers: convolution, pooling, and dropout\n    model.add(Conv2D(32, kernel_size=(5, 5), activation=conv_activation))  \n    if pool_type == 'max':\n        model.add(MaxPooling2D(pool_size=(2, 2)))\n    if pool_type == 'average':\n        model.add(AveragePooling2D(pool_size=(2, 2)))\n    if dropout_rate != 0:\n        model.add(Dropout(rate=dropout_rate))     \n    \n    # third series of layers: convolution, pooling, and dropout    \n    model.add(Conv2D(64, kernel_size=(3, 3), activation=conv_activation))   # 32   \n    if pool_type == 'max':\n        model.add(MaxPooling2D(pool_size=(2, 2)))\n    if pool_type == 'average':\n        model.add(AveragePooling2D(pool_size=(2, 2)))\n    if dropout_rate != 0:\n        model.add(Dropout(rate=dropout_rate))     \n      \n    # fourth series\n    model.add(Flatten())         \n    model.add(Dense(64, activation='sigmoid')) # 64\n    # add a dropout layer if rate is not null    \n    if dropout_rate != 0:\n        model.add(Dropout(rate=dropout_rate)) \n        \n    model.add(Dense(10, activation='softmax'))\n    \n    # Compile model\n    model.compile( \n        optimizer='adam',\n        loss='categorical_crossentropy',\n        metrics=['accuracy'],\n        )    \n    return model\n\ncnn = create_cnn_model()\n\ncnn.compile(\n  optimizer='adam',\n  loss='categorical_crossentropy',  \n  metrics=['accuracy'],\n)\n\ncnn.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check the number of parameters:\n* the first layer is a convolution layer with 16 kernels of size 5x5. Each kernel has 26 parameters (25 weigths plus bias). Total is 16 x 26 = 416\n* the second layer is also a convolution layer with 32 kernels of size 5x5. There are 16 input images from the output of the previous layer. For every kernel, the 16 images are combined together, with 16 (input size) x 25 (kernel size) weights plus bias (that is 401 parameters for each kernel). Total number of parameters for this layers is 401 x 32 = 12,832\n* using same calculation logic, the number of parameters for the 3rd layer is ( 32 (input size) x 9 (kernel size) + 1 ) x 64 (number of kernels) = 18,4896\n* for the dense layer, its input is the result of the flatten layer (taking as input the 64 x 16 images from the pooling layer, and mapping them to a flat array of 64 x 16 = 1024). Output size is 64, therefore number of parameters is 64 x (1024 weights + bias) = 65,600\n* same way, the final layer has (64 x 1) x 10 = 650 parameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Train the default CNN model\nhistory = cnn.fit(\n    X_train,\n    to_categorical(y_train),\n    epochs=n_epochs,  \n    validation_data=(X_val, to_categorical(y_val)), \n    batch_size=32,\n    callbacks = [early_stop]\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We gained about 1% of accuracy with the CNN! But there is a cost: training the CNN model takes 45 seconds per epoch, about 10 times what the MLP  required!  \nLet's optimize now some hyperparameters (pooling type and activation function) to further improve the performance of this CNN. "},{"metadata":{},"cell_type":"markdown","source":"## 5. Optimize the hyperparameters of the CNN thru GridSearch"},{"metadata":{"trusted":true},"cell_type":"code","source":"# optimize model \nstart = time()\n\n# create model\nmodel = KerasClassifier(build_fn=create_cnn_model, verbose=1)\n# define parameters and values for grid search \nparam_grid = {\n    'pool_type': ['max', 'average'],\n    'conv_activation': ['sigmoid', 'tanh'],    \n    'epochs': [n_epochs_cv],\n}\n\ngrid = GridSearchCV(estimator=model, param_grid=param_grid, n_jobs=-1, cv=n_cv)\ngrid_result = grid.fit(X_train, to_categorical(y_train))\n\n# summarize results\nprint('time for grid search = {:.0f} sec'.format(time()-start))\ndisplay_cv_results(grid_result)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 6. Use data augmentation to improve performance"},{"metadata":{"trusted":true},"cell_type":"code","source":"# optimize parameters of the fit method \ncnn_model = create_cnn_model(pool_type = grid_result.best_params_['pool_type'],\n                             conv_activation = grid_result.best_params_['conv_activation'])\n\n# With data augmentation \ndatagen = ImageDataGenerator(\n        rotation_range=10,  # randomly rotate images in the range (degrees, 0 to 180)\n        zoom_range = 0.1, # Randomly zoom image \n        width_shift_range=0.1,  # randomly shift images horizontally (fraction of total width)\n        height_shift_range=0.1,  # randomly shift images vertically (fraction of total height)\n        horizontal_flip=False,  # randomly flip images\n        vertical_flip=False,  # randomly flip images\n        fill_mode='constant', cval = 0.0)\n\ndatagen.fit(X_train)\n\nhistory = cnn_model.fit_generator(datagen.flow(X_train,to_categorical(y_train), batch_size=32),\n                                  epochs = n_epochs, \n                                  validation_data = (X_val,to_categorical(y_val)),\n                                  verbose = 1, \n                                  steps_per_epoch = X_train.shape[0] / 32,\n                                  callbacks = [early_stop])\n\n# plot accuracy on training and validation data\ndf_history = pd.DataFrame(history.history)\nsns.lineplot(data=df_history[['accuracy','val_accuracy']], palette=\"tab10\", linewidth=2.5);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The accuracy on the validation is consistently higher than on the training set. The model is too much constrained: let's retrain it with data augmentation but without dropout and see how it performs. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# optimize parameters of the fit method \ncnn_model = create_cnn_model(pool_type = grid_result.best_params_['pool_type'],\n                             conv_activation = grid_result.best_params_['conv_activation'], \n                            dropout_rate=0.0)\n\n#define early stop on the accuracy as this is the metric we want to improve\nearly_stop = EarlyStopping(monitor = 'accuracy', mode = 'max', patience=5, restore_best_weights=True)\nhistory = cnn_model.fit_generator(datagen.flow(X_train,to_categorical(y_train), batch_size=32),\n                                  epochs = n_epochs, \n                                  validation_data = (X_val,to_categorical(y_val)),\n                                  verbose = 1, \n                                  steps_per_epoch = X_train.shape[0] / 32,\n                                  callbacks = [early_stop])\n\n# plot accuracy on training and validation data\ndf_history = pd.DataFrame(history.history)\nsns.lineplot(data=df_history[['accuracy','val_accuracy']], palette=\"tab10\", linewidth=2.5);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# save weights\ncnn_model.save_weights('mnist_cnn.h5')\ntest = pd.read_csv('/kaggle/input/digit-recognizer/test.csv')\n\nX_test = test.values.reshape(-1, img_size, img_size, 1)\ny_test = cnn_model.predict(X_test)\n\n# convert to df\nsubmission = pd.DataFrame(data=y_test)\n\n# set label as the 0-9 class with highest value \nsubmission['Label'] = submission.idxmax(axis=1)\nsubmission['ImageId'] = np.asarray([i+1 for i in range(submission.shape[0])])\n\nsubmission.to_csv('submission-cnn.csv', \n                  columns=['ImageId','Label'],\n                  header=True,\n                  index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}