<jupyter_start><jupyter_text>Bayesian Optimization of CNN with Scikit-OptimizeIn this notebook, we will use **Bayesian Optimization** to select the best **hyperparameters** for a CNN that recognizes digits in images, using the MNIST dataset and the open source Python package [Scikit-Optimize](https://scikit-optimize.github.io/stable/index.html).The MNIST dataset is availale in [Kaggle](https://www.kaggle.com/c/digit-recognizer/data). Download dataset- Navigate to the [MNIST website in Kaggle](https://www.kaggle.com/c/digit-recognizer/data)- Download the train.csv file- Unzip and copy the train.csv file to where you see the SAVE_DATASETS-HERE.txt file- Rename to mnist.csv**Remember that you need to be logged in to be able to download the dataset** Notebook content- Data Preparation- Set up a simple CNN- Set up the hyperparameter search shape- Set up the objective function- Perform Bayesian Optimization- Evaluate Model Performance<jupyter_code>import itertools

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from keras.utils.np_utils import to_categorical
from keras.models import Sequential, load_model
from keras.layers import Dense, Flatten, Conv2D, MaxPool2D
from keras.optimizers import Adam
from keras.callbacks import ReduceLROnPlateau
from skopt import gp_minimize
from skopt.space import Real, Categorical, Integer
from skopt.plots import plot_convergence
from skopt.plots import plot_objective, plot_evaluations
from skopt.utils import use_named_args<jupyter_output><empty_output><jupyter_text>Data PreparationThe dataset contains information about images, each image is a hand-written digit. The aim is to have the computer predict which digit was written by the person, automatically, by "looking" at the image. Each image is 28 pixels in height and 28 pixels in width (28 x 28), making a total of 784 pixels. Each pixel value is an integer between 0 and 255, indicating the darkness in a gray-scale of that pixel.The data is stored in a dataframe where each each pixel is a column (so it is flattened and not in the 28 x 28 format). The data set the has 785 columns. The first column, called "label", is the digit that was drawn by the user. The rest of the columns contain the pixel-values of the associated image.<jupyter_code># Load the data

data = pd.read_csv("../mnist.csv")

# first column is the target, the rest of the columns
# are the pixels of the image

# each row is 1 image
data.head()
# split dataset into a train and test set

X_train, X_test, y_train, y_test = train_test_split(
    data.drop(['label'], axis=1), # the images
    data['label'], # the target
    test_size = 0.1,
    random_state=0)

X_train.shape, X_test.shape
# number of images for each digit

g = sns.countplot(x=y_train)
plt.xlabel('Digits')
plt.ylabel('Number of images')<jupyter_output><empty_output><jupyter_text>There are roughly the same amount of images for each of the 10 digits.  Image re-scalingWe re-scale data for the CNN, between 0 and 1.<jupyter_code># Re-scale the data

# 255 is the maximum value a pixel can take

X_train = X_train / 255
X_test = X_test / 255<jupyter_output><empty_output><jupyter_text>ReshapeThe images were stored in a pandas dataframe as 1-D vectors of 784 values. For a CNN with Keras, we need tensors with the following dimensions: width x height x channel. Thus, we reshape all data to 28 x 2 8 x 1, 3-D matrices. The 3rd dimension corresponds to the channel. RGB images have 3 channels. MNIST images are in gray-scale, thus they have only one channel in the 3rd dimension.<jupyter_code># Reshape image in 3 dimensions:
# height: 28px X width: 28px X channel: 1 

X_train = X_train.values.reshape(-1,28,28,1)
X_test = X_test.values.reshape(-1,28,28,1)<jupyter_output><empty_output><jupyter_text>Target encoding<jupyter_code># the target is 1 variable with the 9 different digits
# as values

y_train.unique()
# For Keras, we need to create 10 dummy variables,
# one for each digit

# Encode labels to one hot vectors (ex : digit 2 -> [0,0,1,0,0,0,0,0,0,0])

y_train = to_categorical(y_train, num_classes = 10)
y_test = to_categorical(y_test, num_classes = 10)

# the new target
y_train<jupyter_output><empty_output><jupyter_text>Let's print some example images.<jupyter_code># Some image examples 

g = plt.imshow(X_train[0][:,:,0])
# Some image examples 

g = plt.imshow(X_train[10][:,:,0])<jupyter_output><empty_output><jupyter_text>Define the CNNWe will create a CNN, with 2 Convolutional layers followed by Pooling, and varying number of fully-connected Dense layers.In fact, the number of fully-connected Dense layers and the number of Neurons in each one of the Dense layers, are some of the hyperparameters we want to optimize.We could also optimize the number of Convolutional layers. But we will keep that for later, and here we keep it a bit simpler.<jupyter_code># function to create the CNN

def create_cnn(
    learning_rate,
    num_dense_layers,
    num_dense_nodes,
    activation,
):
    """
    Hyper-parameters:
    learning_rate:     Learning-rate for the optimizer.
    num_dense_layers:  Number of dense layers.
    num_dense_nodes:   Number of nodes in each dense layer.
    activation:        Activation function for all layers.
    """

    # Start construction of a Keras Sequential model.
    model = Sequential()

    # First convolutional layer.
    # There are many hyper-parameters in this layer
    # For this demo, we will optimize the activation function only.
    model.add(Conv2D(kernel_size=5, strides=1, filters=16, padding='same',
                     activation=activation, name='layer_conv1'))
    model.add(MaxPool2D(pool_size=2, strides=2))

    # Second convolutional layer.
    # Again, we will only optimize the activation function.
    model.add(Conv2D(kernel_size=5, strides=1, filters=36, padding='same',
                     activation=activation, name='layer_conv2'))
    model.add(MaxPool2D(pool_size=2, strides=2))

    # Flatten the 4-rank output of the convolutional layers
    # to 2-rank that can be input to a fully-connected Dense layer.
    model.add(Flatten())

    # Add fully-connected Dense layers.
    # The number of layers is a hyper-parameter we want to optimize.
    # We add the different number of layers in the following loop:
    
    for i in range(num_dense_layers):
        
        # Add the dense fully-connected layer to the model.
        # This has two hyper-parameters we want to optimize:
        # The number of nodes (neurons) and the activation function.
        model.add(Dense(num_dense_nodes,
                        activation=activation,
                        ))

    # Last fully-connected dense layer with softmax-activation
    # for use in classification.
    model.add(Dense(10, activation='softmax'))

    # Use the Adam method for training the network.
    # We want to find the best learning-rate for the Adam method.
    optimizer = Adam(lr=learning_rate)

    # In Keras we need to compile the model so it can be trained.
    model.compile(optimizer=optimizer,
                  loss='categorical_crossentropy',
                  metrics=['accuracy'])

    return model<jupyter_output><empty_output><jupyter_text>Define the Hyperparameter SpaceScikit-optimize provides an utility function to create the range of values to examine for each hyperparameters. More details in [skopt.Space](https://scikit-optimize.github.io/stable/modules/generated/skopt.Space.html)We want to find the following hyper-parameters:- The learning rate of the optimizer.- The number of fully-connected Dense layers.- The number of nodes (neurons) for each of the dense layers.- Whether to use 'sigmoid' or 'relu' activation in all the layers.<jupyter_code>dim_learning_rate = Real(
    low=1e-6, high=1e-2, prior='log-uniform', name='learning_rate',
)

dim_num_dense_layers = Integer(low=1, high=5, name='num_dense_layers')

dim_num_dense_nodes = Integer(low=5, high=512, name='num_dense_nodes')


dim_activation = Categorical(
    categories=['relu', 'sigmoid'], name='activation',
)
# the hyperparameter space grid

param_grid = [dim_learning_rate,
              dim_num_dense_layers,
              dim_num_dense_nodes,
              dim_activation]<jupyter_output><empty_output><jupyter_text>Define the Objective Function<jupyter_code># we will save the model with this name
path_best_model = 'cnn_model.h5'

# starting point for the optimization
best_accuracy = 0
@use_named_args(param_grid)
def objective(
    learning_rate,
    num_dense_layers,
    num_dense_nodes,
    activation,
):
    
    """
    Hyper-parameters:
    learning_rate:     Learning-rate for the optimizer.
    num_dense_layers:  Number of dense layers.
    num_dense_nodes:   Number of nodes in each dense layer.
    activation:        Activation function for all layers.
    """

    # Print the hyper-parameters.
    print('learning rate: {0:.1e}'.format(learning_rate))
    print('num_dense_layers:', num_dense_layers)
    print('num_dense_nodes:', num_dense_nodes)
    print('activation:', activation)
    print()
    
    # Create the neural network with the hyper-parameters.
    # We call the function we created previously.
    model = create_cnn(learning_rate=learning_rate,
                       num_dense_layers=num_dense_layers,
                       num_dense_nodes=num_dense_nodes,
                       activation=activation)

   
    # Set a learning rate annealer
    # this reduces the learning rate if learning does not improve
    # for a certain number of epochs
    learning_rate_reduction = ReduceLROnPlateau(monitor='val_accuracy', 
                                                patience=2, 
                                                verbose=1, 
                                                factor=0.5, 
                                                min_lr=0.00001)
   
    # train the model
    # we use 3 epochs to be able to run the notebook in a "reasonable"
    # time. If we increase the epochs, we will have better performance
    # this could be another parameter to optimize in fact.
    history = model.fit(x=X_train,
                        y=y_train,
                        epochs=3,
                        batch_size=128,
                        validation_split=0.1,
                        callbacks=learning_rate_reduction)

    # Get the classification accuracy on the validation-set
    # after the last training-epoch.
    accuracy = history.history['val_accuracy'][-1]

    # Print the classification accuracy.
    print()
    print("Accuracy: {0:.2%}".format(accuracy))
    print()

    # Save the model if it improves on the best-found performance.
    # We use the global keyword so we update the variable outside
    # of this function.
    global best_accuracy

    # If the classification accuracy of the saved model is improved ...
    if accuracy > best_accuracy:
        # Save the new model to harddisk.
        # Training CNNs is costly, so we want to avoid having to re-train
        # the network with the best found parameters. We save it instead
        # as we search for the best hyperparam space.
        model.save(path_best_model)
        
        # Update the classification accuracy.
        best_accuracy = accuracy

    # Delete the Keras model with these hyper-parameters from memory.
    del model

    
    # Remember that Scikit-optimize always minimizes the objective
    # function, so we need to negate the accuracy (because we want
    # the maximum accuracy)
    return -accuracy<jupyter_output><empty_output><jupyter_text>Test run<jupyter_code># Before we run the hyper-parameter optimization, 
# let's first check that the everything is working
# by passing some default hyper-parameters.

default_parameters = [1e-5, 1, 16, 'relu']

objective(x=default_parameters)<jupyter_output>learning rate: 1.0e-05
num_dense_layers: 1
num_dense_nodes: 16
activation: relu

Epoch 1/3
266/266 [==============================] - 33s 121ms/step - loss: 2.2890 - accuracy: 0.1089 - val_loss: 2.2262 - val_accuracy: 0.2344
Epoch 2/3
266/266 [==============================] - 30s 113ms/step - loss: 2.1948 - accuracy: 0.2784 - val_loss: 2.0903 - val_accuracy: 0.4019
Epoch 3/3
266/266 [==============================] - 33s 122ms/step - loss: 2.0400 - accuracy: 0.4242 - val_loss: 1.8993 - val_accuracy: 0.4860

Accuracy: 48.60%<jupyter_text>We obtained a mediocre accuracy, but all our code is working. So let's get started with the Optimization now!!  Bayesian Optimization<jupyter_code># gp_minimize performs by default GP Optimization 
# using a Marten Kernel

gp_ = gp_minimize(
    objective, # the objective function to minimize
    param_grid, # the hyperparameter space
    x0=default_parameters, # the initial parameters to test
    acq_func='EI', # the acquisition function
    n_calls=30, # the number of subsequent evaluations of f(x)
    random_state=0, 
)<jupyter_output>learning rate: 1.0e-05
num_dense_layers: 1
num_dense_nodes: 16
activation: relu

Epoch 1/3
266/266 [==============================] - 30s 109ms/step - loss: 2.2821 - accuracy: 0.1131 - val_loss: 2.2096 - val_accuracy: 0.2550
Epoch 2/3
266/266 [==============================] - 29s 109ms/step - loss: 2.1760 - accuracy: 0.2998 - val_loss: 2.0578 - val_accuracy: 0.4373
Epoch 3/3
266/266 [==============================] - 29s 109ms/step - loss: 2.0070 - accuracy: 0.4580 - val_loss: 1.8437 - val_accuracy: 0.5341

Accuracy: 53.41%

learning rate: 2.4e-04
num_dense_layers: 4
num_dense_nodes: 440
activation: sigmoid

Epoch 1/3
266/266 [==============================] - 36s 134ms/step - loss: 2.3177 - accuracy: 0.1025 - val_loss: 2.3125 - val_accuracy: 0.1146
Epoch 2/3
266/266 [==============================] - 37s 138ms/step - loss: 2.0656 - accuracy: 0.2412 - val_loss: 0.7723 - val_accuracy: 0.7500
Epoch 3/3
266/266 [==============================] - 35s 130ms/step - loss: 0.6654 - accuracy: [...]<jupyter_text>Analyze results<jupyter_code># function value at the minimum.
# note that it is the negative of the accuracy

"Best score=%.4f" % gp_.fun
gp_.x
gp_.space
print("""Best parameters:
=========================
- learning rate=%.6f
- num_dense_layers=%d
- num_nodes=%d
- activation = %s""" % (gp_.x[0], 
                gp_.x[1],
                gp_.x[2],
                gp_.x[3]))<jupyter_output>Best parameters:
=========================
- learning rate=0.001810
- num_dense_layers=1
- num_nodes=303
- activation = relu<jupyter_text>Convergence<jupyter_code>plot_convergence(gp_)<jupyter_output><empty_output><jupyter_text>Partially dependency plots[plot_objective](https://scikit-optimize.github.io/stable/modules/generated/skopt.plots.plot_objective.htmlskopt.plots.plot_objective)<jupyter_code>dim_names = ['learning_rate', 'num_dense_nodes', 'num_dense_layers']
plot_objective(result=gp_, plot_dims=dim_names)
plt.show()<jupyter_output><empty_output><jupyter_text>Evaluation order[plot_evaluations](https://scikit-optimize.github.io/stable/modules/generated/skopt.plots.plot_evaluations.html)<jupyter_code>plot_evaluations(result=gp_, plot_dims=dim_names)
plt.show()<jupyter_output><empty_output><jupyter_text>Evaluate the model<jupyter_code># load best model

model = load_model(path_best_model)
# make predictions in test set

result = model.evaluate(x=X_test,
                        y=y_test)
# print evaluation metrics

for name, value in zip(model.metrics_names, result):
    print(name, value)<jupyter_output>loss 0.04718650504946709
accuracy 0.9833333492279053<jupyter_text>Confusion matrix<jupyter_code># Predict the values from the validation dataset
y_pred = model.predict(X_test)

# Convert predictions classes to one hot vectors 
y_pred_classes = np.argmax(y_pred, axis = 1)

# Convert validation observations to one hot vectors
y_true = np.argmax(y_test, axis = 1)

# compute the confusion matrix
cm = confusion_matrix(y_true, y_pred_classes) 

cm
# let's make it more colourful
classes = 10

plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
plt.title('Confusion matrix')
plt.colorbar()
tick_marks = np.arange(classes)
plt.xticks(tick_marks, range(classes), rotation=45)
plt.yticks(tick_marks, range(classes))

for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
    plt.text(j, i, cm[i, j],
             horizontalalignment="center",
             color="white" if cm[i, j] > 100 else "black",
            )

plt.tight_layout()
plt.ylabel('True label')
plt.xlabel('Predicted label')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10435620895342532, 'fraction_numerical': 0.03241549677652588, 'mean_word_length': 3.6807670614777215, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 68, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '159327', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5544, 'n_tokens_neox': 4988, 'n_words': 1637, 'path': 'Section-10-Scikit-Optimize/10-Bayesian-Optimization-CNN.ipynb', 'repo_name': 'ankitario/hyperparameter-optimization', 'size': '16597'}