<jupyter_start><jupyter_text>Train NCDOT flood detection model This is a jupyter notebook for training a supervised convolutional neural network model for flood detection using traffic webcams in the NCDOT networkPart of the Flood CAM ML projecthttps://github.com/FloodCamML Code contributions:* <NAME>, Marda Science / USGS* <NAME>, UNCG Overview:> Classification targets: > 0 = not usable > 1 = not flooded> 2 = not sure> 3 = floodedApproach:1. Use transfer learning to initialize a mobilenet-v2 model with pre-trained weights (from imagenet) as feature extractor2. Add a classification head with dropout regularization3. Unfreeze all layers and train from scratch (both the feature extractor and the classification head)4. Train with 50% of the data for validation, 50% for training5. Use keras-tuner for hyperparameter (learning rate, dropout rate, and number of classifying densely connected neurons) optimization6. Use of more than 2 classes, and one-hot encoded labels, enables use of categorical cross-entropy softmax scores to be used as independent probabilities of prediction7. models are saved as human-readable json files, with h5 weights, plus a pb version of the model for deployment using R/shinyUse of automated hyperparameter tuning means this approach may be more adaptable to different datasets and classification targets Import libraries<jupyter_code>import os
import matplotlib.pyplot as plt
import numpy as np
from glob import glob
import matplotlib.cm as cm
from skimage.transform import resize
from skimage.io import imread
from sklearn.metrics import confusion_matrix
import seaborn as sns #extended functionality / style to matplotlib plots
from random import shuffle
from datetime import datetime

#import the tf stuff
import tensorflow as tf
from tensorflow.keras import layers
from tensorflow.keras import Model
from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint
from tensorflow.keras.models import model_from_json
import keras_tuner as kt<jupyter_output><empty_output><jupyter_text>Print some important things to screen. We need tf version 2.3 or greater, and a GPU for training (ideally)<jupyter_code>print("Version: ", tf.__version__)
print("Eager mode: ", tf.executing_eagerly())
print('GPU name: ', tf.config.experimental.list_physical_devices('GPU'))
print("Num GPUs Available: ", len(tf.config.experimental.list_physical_devices('GPU')))<jupyter_output><empty_output><jupyter_text>Hard-code the path to the data like I have done. Each image file name is modified to contain the class string ('Water', 'No_water', 'Not_sure', and 'Not_usable'), and placed into a subfolder of the same name<jupyter_code>root = '../../../../HX_Ted_2020_NCTC/All_Photos/Recoded'

mode='categorical'
# mode = 'binary'

if mode == 'binary':
    #put names in alphabetical order
    CLASSES = ['No_water','Water']
    class_dict={'No_water':0, 'Water':1}
    all_files = glob(root+'/Water/*.jpg')+glob(root+'/No_water/*.jpg')
    total_files = len(all_files)
else:
    #put names in alphabetical order
    CLASSES = ['Not_usable','No_water','Not_sure','Water']
    class_dict={'Not_usable':0, 'No_water':1, 'Not_sure':2,  'Water':3}
    all_files = glob(root+'/Water/*.jpg')+glob(root+'/No_water/*.jpg')+glob(root+'/Not_sure/*.jpg')+glob(root+'/Not_usable/*.jpg')
    total_files = len(all_files)

NCLASSES = len(CLASSES)

train_files = all_files[::2]
test_files = all_files[1::2]

#randomize
shuffle(train_files)
shuffle(test_files)

print('%i train files' % (len(train_files)))
print('%i test files' % (len(test_files)))<jupyter_output><empty_output><jupyter_text>These are the hyperparameters you need to specify<jupyter_code>SCRATCH = True # or False for no training from scratch
height_width = 224 # size of images (must be square)
BS = 16 #batch size
EPOCHS = 200 #maximum number of training epochs
PATIENCE = 25 #number of training epochs over which no improvement in validation loss, before terminating training
VAL_SPLIT = 0.5 #proportion of dataset to use for validation (>= 0.5 is recommended)<jupyter_output><empty_output><jupyter_text>Function for building the model and testing a range of hyperparameters<jupyter_code>##==============================================
def model_builder(hp):

    imshape = (height_width,height_width,3)

    # #base_model, no top layer, w/ imagenet weights
    # base_model = tf.keras.applications.MobileNetV2(input_shape = imshape,
    #                                              include_top = False,
    #                                              weights = None)

    #base_model, no top layer, w/ imagenet weights
    base_model = tf.keras.applications.MobileNetV2(input_shape = imshape,
                                                 include_top = False,
                                                 weights = 'imagenet')

    if SCRATCH:
        base_model.trainable = True
    else:
        base_model.trainable = False
    # base_model.summary()

    # add a new classifcation head
    final_layer = base_model.get_layer('out_relu')
    #print('shape of last layer is ', final_layer.output_shape)
    final_base_output = final_layer.output

    #add the last layer
    # GL.Av Pool
    x = layers.GlobalAveragePooling2D()(final_base_output)
    
    # Add a fully connected layer with hidden units and ReLU activation

    hp_units = hp.Int('units', min_value=128, max_value=512, step=32)

    x = layers.Dense(hp_units, activation='relu')(x)

    dropout_rate = hp.Choice('dropout_rate', values=[.4,.5,.6])

    x = layers.Dropout(dropout_rate)(x)
    
    # Add a final sigmoid layer for classification
    if mode=='binary':
        x = layers.Dense(1, activation='sigmoid')(x)
    else:
        x = layers.Dense(NCLASSES, activation='sigmoid')(x)

    model = Model(base_model.input, x)

    hp_learning_rate = hp.Choice('learning_rate', values=[1e-4, 1e-5, 1e-6])

    if mode is 'binary':
        acc_metric = tf.keras.metrics.BinaryAccuracy(name='accuracy')

        #build the model
        model.compile(loss = 'binary_crossentropy',
                      optimizer = tf.keras.optimizers.Adam(lr = hp_learning_rate),
                      metrics = acc_metric)
    else:

        loss_fn = tf.keras.losses.CategoricalCrossentropy(
            from_logits=False,
            name='categorical_crossentropy')

        model.compile(optimizer=tf.keras.optimizers.Adam(lr= hp_learning_rate), 
                      loss=loss_fn,
                      metrics=['accuracy'])

    return model<jupyter_output><empty_output><jupyter_text>Function for training the model<jupyter_code>##==============================================
def model_trainer(model,aug=False):

    # set checkpoint file
    model_checkpoint = ModelCheckpoint(weights_file, monitor='val_loss',
                                    verbose=0, save_best_only=True, mode='min',
                                    save_weights_only = True)

    callbacks = [earlystop,model_checkpoint]

    if mode is 'binary':

        if aug:
            history = model.fit(img_generator, epochs=EPOCHS, batch_size=BS, callbacks=callbacks, validation_data=val_generator)

        else:
            history = model.fit(
                x_train,
                y_train,
                batch_size=BS,
                epochs=EPOCHS,
                # We pass some validation for
                # monitoring validation loss and metrics
                # at the end of each epoch
                validation_data=(x_test, y_test),
                callbacks =[callbacks])
    else:

        if aug:
            history = model.fit(img_generator, epochs=EPOCHS, batch_size=BS, callbacks=callbacks, validation_data=val_generator)

        else:
            history = model.fit(
                x_train,
                tf.one_hot(y_train,NCLASSES),
                batch_size=BS,
                epochs=EPOCHS,
                # We pass some validation for
                # monitoring validation loss and metrics
                # at the end of each epoch
                validation_data=(x_test, tf.one_hot(y_test,NCLASSES)),
                callbacks =[callbacks])

    #look at the metrics from training
    acc = history.history['accuracy']
    val_acc = history.history['val_accuracy']
    loss = history.history['loss']
    val_loss = history.history['val_loss']

    epochs = range(len(acc))

    plt.figure()
    plt.plot(epochs, acc, 'r', label='Training accuracy')
    plt.plot(epochs, val_acc, 'b', label='Validation accuracy')
    plt.title('Training and validation accuracy')
    plt.legend(loc=0)
    # plt.show()
    plt.savefig(weights_file.replace('.h5','_acc.png'), dpi=200, bbox_inches='tight')
    plt.close()

    plt.figure()
    plt.plot(epochs, loss, 'r--', label='Training loss')
    plt.plot(epochs, val_loss, 'b--', label='Validation loss')
    plt.title('Training and validation loss')
    plt.legend(loc=0)
    # plt.show()
    plt.savefig(weights_file.replace('.h5','_loss.png'), dpi=200, bbox_inches='tight')
    plt.close()

    # history = model.fit(img_train, label_train, epochs=50, validation_split=0.2)

    val_acc_per_epoch = history.history['val_accuracy']
    best_epoch = val_acc_per_epoch.index(max(val_acc_per_epoch)) + 1
    print('Best epoch: %d' % (best_epoch,))
    return model<jupyter_output><empty_output><jupyter_text>Function for 1. reading all train and test images into memory, and standardizing them for model training2. stripping names from filenames and looking up class integer code from class dictionary<jupyter_code>##==============================================
def get_data(train_files,test_files,height_width):
    '''
    read and resize and standardize imager from files to numpy arrays
    '''
    x_train = np.zeros((len(train_files),height_width,height_width,3))
    for counter,f in enumerate(train_files):
        im = resize(imread(f), (height_width,height_width))
        x_train[counter]=standardize(im)
    x_train = x_train.astype("float32") #/ 255.0

    x_test = np.zeros((len(test_files),height_width,height_width,3))
    for counter,f in enumerate(test_files):
        im = resize(imread(f), (height_width,height_width))
        x_test[counter]=standardize(im)

    x_test = x_test.astype("float32") #/ 255.0

    y_train = []
    for f in train_files:
        y_train.append(class_dict[f.split(os.sep)[-1].split('_X_')[0]])

    y_train = np.expand_dims(y_train,-1).astype('uint8')
    y_train = np.squeeze(y_train)

    y_test = []
    for f in test_files:
        y_test.append(class_dict[f.split(os.sep)[-1].split('_X_')[0]])
    y_test = np.expand_dims(y_test,-1).astype('uint8')

    y_test = np.squeeze(y_test)
    return x_train, y_train, x_test, y_test<jupyter_output><empty_output><jupyter_text>Functions for image standardization<jupyter_code># #IMG UTILS
##==============================================
def standardize(img):
    '''
    standardization using adjusted standard deviation,
    then rescales an input dat between mn and mx
    '''
    N = np.shape(img)[0] * np.shape(img)[1]
    s = np.maximum(np.std(img), 1.0/np.sqrt(N))
    m = np.mean(img)
    img = (img - m) / s
    img = rescale(img, 0, 1)
    del m, s, N
    if np.ndim(img)!=3:
        img = np.dstack((img,img,img))

    return img

##==============================================
def rescale(dat,mn,mx):
    '''
    rescales an input dat between mn and mx
    '''
    m = min(dat.flatten())
    M = max(dat.flatten())
    return (mx-mn)*(dat-m)/(M-m)+mn<jupyter_output><empty_output><jupyter_text>Functions to compute and apply gradCAM feature importances mapping<jupyter_code>#from https://keras.io/examples/vision/grad_cam/
##==============================================
# #GRADCAM
##==============================================
def make_gradcam_heatmap(
    img_array, model, last_conv_layer_name, classifier_layer_names
):
    # First, we create a model that maps the input image to the activations
    # of the last conv layer
    last_conv_layer = model.get_layer(last_conv_layer_name)
    last_conv_layer_model = Model(model.inputs, last_conv_layer.output)

    # Second, we create a model that maps the activations of the last conv
    # layer to the final class predictions
    classifier_input = tf.keras.Input(shape=last_conv_layer.output.shape[1:])
    x = classifier_input
    for layer_name in classifier_layer_names:
        x = model.get_layer(layer_name)(x)
    classifier_model = tf.keras.Model(classifier_input, x)

    # Then, we compute the gradient of the top predicted class for our input image
    # with respect to the activations of the last conv layer
    with tf.GradientTape() as tape:
        # Compute activations of the last conv layer and make the tape watch it
        last_conv_layer_output = last_conv_layer_model(img_array)
        tape.watch(last_conv_layer_output)
        # Compute class predictions
        preds = classifier_model(last_conv_layer_output)
        top_pred_index = tf.argmax(preds[0])
        top_class_channel = preds[:, top_pred_index]

    # This is the gradient of the top predicted class with regard to
    # the output feature map of the last conv layer
    grads = tape.gradient(top_class_channel, last_conv_layer_output)
    #grads = tape.gradient(bottom_class_channel, last_conv_layer_output)

    # This is a vector where each entry is the mean intensity of the gradient
    # over a specific feature map channel
    pooled_grads = tf.reduce_mean(grads, axis=(0, 1, 2))

    # We multiply each channel in the feature map array
    # by "how important this channel is" with regard to the top predicted class
    last_conv_layer_output = last_conv_layer_output.numpy()[0]
    pooled_grads = pooled_grads.numpy()
    for i in range(pooled_grads.shape[-1]):
        last_conv_layer_output[:, :, i] *= pooled_grads[i]

    # The channel-wise mean of the resulting feature map
    # is our heatmap of class activation
    heatmap = np.mean(last_conv_layer_output, axis=-1)

    # For visualization purpose, we will also normalize the heatmap between 0 & 1
    heatmap = np.maximum(heatmap, 0) / np.max(heatmap)
    return heatmap

##==============================================
def do_gradcam_viz(img, model, outfile):
    #define last conv layers
    last_conv_layer_name = "out_relu"
    classifier_layer_names = [
        "global_average_pooling2d",
        "dense",
        "dropout",
        "dense_1",
    ]
    
    alpha=0.4

    img_array = np.expand_dims(img,axis=0)

    #from https://keras.io/examples/vision/grad_cam/
    # Make the heatmap
    try:
        heatmap = make_gradcam_heatmap(img_array, model, last_conv_layer_name, classifier_layer_names)
    except:

        classifier_layer_names = [
            "global_average_pooling2d_1",
            "dense_2",
            "dropout_1",
            "dense_3",
        ]
        heatmap = make_gradcam_heatmap(img_array, model, last_conv_layer_name, classifier_layer_names)

    # rescale image (range 0-255)
    heatmap = np.uint8(255 * heatmap)

    # use viridis for heatmap
    vir = cm.get_cmap("viridis")
    vir_colors = vir(np.arange(256))[:, :3]
    vir_heatmap = vir_colors[heatmap]

    # make the heatmap
    vir_heatmap = tf.keras.preprocessing.image.array_to_img(vir_heatmap)
    vir_heatmapI = vir_heatmap.resize((Orimg.shape[1], Orimg.shape[0]))
    vir_heatmap = tf.keras.preprocessing.image.img_to_array(vir_heatmapI)

    #put heatmpa on image
    superimposed_img = vir_heatmap * alpha + Orimg
    superimposed_img = tf.keras.preprocessing.image.array_to_img(superimposed_img)

    # Display Image, heatmap and overlay
    # Display heatmap
    plt.figure(figsize=(20,10))
    plt.subplot(121)#31)
    plt.imshow(Orimg) ; plt.title('Orig image')
    #tf.keras.preprocessing.image.load_img(impath, target_size = imsize))
    # plt.subplot(132)
    # plt.imshow(vir_heatmapI); plt.title('heatmap')
    plt.subplot(122)#33)
    plt.imshow(superimposed_img); plt.title('Superimposed GRADCAM')

    plt.savefig(outfile, dpi=200, bbox_inches='tight')
    plt.close()<jupyter_output><empty_output><jupyter_text>Function to read an image from file and use a trained model for gradCAM estimation and prediction<jupyter_code>def predict_flood(f,FloodCAMML,height_width=224):
    '''
    given file path string, f, and height_width (default=224)
    return 0 for no flood and 1 for flood
    '''
    img = resize(imread(f), (height_width,height_width))
    Orimg = imread(f)
    img_array = np.expand_dims(img,axis=0)

    try:
        do_gradcam_viz(img, FloodCAMML, outfile=f.replace('.jpg','_gradcam.png'))
    except:
        pass
    return FloodCAMML.predict(img_array, batch_size = 1, verbose = False)

##==============================================
def get_model_load_weights(weights_file):
    '''
    given file path string, weights file, loads keras model and assigns weights
    '''
    json_file = open(weights_file.replace('.h5','.json'), 'r')
    loaded_model_json = json_file.read()
    json_file.close()
    FloodCAMML = model_from_json(loaded_model_json)
    #print("Loading weights into model")
    # load weights into new model
    FloodCAMML.load_weights(weights_file)
    return FloodCAMML<jupyter_output><empty_output><jupyter_text>Training begins here Get train and test datasets (read into memory - could replace with generators later)<jupyter_code>x_train, y_train, x_test, y_test = get_data(train_files,test_files,height_width)<jupyter_output><empty_output><jupyter_text>Grab and plot a batch to visualize / check the input data are good<jupyter_code>x = x_train[:BS]
y = y_train[:BS]

plt.figure(figsize=(15, 15))
counter = 1
for im, l in zip(x,y):
    ax = plt.subplot(4, 4,counter)
    plt.imshow(im)

    l=np.int(l)
    plt.title(CLASSES[l], fontsize=9)

    plt.axis('off')
    counter +=1

# plt.show()
plt.savefig('results/NCDot_example_train_batch_'+mode+'.png',dpi=300,bbox_inches='tight')
plt.close()
del x, y<jupyter_output><empty_output><jupyter_text>Use a hyperband tuner to search the hyperparameter space<jupyter_code>tuner = kt.Hyperband(model_builder,
                     objective='val_accuracy',
                     max_epochs=10,
                     factor=3,
                     directory='tuner',
                     project_name='CamML_'+datetime.now().strftime("%Y-%m-%d-%H-%M-%S"))

earlystop = EarlyStopping(monitor="val_loss",
                              mode="min", patience=PATIENCE)

tuner.search(x_train, 
             tf.one_hot(y_train,NCLASSES), 
             epochs=EPOCHS, 
             validation_data=(x_test, tf.one_hot(y_test,NCLASSES)), 
             callbacks=[earlystop])

# models = tuner.get_best_models(num_models=2)

tuner.results_summary()<jupyter_output><empty_output><jupyter_text>Get the optimal hyperparameters and print to screen<jupyter_code>best_hps=tuner.get_best_hyperparameters(num_trials=1)[0]

print(f"""
optimal number of units in the first densely-connected layer = {best_hps.get('units')} \n
optimal learning rate for the optimizer is {best_hps.get('learning_rate')}\n
optimal dropout rate for the optimizer is {best_hps.get('dropout_rate')}\n
""")<jupyter_output><empty_output><jupyter_text>Build the best model and create a weights file name<jupyter_code>model = tuner.hypermodel.build(best_hps)

DROPOUT = best_hps.get('dropout_rate')
NUM_NEURONS = best_hps.get('units')
lr = best_hps.get('learning_rate')

if SCRATCH:
    #tl=transferlarning, mv2=MobileNetV2 feature extractor
    weights_file = 'results/scratch_mv2_bs'+str(BS)+'_drop'+str(DROPOUT)+'_nn'+str(NUM_NEURONS)+'_sz'+str(height_width)+'_lr'+str(lr)+'_val'+str(VAL_SPLIT)+'.h5'

else:
    #tl=transferlarning, mv2=MobileNetV2 feature extractor
    weights_file = 'results/tl_mv2_bs'+str(BS)+'_drop'+str(DROPOUT)+'_nn'+str(NUM_NEURONS)+'_sz'+str(height_width)+'_lr'+str(lr)+'_val'+str(VAL_SPLIT)+'.h5'

print(weights_file)<jupyter_output><empty_output><jupyter_text>Train the model, and save to json format, and finally to pb format<jupyter_code>model = model_trainer(model)

model_json = model.to_json()
with open(weights_file.replace('.h5','.json'), "w") as json_file:
    json_file.write(model_json)

# save as TF.js
# tfjs.converters.save_keras_model(model, './JSmodel')

model.save('Rmodel_'+weights_file.replace('.h5','')+datetime.now().strftime("%d-%m-%Y-%H-%M-%S"))<jupyter_output><empty_output><jupyter_text>Test the model on some sample files<jupyter_code>notsure_files = glob(root+'/Not_sure/*.jpg')

for k in np.arange(BS):
    use = np.random.randint(100)
    f = notsure_files[use]
    img = resize(imread(f), (height_width,height_width))
    Orimg = imread(f)
    outfile = 'predict'+os.sep+f.split(os.sep)[-1].replace('.jpg','_gradcam.png')

    try:
        do_gradcam_viz(img, model, outfile)
    except:
        pass<jupyter_output><empty_output><jupyter_text>Predict on each test image and construct a confusion matrix shoing one-on-one correspondences across all 4 classes<jupyter_code>y_pred = model.predict(x_test).squeeze()

if mode=='binary':
    y_pred = (y_pred>.5).astype(int)
else:
    y_pred = np.argmax(y_pred, -1)

cmat = confusion_matrix(y_test, y_pred)
cmat = cmat.astype('float') / cmat.sum(axis=1)[:, np.newaxis]

plt.figure(figsize=(8,8))
sns.heatmap(cmat,
  annot=True,
  cmap = sns.cubehelix_palette(dark=0, light=1, as_cmap=True))

tick_marks = np.arange(len(CLASSES))+.5
plt.xticks(tick_marks, [c for c in CLASSES], rotation=90,fontsize=12)
plt.yticks(tick_marks, [c for c in CLASSES],rotation=0, fontsize=12)
plt.title('N = '+str(len(y_test)), fontsize=12)

plt.savefig(weights_file.replace('.h5','_cm.png'), dpi=200, bbox_inches='tight')
plt.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12937644341801385, 'fraction_numerical': 0.010900692840646651, 'mean_word_length': 3.78475138121547, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 82, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '423979', 'license': 'MIT', 'n_tokens_mistral': 7344, 'n_tokens_neox': 6975, 'n_words': 1619, 'path': 'model/src/NCDOT_flood_detector_training.ipynb', 'repo_name': 'CoPeCOMET/NCTrafficCameras', 'size': '21667'}