Download and load Raman spectra data

You can download data from:

Bocklitz-Lab. Example-Raman-spectral-analysis. https://github.com/Bocklitz-Lab/Example-Raman-spectral-analysis. 2021. GitHub

The file DATA_dp_wc_bc.csv contains Raman spectra from mice tissues, averaged from each Raman map after despiking, wavenumber calibration and baseline correction. The ID, annotation, and type of the tissues are included as separate columns in the file.
Load Data: df = pd.read_csv("data/DATA_dp_wc_bc.csv")

Divide the dataset into training and testing.

It usually happens that a data set contains multiple samples from a single organism, which could lead to errors and the concept of testing with training.
Therefore, to avoid bias, samples from the same mouse must be in the same subset either CV or T_V_T.

• T_V_T: training , validation, testing sets.
• CV: cross-validation with training and testing folds.
• Samples from the same subset start with a similar name, i.e., M10
• We select the first component of the name and transform them into a categorical class to split the dataset.

Code

# insert in postion 2 new column containing only the string before the first "_"
df.insert(2, 'name_split', df['Name'].str.split('_').str[:1].str.join(""))

# insert in postion 3 new column transfoming name_split to categorical 
df.insert(3, 'name_split_cat', df['name_split'].astype('category').cat.codes)

# Find number of independet mice
number_samples = np.max(df['name_split_cat'])

print("number of independent samples : %d" %number_samples)

# print a section of the dataset, name_split_cat show the mouse index
print(df.iloc[15:25, 1:4])

Result

number of independent samples : 143 

    Name      name_split  name_split_cat
M10_R1_I_007     M10            0      
M10_R1_I_007     M10            0      
M10_R1_I_008     M10            0      
M10_R1_I_008     M10            0      
M10_R1_II_001    M10            0      
M101_B2_I_001   M101            1      
M106_B2_II_001  M106            2      
M11_C2_I_001     M11            3      
M11_C2_I_002     M11            3      
M11_C2_I_003     M11            3  

We can observe some data examples. The samples from mice M10, M101, M106, and M11 have now a categorical names 0, 1, 2, and 3 respectively.
It will be used later to split the dataset on the mouse level.

T_V_T

A dataset can be divide into training, testing and validation.
In this example we randomly select training (70%), validation (15%) and testing (15%)
Code

# Shuffle mice index, as example we do not shuffle to visualize the result
# idx = np.random.permutation(np.arange(number_samples))
idx = np.arange(number_samples)

# select first 70% data for training
train_index = idx[:int(number_samples*.7)]

# select last 30% data to split into validation and test
temporal_index = idx[int(number_samples*.7):]

# select 15% data for validation
validation_index = temporal_index[:int(number_samples*.15)]
# select 15% data for testing
test_index = temporal_index[int(number_samples*.15):]

print("train idex : ", train_index)
print("validation idext : ", validation_index)
print("test index : ", test_index)

Result

train idex :  [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99]
validation idext :  [ 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120]
test index :  [ 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142]

K-fold Cross validation

We can use sklearn to select the index to train and test K different models on different data folds.
In this example, we use five splits.
Note:

• A new model must be generated for each fold. Retraining a single model should be avoided.
• Auxiliar functions must be defined outside the cross-validation for the loop.
• In the following, we compute using the last fold's `train_index` and `test_index` to find the learning rate that best suit the model.
• We run cross-validation later.

Code

# import sklearn module
from sklearn.model_selection import KFold

# functions definitions to train a model here!

# define X as an array containing all mice index (0-142)
X = np.arange(number_samples)

# select the number of folds
kf = KFold(n_splits=5, shuffle=True)

# generate the splits
kf.get_n_splits(X)

# get the index for each split
for fold, (train_index, test_index) in enumerate(kf.split(X)):
    
    # model definition here!
    # model.fit here!
    
    # print index for testing, the index for training are the remainder
    print("Fold ", fold ," - Test index:", test_index)

Result

Fold  0  - Test index: [  1  17  33  36  38  39  56  62  64  67  81  82  83  86  88  95  97  98 100 105 106 111 113 119 120 122 124 129 136]
Fold  1  - Test index: [  5  12  15  16  20  21  22  31  32  34  37  41  45  51  57  59  61  63  69  73  89  91  93  94 115 128 134 137 138]
Fold  2  - Test index: [  8   9  11  14  24  26  27  49  50  60  65  74  76  78  79  87  92 101 102 103 104 108 109 114 117 126 130 139 140]
Fold  3  - Test index: [  0   3   4   7  13  19  35  42  46  47  53  54  55  58  68  72  75  85 99 107 110 112 121 123 125 127 141 142]
Fold  4  - Test index: [  2   6  10  18  23  25  28  29  30  40  43  44  48  52  66  70  71  77 80  84  90  96 116 118 131 132 133 135]

Examine the distribution of the classes

We can examine the distribution of the classes of the training subset by ploting the histogram over the classes.
Code

df[['class']].hist(figsize = (10, 5))
plt.title("Distribution of the classes")
plt.xticks(np.arange(4), classes)
plt.show()

Result

The dataset has four classes : Normal, HP, Adenom and Karzinom.
It is unbalanced. Nevertheless, it is not a critical situation, and we can continue.

Split dataset using the fold index

Prepare data to train: Tensorflow requieres X and Y sets.
X (array): samples x dimension, we need to remove extra columns.
Y (array): samples x label.
Code

# list the dataset keys (column names)
keys_df = list(df.keys())

# elements to remove
remove_elements = ["Unnamed: 0", "Name", 'Annotation', 'Tissue', 'name_split', 'name_split_cat', 'class']

# remove keys from dataset keys list
for element in remove_elements:
    keys_df.remove(element)
    
# divide data into train and validation    
Y_train = df.loc[df["name_split_cat"].isin(train_index), "class"]
X_train = df.loc[df["name_split_cat"].isin(train_index), keys_df]
    
Y_validate = df.loc[df["name_split_cat"].isin(test_index), "class"]
X_validate = df.loc[df["name_split_cat"].isin(test_index), keys_df]
        
# transform data into numpy arrays, cast Y to numpy int 32 dytpe    
Y_train = np.array(Y_train, dtype=np.int32)
X_train = np.array(X_train)
    
Y_validate = np.array(Y_validate, dtype=np.int32)
X_validate = np.array(X_validate)

Normalization methods

We define three normalization methods.

• Normalize spectra by the norm.
• Standard normal variate.
• Normalize spectra by the integral between a range of wavenumbers .

def vector_normalization(v):
    # normalize spectra by the L2-norm
    norm = np.linalg.norm(v)
    if norm == 0: 
        return v
    return v / norm


def snv(v):
    # SNV: Standard normal variate
    if  np.std(v) ==0:
        return v
    return (v - np.mean(v))/ np.std(v)
    

def normalize_integrated(v, wavenumber, wavenumber_min=1630, wavenumber_max=1690):
    # normalize spectra by the integral between the interval Wavenumbers/cm^{-1} wavenumber_min-wavenumber_max 
    valid_Wavenumber = np.multiply(np.array(wavenumber)>=wavenumber_min, np.array(wavenumber)<=wavenumber_max).astype(float)
    integral = np.sum(valid_Wavenumber*v)

    if integral == 0: 
        return v
    return v / integra

We plot the original Raman spectra and three normalization methods:

Ploting training and testing data

We can generate a plot of the mean and standard deviation of the training and testing dataset to observe similarities. The code for the plot can be found in the Jupyter notebook here .

Model definition

We define a simple model, the input has a shape of 691 Raman numbers, and the output has four neurons corresponding to the number of classes. We include three convolutional layers with bias and activations and two dense layers.

model = tf.keras.Sequential([
keras.Input(shape=(696, 1)), 
keras.layers.Conv1D(32, 5, padding="same", use_bias=True, activation='relu'),
keras.layers.Conv1D(64, 5, padding="same", use_bias=True, activation='relu'),
keras.layers.Conv1D(128, 5, padding="same", use_bias=True, activation='relu'),
keras.layers.Flatten(),
keras.layers.Dense(256, activation='relu'),
keras.layers.Dense(4, activation='softmax')])

Neither the network configuration nor the final classification is our primary interest in this practice. Our goal is to recognize when the learning rate is adequate or not.

Define a number of epochs

There is no rule to define the number of necessary epochs. Some conditions require more epochs:

• CNN with several layers to train.
• A small learning rate.
• A large dataset.

We can start with a lengthy number of epochs epochs = 1000 to see how the training parameters respond (loss and accuracy).

Learning rate too large, learning_rate = 0.01

What can we observe on the plots?

• The training Accuracy (TA) value is set to 0 before the first iteration. After the first epoch, TA has a value of 0.52. This value is obtained after randomly initializing the kernels. That value is not relevant.
• Balance Accuracy shows a wide variation between two values, 0.5 and 0.7. However, the behavior of the model is better described with the loss and the accuracy.
• The accuracies and loses maintain invariable values. Thus, the model does not learn.
• The learning rate is too high that it can not converge.
• When the learning rate is very low, the same phenomena are observed.

Learning rate large, learning_rate = 0.001

What can we observe on the plots?

• The model reach high training accuracy and low loss (1e-5) after 300 epochs.
• Around 300 epochs, the model reaches 100% accuracy, and the loss is too small and reduces over time. The model ceases to find a general solution, fits excessively to the input data (overfitting), and becomes more restrictive. .
• Overfitting is the most common reason why the loss for testing increases.
• The learning rate is high, it produces large fluctations in the training accuracy in the first epochs.
• It is possible to fall into local minima, and the model can get stuck with a high learning rate.
• Sometimes, the optimizer gets an good solution, but it jupms.
• The validation loss increases indicating overfitting.
• There was no need to train 500 epochs. Indeed it is better to stop the training earlier.
• In the case of cross-validation, it is recommended to test a fold to select the number of epochs.
• The results for testing are:
  loss: 9.013
Accuracy: 0.572
Balance Accuracy: 0.365

Learning rate too small, learning_rate = 1e-6

What can we observe on the plots?

• The training accuracy (TA) value is set to 0 before the first iteration. After the first epoch, TA has a value of 0.52. This value is obtained after randomly initializing the kernels. That value is not relevant.
• The Balance accuracy curve oscillates a lot, the model behavior is better observed using Accuracy and loss.
• The accuracies maintain invariable values. However the losses decrease.
• It needs a lot of time to learn, the loss decreases very slowly.

Learning rate small, learning_rate = 1e-4

What can we observe on the plots?

• The model reaches training accuracy = 100% and loss = 0.0011 after 500 iterations.
• It requires longer time that with a larger learning rate and the loss is bigger.
• Test loss increases.
• The result of testing varies in a small range.
• On two occasions, the model fell into the wrong minima and had an inadequate result.
• The results are:
  loss: 6.303
Accuracy: 0.559
Balance Accuracy: 0.456

Learning rate OK, learning_rate = 1e-3

What can we observe on the plots?

• The model reaches training accuracy = 100% and loss = 1e-6 after 300 iterations.
• It requires less time to converge, around 150 epochs.
• Test loss increases.
• In one occasions, the model fell into the wrong minima and had an inadequate result.
• The model learns for the training data, which means the learning rate is adequate.
• Around epoch 250, the model fell into the wrong minima and had an inadequate result; fortunately, it could escape.
• The results are:
  loss: 9.819
Accuracy: 0.591
Balance Accuracy: 0.420

Learning rate scheduler

We can define a function that decreases a initial learning rate during the training process.

Exponential scheduler

This function keeps the initial learning rate for the first N epochs and decreases it exponentially.
We can fix the learning rate until EPOCH_N = 10, and Initial learning_rate = 0.01.
Our network already has a low loss at this value, and then we decrease the learning rate. The exponent of "e" controls the rate of change.

def scheduler(epoch, lr):
    EPOCH_N = 10
    if epoch < EPOCH_N: 
        return lr
    else:
        return lr * tf.math.exp(-0.01)

What can we observe on the plots?

• The model reaches a high accuracy after 200 epochs.
• The minimum loss is loss=0.01, much higher than in previous experiments.
• The fluctuations for testing accuracy and loss are reduced.
• We can observe how the learning rate decreases (green line).
• The fluctuations in testing accuracy and loss reduce when the learning rate decreases.
• The learning rate become too low that there are no variations soon after 200 epochs.
• A large LR at the beginning speeds up the training, and lowering the value prevents the learning from fluctuating.
• Testing loss increases, and testing accuracy does not improve.
• The final result for testing:
  loss: 3.429
Accuracy: 0.578
Balance Accuracy: 0.419

Ladder scheduler

This function drops the learning rate, every epochs_drop epochs.
We can set the epochs_drop=20, drop_rate=0.9 and Initial learning_rate = 0.01. The drop rate controls the rate of change.

def lr_step_decay(epoch, lr):
    drop_rate = 0.9
    epochs_drop = 20.0
    return initial_learning_rate * math.pow(drop_rate, math.floor(epoch/epochs_drop))

What can we observe on the plots?

• The model reaches a high accuracy before 150 epochs.
• The minimum loss is loss 1e-4.
• We can observe how the learning rate decreases (green line).
• The fluctuations in testing accuracy and loss reduce when the learning rate decreases.
• The learning rate decreases in steps, we can modify the function lr_step_decay for a faster or slower decay.
• Testing loss increases, and testing accuracy does not improve.
• The final result for testing is better than in previous cases:
  loss: 5.017
Accuracy: 0.572
Balance Accuracy: 0.389

def get_new_model_v2():
    model = tf.keras.Sequential([
            keras.Input(shape=(696, 1)),  
            keras.layers.Conv1D(32, 5, padding="valid", use_bias=True, activation='relu'),
            keras.layers.MaxPooling1D(pool_size=2),   
            keras.layers.Conv1D(64, 5, padding="valid", use_bias=True, activation='relu'),
            keras.layers.MaxPooling1D(pool_size=2),   
            keras.layers.Conv1D(128, 5, padding="valid", use_bias=True, activation='relu'),
            keras.layers.MaxPooling1D(pool_size=2),
            keras.layers.Flatten(),
            keras.layers.Dense(256, activation='relu'),
            keras.layers.Dropout(.2),
            keras.layers.Dense(4, activation='softmax')
        ])
    return model

Download and load the data

We will use the Diabetic Retinopathy Detection challenge dataset. It is a large set of high-resolution retina images.
Visit session3.2_get_data_from_the_internet for more information about how to download datasets from kaggle.
The dataset is divided into training, validation, and testing so that no cross-validation is needed if we want to participate in the challenge.

The dataset has four classes rated by a clinician. The presence of diabetic retinopathy is graded from 0 to 4 according to the following scale:

As is typical in challengers, we do not have the test data values, so we are going to split the training set to test our models.
We first visualize the distribution of the classes.

The dataset is highly unbalanced. There are some alternatives to handle this type of problem. We will use only a part of the data of class 0 for training (5000 images).
We can not use global accuracy metrics. However, we can use sensitivity which gives us a metric regarding the classifications of each class.
During training, we evaluate our model on validation data every n epochs to speed up the training process.

Prepare raw data

Follow this Jupyter notebook to prepare the data.
In preprocessing, we crop the region of interest from the images. Then we resize it to a square image by padding the smaller dimension with zeros.
We use the Contrast Limited Adaptive Histogram (CLAHE) algorithm to enhance local contrast.
Finally, we save the images into folders corresponding to each class.

Model definition

The Jupyter notebook can be found here .

We define a simple model using a pretrained model (MobileNetV2). We set the image size to 224,224,3. The network freezes the base model parameters and uses the GlobalMaxPooling2D to generate a global feature that passes through three dense layers.

def get_pretrain_model(network='VGG16', nClass=2, img_width=32, img_height=32):
    """
    This function returns  a simple 2D CNN model.
    It use as base model a pretrain model
    network (string) : select an option 'VGG16', 'MobileNetV2', 'ResNet50' or'InceptionV3'
    nClass (int) : number of classes
    img_widht (int) : input image widht
    img_height (int) : input image height

    return a tensorflow model
    """
    # Define model
    if network=='VGG16':
        base_model = vgg16.VGG16(weights='imagenet',include_top=False, input_shape=(img_height, img_width,3))
    elif network =='MobileNetV2':
        base_model = mobilenet_v2.MobileNetV2(weights='imagenet',include_top=False, input_shape=(img_height, img_width, 3))
    elif network == 'ResNet50':
        base_model = resnet50.ResNet50(weights='imagenet',include_top=False, input_shape=(img_height, img_width, 3))
    elif network == 'InceptionV3':
        base_model = inception_v3.InceptionV3(weights='imagenet',include_top=False, input_shape=(img_height, img_width, 3))

    base_model.trainable = False
    x = base_model.output
    x = GlobalMaxPooling2D()(x)
    x = Dense(2048, activation='relu')(x)
    x = Dense(256, activation="relu" )(x)
    x = Dense(nClass, activation='softmax')(x)

    return Model(base_model.input, x)

Neither the network configuration nor the final classification is our primary interest in this practice. Our goal is to recognize when the learning rate is adequate or not.

Define a number of epochs

There is no rule to define the number of necessary epochs. Some conditions require more epochs:

• CNN with several layers to train.
• A small learning rate.
• A large dataset.

In this case, the network is extensive. However, the number of learnable parameters is small, and we do not require many epochs. We set epochs = 1000 to check how the training parameters respond (loss and accuracy). Later we reduce this number.

Learning rate too Large, learning_rate = 0.01

What can we observe on the plots?

• The model accuracy decreases over time, and the loss increases.
• The learning rate is too high, producing large fluctuations in the model accuracies and losses.
• Around the epoch 400, the model diverges.
• It is a better idea to set early stopping to avoid running unnecessary iterations
• The results:
  loss: 1.269
Sensitivity: 0.0