<jupyter_start><jupyter_text>Deep Learning  for Satellite Image Classification (Manning Publications)by <NAME>        > Chapter 5: Model Optimization                            What you learned in Part 4.      In Part 4, you trained U-Net models to segment water pixels in imagery            What you will learn in this Part.      In Part 5, you will optimize the models you made in Part 4. You will first learn how to tune hyper-parameters, add regularization and modify model architectures in several ways, and deal with class imbalance. Then, you will learn how to use a machine learning model known as a ‘fully-connected conditional random field” to refine label images, so labels are as accurate as possible.          Preliminaries for ColabLike Part 3 and 4 previously, below are some convenience functions for those working on Google Colab with a GPU runtime<jupyter_code>colab = 0
#colab = 1

if colab==1:
    %tensorflow_version 2.x
    #!pip install --default-timeout=1000 tensorflow-gpu==2.0<jupyter_output><empty_output><jupyter_text>You may have to restart the runtime and/or change runtime type here, if the following doesn't show Tensorflow version 2, and a GPU available<jupyter_code>import tensorflow as tf
print(tf.__version__)
print(tf.test.is_gpu_available())<jupyter_output><empty_output><jupyter_text>Convenience functions if you need to download example (minimal) imagery sets derived from NWPU and Sentinel-2 cloudless:<jupyter_code># from https://stackoverflow.com/questions/38511444/python-download-files-from-google-drive-using-url
import requests

def download_file_from_google_drive(id, destination):
    URL = "https://docs.google.com/uc?export=download"

    session = requests.Session()

    response = session.get(URL, params = { 'id' : id }, stream = True)
    token = get_confirm_token(response)

    if token:
        params = { 'id' : id, 'confirm' : token }
        response = session.get(URL, params = params, stream = True)

    save_response_content(response, destination)    

def get_confirm_token(response):
    for key, value in response.cookies.items():
        if key.startswith('download_warning'):
            return value

    return None

def save_response_content(response, destination):
    CHUNK_SIZE = 32768

    with open(destination, "wb") as f:
        for chunk in response.iter_content(CHUNK_SIZE):
            if chunk: # filter out keep-alive new chunks
                f.write(chunk)


#s2 cloudless imagery
file_id = '1iMfIjr_ul49Ghs2ewazjCt8HMPfhY47h'
destination = 's2cloudless_imagery.zip'
if colab==1:
    download_file_from_google_drive(file_id, destination)

#s2 cloudless labels
file_id = '1c7MpwKVejoUuW9F2UaF_vps8Vq2RZRfR'
destination = 's2cloudless_label_imagery.zip'
if colab==1:
    download_file_from_google_drive(file_id, destination)

#nwpu imagery
file_id = '1gtuqy1VlU8-M5IEMnmiSuTlI5PxQPnGB'
destination = 'nwpu_images.zip'
if colab==1:
    download_file_from_google_drive(file_id, destination)

#nwpu labels
file_id = '1W5LGbcYAcFbG5YjLgX_ekBn0u5Rno35x'
destination = 'nwpu_label_images.zip'
if colab==1:
    download_file_from_google_drive(file_id, destination)                        
import zipfile
def unzip(f):
    """
    f = file to be unzipped
    """    
    with zipfile.ZipFile(f, 'r') as zip_ref:
        zip_ref.extractall()
        
if colab==1:
    unzip('s2cloudless_imagery.zip')
    unzip('s2cloudless_label_imagery.zip')   
    unzip('nwpu_images.zip')
    unzip('nwpu_label_images.zip')<jupyter_output><empty_output><jupyter_text>Import the libraries we will need:<jupyter_code>%matplotlib inline
import tensorflow as tf
from tensorflow.keras.layers import Input, Conv2D, MaxPooling2D, UpSampling2D
from tensorflow.keras.layers import Concatenate, Conv2DTranspose
from tensorflow.keras.models import Model
import numpy as np
import json, os
from random import shuffle
from PIL import Image
import matplotlib
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Training tuning strategiesWe will set up a baseline model with no optimization using the NWPU imagery, then explore the following training optimization strategies:  Building a bigger model with more layers      Using Early Stopping and Adaptive Learning Rates   Using a bigger model (and dropout)   Using regularization (Batch Normalization)   Using residual connections                   Getting things set up with a baseline modelZip through all these functions that we defined in the last Part. Like in the last Part, we will define the ```unet``` model, create model training callbacks, and generate augmented imagery<jupyter_code>def mean_iou(y_true, y_pred):
    yt0 = y_true[:,:,:,0]
    yp0 = tf.keras.backend.cast(y_pred[:,:,:,0] > 0.5, 'float32')
    inter = tf.math.count_nonzero(tf.logical_and(tf.equal(yt0, 1), tf.equal(yp0, 1)))
    union = tf.math.count_nonzero(tf.add(yt0, yp0))
    iou = tf.where(tf.equal(union, 0), 1., tf.cast(inter/union, 'float32'))
    return iou

def unet(sz = (512, 512, 3)):
  inputs = Input(sz)
  _ = inputs
  
  #down sampling 
  f = 8
  layers = []
  
  for i in range(0, 6):
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    layers.append(_)
    _ = MaxPooling2D() (_)
    f = f*2
  ff2 = 64 
  
  #bottleneck 
  j = len(layers) - 1
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2DTranspose(ff2, 2, strides=(2, 2), padding='same') (_)
  _ = Concatenate(axis=3)([_, layers[j]])
  j = j -1 
  
  #upsampling 
  for i in range(0, 5):
    ff2 = ff2//2
    f = f // 2 
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2DTranspose(ff2, 2, strides=(2, 2), padding='same') (_)
    _ = Concatenate(axis=3)([_, layers[j]])
    j = j -1 
    
  #classification 
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  outputs = Conv2D(1, 1, activation='sigmoid') (_)
  
  #model creation 
  model = Model(inputs=[inputs], outputs=[outputs])
  model.compile(optimizer = 'rmsprop', loss = 'binary_crossentropy', metrics = [mean_iou])
  
  return model 

model = unet()      
# inheritance for training process plot 
class PlotLearning(tf.keras.callbacks.Callback):

    def on_train_begin(self, logs={}):
        self.i = 0
        self.x = []
        self.losses = []
        self.val_losses = []
        self.acc = []
        self.val_acc = []
        #self.fig = plt.figure()
        self.logs = []
    def on_epoch_end(self, epoch, logs={}):
        self.logs.append(logs)
        self.x.append(self.i)
        self.losses.append(logs.get('loss'))
        self.val_losses.append(logs.get('val_loss'))
        self.acc.append(logs.get('mean_iou'))
        self.val_acc.append(logs.get('val_mean_iou'))
        self.i += 1
        print('i=',self.i,'loss=',logs.get('loss'),'val_loss=',logs.get('val_loss'),'mean_iou=',logs.get('mean_iou'),'val_mean_iou=',logs.get('val_mean_iou'))
        
        #choose a random test image and preprocess
        path = np.random.choice(test_files)
        infile = f'nwpu_images/data/{path}'
        raw = Image.open(infile)
        raw = np.array(raw.resize((512, 512)))/255.
        raw = raw[:,:,0:3]
        
        #predict the mask 
        pred = 255*model.predict(np.expand_dims(raw, 0)).squeeze()
        print(np.max(pred))
                
        #mask post-processing 
        msk  = (pred>60).astype('int') #100       
        msk = np.stack((msk,)*3, axis=-1)
        #msk[msk >= 0.5] = 1 
        #msk[msk < 0.5] = 0 
        
        #show the mask and the segmented image 
        combined = np.concatenate([raw, msk, raw* msk], axis = 1)
        plt.axis('off')
        plt.imshow(combined)
        plt.show()<jupyter_output><empty_output><jupyter_text>We use Keras callbacks to implement learning rate decay if the validation loss does not improve for 5 continues epochs. Called "reduce loss on plateau"Also, we implement early stopping if the validation loss does not improve for 5 continuous epochs.<jupyter_code>def image_batch_generator(files, batch_size = 32, sz = (512, 512)):
  
  while True: # this is here because it will be called repeatedly by the training function
    
    #extract a random subset of files of length "batch_size"
    batch = np.random.choice(files, size = batch_size)    
    
    #variables for collecting batches of inputs (x) and outputs (y)
    batch_x = []
    batch_y = []
    
    #cycle through each image in the batch
    for f in batch:

        #preprocess the raw images 
        rawfile = f'nwpu_images/data/{f}'
        raw = Image.open(rawfile)
        raw = raw.resize(sz)
        raw = np.array(raw)

        #check the number of channels because some of the images are RGBA or GRAY
        if len(raw.shape) == 2:
            raw = np.stack((raw,)*3, axis=-1)

        else:
            raw = raw[:,:,0:3]
            
        #get the image dimensions, find the min dimension, then square the image off    
        nx, ny, nz = np.shape(raw)
        n = np.minimum(nx,ny)
        raw = raw[:n,:n,:] 
            
        batch_x.append(raw)
        
        #get the masks. 
        maskfile = rawfile.replace('nwpu_images','nwpu_label_images')+'_mask.jpg'
        mask = Image.open(maskfile)
        # the mask is 3-dimensional so get the max in each channel to flatten to 2D
        mask = np.max(np.array(mask.resize(sz)),axis=2)
        # water pixels are always greater than 100
        mask = (mask>200).astype('int')
        
        mask = mask[:n,:n]

        batch_y.append(mask)

    #preprocess a batch of images and masks 
    batch_x = np.array(batch_x)/255. #divide image by 255 to normalize
    batch_y = np.array(batch_y)
    batch_y = np.expand_dims(batch_y,3) #add singleton dimension to batch_y

    yield (batch_x, batch_y) #yield both the image and the label together<jupyter_output><empty_output><jupyter_text>We've seen code like the below in the previous Part, setting up batch size, proportion of the dataset to train with, getting randomized lists of test and train file names, and finally setting up generators for model training and testing<jupyter_code>batch_size = 8

prop_train = 0.6

all_files = os.listdir('nwpu_images/data')
shuffle(all_files)

split = int(prop_train * len(all_files))

#split into training and testing
train_files = all_files[0:split]
test_files  = all_files[split:]

train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)<jupyter_output><empty_output><jupyter_text>A customary check that things worked<jupyter_code>x, y = next(train_generator)
plt.imshow(x[0])
plt.imshow(y[0].squeeze(), cmap='gray', alpha=0.5)<jupyter_output><empty_output><jupyter_text>Optimization strategy: Using Early Stoppingand Adaptive Learning Rates    Neural networks can overfit if they train for too long, but it's hard to know how many epochs is too many without a lot of trial and error, which is not efficient. Automatically computing how many training epochs a model needs or, to put it another way, evaluating when to stop training or early stopping, is achieved by monitoring validation loss during training. If loss doesn't improve for a certain number of epochs, called the patience, model training is terminated. Below patience is set to 5. Another common training strategy to prevent overly long training times is to adaptively change the learning rate during model training. This is achieved by starting with a relatively large number and then decreasing it after every training epoch until a specified minimum min_lr at a certain rate of reduction min_delta.<jupyter_code>from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping, ModelCheckpoint
def build_callbacks(filepath, min_delta, min_lr, factor):

    earlystop = EarlyStopping(monitor="val_loss", 
                                  mode="min", patience=5) 
    
    # reduction of learning rate if and when validation scores plateau upon successive epochs
    reduceloss_plat = ReduceLROnPlateau(monitor='val_loss', factor=factor, patience=5, 
                                    verbose=1, mode='auto', min_delta=min_delta, 
                                    cooldown=5, min_lr=min_lr)

    # set checkpoint file 
    model_checkpoint = ModelCheckpoint(filepath, monitor='val_loss', 
                                   verbose=0, save_best_only=True, mode='min', 
                                   save_weights_only = True)
        
    callbacks = [model_checkpoint, reduceloss_plat, earlystop, PlotLearning()]

    return callbacks  
# a tolerance for the training.
min_delta = 0.0001

# minimum learning rate (lambda)
min_lr = 0.0001

# the factor applied to the learning rate when the appropriate triggers are made
factor = 0.8

batch_size = 8

filepath = 'unet'+str(batch_size)+'.h5'
train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)
train_steps = len(train_files) //batch_size
test_steps = len(test_files) //batch_size
print(train_steps)
print(test_steps)
history = model.fit_generator(train_generator, 
                    epochs = 50, steps_per_epoch = train_steps,
                    validation_data = test_generator, validation_steps = test_steps,
                    callbacks = build_callbacks(filepath, min_delta, min_lr, factor), verbose = 0,
                    use_multiprocessing=True)
del model<jupyter_output><empty_output><jupyter_text>Optimization strategy: changing model architecture use a bigger rocket    We'll create a new `unet` function that includes initial filter size f and "bottleneck-creation filter size" ff2. The the previous models, f=8 and ff2=64. This time we make that a user-decision. We'll use f=16 and ff2=128. Because this will massively increase the number of tunable parameters, we will implement a regularization strategy called Dropout. Regularization refers to a training strategy that prevents the model from overfitting the data. This ensures the model generalizes well and is useful for prediction on unseen test data. Dropout is a popular regularization trick, which involves randomly dropping neurons of output layers. The hyperparameter that governs its behavior is a rate. 0.25 means that 25% of neurons, randomly selected, will be dropped. Dropout forces the model to learn robust features, because dropout may deactivate a neuron heavily weighted for a specific feature seen only in training imagery. Dropout thus ensures that more neurons share feature extraction tasks.<jupyter_code>from tensorflow.keras.layers import Dropout
def unet(sz, f, ff2):
  inputs = Input(sz)
  _ = inputs
  
  #down sampling 
  #f = 8
  layers = []
  
  for i in range(0, 6):
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    layers.append(_)
    _ = MaxPooling2D() (_)
    _ = Dropout(0.25)(_) #added Dropout layer
    f = f*2
  #ff2 = 64 
  
  #bottleneck 
  j = len(layers) - 1
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2DTranspose(ff2, 2, strides=(2, 2), padding='same') (_)
  _ = Concatenate(axis=3)([_, layers[j]])
  j = j -1 
      
  #upsampling 
  for i in range(0, 5):
    ff2 = ff2//2
    f = f // 2 
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2DTranspose(ff2, 2, strides=(2, 2), padding='same') (_)
    _ = Concatenate(axis=3)([_, layers[j]])
    _ = Dropout(0.25)(_) #added Dropout layer   
    j = j -1 
    
  #classification 
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  outputs = Conv2D(1, 1, activation='sigmoid') (_)
  
  #model creation 
  model = Model(inputs=[inputs], outputs=[outputs])
  model.compile(optimizer = 'rmsprop', loss = 'binary_crossentropy', metrics = [mean_iou])
  
  return model 

model = unet((512, 512, 3), 16, 128)
model.summary()<jupyter_output><empty_output><jupyter_text>From 6 to 26 million parameters!<jupyter_code># since the model is much larger the batch size should be smaller so training fits in GPU memory
batch_size = 2
# define file path
filepath = 'large_unet'+str(batch_size)+'.h5'

train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)
train_steps = len(train_files) //batch_size
test_steps = len(test_files) //batch_size
print(train_steps)
print(test_steps)
large_history = model.fit_generator(train_generator, 
                    epochs = 50, steps_per_epoch = train_steps,
                    validation_data = test_generator, validation_steps = test_steps,
                    callbacks = build_callbacks(filepath, min_delta, min_lr, factor), verbose = 0,
                    use_multiprocessing=True)
del model
# summarize history for iou
plt.figure(figsize=(10,10))
plt.subplot(121)
plt.plot(history.history['mean_iou'],'k',lw=1)
plt.plot(history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Baseline IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.subplot(122)
plt.plot(large_history.history['mean_iou'],'k',lw=1)
plt.plot(large_history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Large model IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>Optimization strategy: changing model architecture using Batch NormalizationThis operation takes the outputs of the previous layer and substracts the mean and divides by the batch standard deviation. In other words, it normalizes each sample according to the standard deviation of the batch. Because batches are randomly drawn, the result is that each individual image is normalized slightly differently each time it is included in a batch. Network training benefits because this is another regularization strategy; the model learsn how to deal with intensity fluctuations in input image batches                 Update the model with batch normalization, so start by importing the relevant layer from keras<jupyter_code>from tensorflow.keras.layers import BatchNormalization<jupyter_output><empty_output><jupyter_text>The dropout layer (with a rate of 0.25) is added after each max pooling layer in the encoder portion of the network, and  then after the concatenation layer in the decoder portion of the net<jupyter_code>def unet(sz, f, ff2):
  inputs = Input(sz)
  _ = inputs
  
  #down sampling 
  #f = 8
  layers = []
  
  for i in range(0, 6):
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = BatchNormalization()(_) # added BN layer
    layers.append(_)
    _ = MaxPooling2D() (_)
    _ = Dropout(0.25)(_) #added Dropout layer
    f = f*2
  #ff2 = 64 
  
  #bottleneck 
  j = len(layers) - 1
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2DTranspose(ff2, 2, strides=(2, 2), padding='same') (_)
  _ = Concatenate(axis=3)([_, layers[j]])
  j = j -1 
      
  #upsampling 
  for i in range(0, 5):
    ff2 = ff2//2
    f = f // 2 
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2D(f, 3, activation='relu', padding='same') (_)
    _ = Conv2DTranspose(ff2, 2, strides=(2, 2), padding='same') (_)
    _ = Concatenate(axis=3)([_, layers[j]])
    _ = Dropout(0.25)(_) #added Dropout layer   
    j = j -1 
    
  #classification 
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  _ = Conv2D(f, 3, activation='relu', padding='same') (_)
  outputs = Conv2D(1, 1, activation='sigmoid') (_)
  
  #model creation 
  model = Model(inputs=[inputs], outputs=[outputs])
  model.compile(optimizer = 'rmsprop', loss = 'binary_crossentropy', metrics = [mean_iou])
  
  return model 
model = unet((512, 512, 3), 8, 64)
model.summary()
batch_size = 2

filepath = 'bn_unet'+str(batch_size)+'.h5'

train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)
train_steps = len(train_files) //batch_size
test_steps = len(test_files) //batch_size
print(train_steps)
print(test_steps)

bn_history = model.fit_generator(train_generator, 
                    epochs = 50, steps_per_epoch = train_steps,
                    validation_data = test_generator, validation_steps = test_steps,
                    callbacks = build_callbacks(filepath, min_delta, min_lr, factor), verbose = 0,
                    use_multiprocessing=True)
del model<jupyter_output><empty_output><jupyter_text>Compare the training histories of the baseline and batch normalization models<jupyter_code># summarize history for iou
plt.figure(figsize=(10,10))
plt.subplot(121)
plt.plot(history.history['mean_iou'],'k',lw=1)
plt.plot(history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Baseline IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.subplot(122)
plt.plot(bn_history.history['mean_iou'],'k',lw=1)
plt.plot(bn_history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Baseline model with batch normalizaton IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>Optimization strategy: changing model architecture using residual connectionsA standard approach is to pass the input image goes through multiple convolutions and obtain high-level features. In a network architecture with residual layers or connections, each layer gets to see both the output from the previous layer (standard) as well as the inputs to that layer. So it not only sees the ouputs but also the data used to learn that output                 ![](https://cdn-images-1.medium.com/max/1000/1*4wx7szWCBse9-7eemGQJSw.png)Import `Activation` and `Add` layers from keras<jupyter_code>from tensorflow.keras.layers import Activation, Add<jupyter_output><empty_output><jupyter_text>Create a new UNet model function. This time we'll use a few convenience functions<jupyter_code>def batchnorm_act(x):
    x = BatchNormalization()(x)
    return Activation("relu")(x)

def conv_block(x, filters, kernel_size=(3, 3), padding="same", strides=1):
    conv = batchnorm_act(x)
    return Conv2D(filters, kernel_size, padding=padding, strides=strides)(conv)

def bottleneck_block(x, filters, kernel_size=(3, 3), padding="same", strides=1):
    conv = Conv2D(filters, kernel_size, padding=padding, strides=strides)(x)
    conv = conv_block(conv, filters, kernel_size=kernel_size, padding=padding, strides=strides)
    
    bottleneck = Conv2D(filters, kernel_size=(1, 1), padding=padding, strides=strides)(x)
    bottleneck = batchnorm_act(bottleneck)
    
    return Add()([conv, bottleneck])

def res_block(x, filters, kernel_size=(3, 3), padding="same", strides=1):
    res = conv_block(x, filters, kernel_size=kernel_size, padding=padding, strides=strides)
    res = conv_block(res, filters, kernel_size=kernel_size, padding=padding, strides=1)
    
    bottleneck = Conv2D(filters, kernel_size=(1, 1), padding=padding, strides=strides)(x)
    bottleneck = batchnorm_act(bottleneck)
    
    return Add()([bottleneck, res])

def upsamp_concat_block(x, xskip):
    u = UpSampling2D((2, 2))(x)
    return Concatenate()([u, xskip])

def res_unet(sz, f):
    inputs = Input(sz)
    
    ## downsample  
    e1 = bottleneck_block(inputs, f); f = int(f*2)
    e2 = res_block(e1, f, strides=2); f = int(f*2)
    e3 = res_block(e2, f, strides=2); f = int(f*2)
    e4 = res_block(e3, f, strides=2); f = int(f*2)
    _ = res_block(e4, f, strides=2)
    
    ## bottleneck
    b0 = conv_block(_, f, strides=1)
    _ = conv_block(b0, f, strides=1)
    
    ## upsample
    _ = upsamp_concat_block(_, e4)
    _ = res_block(_, f); f = int(f/2)
    
    _ = upsamp_concat_block(_, e3)
    _ = res_block(_, f); f = int(f/2)
    
    _ = upsamp_concat_block(_, e2)
    _ = res_block(_, f); f = int(f/2)
    
    _ = upsamp_concat_block(_, e1)
    _ = res_block(_, f)
    
    ## classify
    outputs = Conv2D(1, (1, 1), padding="same", activation="sigmoid")(_)
    
    #model creation 
    model = Model(inputs=[inputs], outputs=[outputs])
    return model
model = res_unet((512, 512, 3), 16)
model.compile(optimizer = 'rmsprop', loss = 'binary_crossentropy', metrics = [mean_iou])
model.summary()<jupyter_output><empty_output><jupyter_text>Train the model with the residual connections:<jupyter_code>batch_size = 2

filepath = 'res_unet'+str(batch_size)+'.h5'

train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)
train_steps = len(train_files) //batch_size
test_steps = len(test_files) //batch_size
print(train_steps)
print(test_steps)

res_history = model.fit_generator(train_generator, 
                    epochs = 100, steps_per_epoch = train_steps,
                    validation_data = test_generator, validation_steps = test_steps,
                    callbacks = build_callbacks(filepath, min_delta, min_lr, factor), verbose = 0,
                    use_multiprocessing=True)
del model<jupyter_output><empty_output><jupyter_text>Compare the training histories of the baseline and U-Net with residual layers:<jupyter_code># summarize history for iou
plt.figure(figsize=(10,10))
plt.subplot(121)
plt.plot(history.history['mean_iou'],'k',lw=1)
plt.plot(history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Baseline IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.subplot(122)
plt.plot(res_history.history['mean_iou'],'k',lw=1)
plt.plot(res_history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Res-UNet IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>Optimization strategy: Dealing with class imbalance  using dice loss    the network tends to ignore smaller classes. A soft dice loss could be used to train a model. Unlike the IoU metric, the numerator is the number of correctly classified pixels, and the denominator is the total number of pixels in a class that is in both estimated and ground truth. Thus it is insensitive to the number of pixels total in each class.<jupyter_code>from tensorflow.keras.layers import Flatten
smooth = 1.

def dice_coef(y_true, y_pred):
    y_true_f = tf.reshape(tf.dtypes.cast(y_true, tf.float32), [-1])
    y_pred_f = tf.reshape(tf.dtypes.cast(y_pred, tf.float32), [-1])
    intersection = tf.reduce_sum(y_true_f * y_pred_f)
    return (2. * intersection + smooth) / (tf.reduce_sum(y_true_f) + tf.reduce_sum(y_pred_f) + smooth)

def dice_coef_loss(y_true, y_pred):
    return 1.0 - dice_coef(y_true, y_pred)
model = res_unet((512, 512, 3), 32)
model.compile(optimizer = 'rmsprop', loss = dice_coef_loss, metrics = [dice_coef])
model.summary()<jupyter_output><empty_output><jupyter_text>Before we start training we need to update out `PlotLearning` class with the new loss metric<jupyter_code># inheritance for training process plot 
class PlotLearning(tf.keras.callbacks.Callback):

    def on_train_begin(self, logs={}):
        self.i = 0
        self.x = []
        self.losses = []
        self.val_losses = []
        self.acc = []
        self.val_acc = []
        #self.fig = plt.figure()
        self.logs = []
    def on_epoch_end(self, epoch, logs={}):
        self.logs.append(logs)
        self.x.append(self.i)
        self.losses.append(logs.get('loss'))
        self.val_losses.append(logs.get('val_loss'))
        self.acc.append(logs.get('dice_coef'))
        self.val_acc.append(logs.get('val_dice_coef'))
        self.i += 1
        print('i=',self.i,'loss=',logs.get('loss'),'val_loss=',logs.get('val_loss'),'dice_coef=',logs.get('dice_coef'),'val_dice_coef=',logs.get('val_dice_coef'))
        
        #choose a random test image and preprocess
        path = np.random.choice(test_files)
        infile = f'nwpu_images/data/{path}'
        raw = Image.open(infile)
        raw = np.array(raw.resize((512, 512)))/255.
        raw = raw[:,:,0:3]
        
        #predict the mask 
        pred = 255*model.predict(np.expand_dims(raw, 0)).squeeze()
        print(np.max(pred))
                
        #mask post-processing 
        msk  = (pred>60).astype('int') #100       
        msk = np.stack((msk,)*3, axis=-1)
        
        #show the mask and the segmented image 
        combined = np.concatenate([raw, msk, raw* msk], axis = 1)
        plt.axis('off')
        plt.imshow(combined)
        plt.show()<jupyter_output><empty_output><jupyter_text>Train the Residual-UNet model again, this time with dice loss<jupyter_code>batch_size = 2

filepath = 'res_dice_unet'+str(batch_size)+'.h5'

train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)
train_steps = len(train_files) //batch_size
test_steps = len(test_files) //batch_size
print(train_steps)
print(test_steps)

res_dice_history = model.fit_generator(train_generator, 
                    epochs = 100, steps_per_epoch = train_steps,
                    validation_data = test_generator, validation_steps = test_steps,
                    callbacks = build_callbacks(filepath, min_delta, min_lr, factor), verbose = 0,
                    use_multiprocessing=True)
res_dice_history.history.keys()
# summarize history for iou
plt.figure(figsize=(10,10))
plt.subplot(121)
plt.plot(res_history.history['mean_iou'],'k',lw=1)
plt.plot(res_history.history['val_mean_iou'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Res-UNet IoU')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.subplot(122)
plt.plot(res_dice_history.history['dice_coef'],'k',lw=1)
plt.plot(res_dice_history.history['val_dice_coef'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.title('Res-UNet dice coefficient')
plt.ylabel('IoU')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>Optimization strategy: Refining label images  using conditional random fieldsIt is common to use some form of post-processing algorithm to refine either labels for training, or model predictions, on an image-by-image basis. A popular approach is called a fully connected or dense conditional random field or CRF. Each image is paired with a label image and the label image is refined based on the image. The CRF makes a model for the label given the image, then assesses the likelihood of each label given the CRF model. The model can be based on relative color differences between a label and the CRF model's idea of what color that label is (within a tolerance), and/or relative spatial differences (i.e. location within the image) between a given label and the CRF model's idea of the location of that label, again within a tolerance.    Tip.     We're using the CRF to correct a pixelwise (dense) label image. It can also be used as an 'inpainting' algorithm if only sparse areas of the image have labels. This was the approach proposed by Buscombe and Ritchie (2018)                  Uncomment below to install the module and its dependencies:<jupyter_code>#!pip install cython
#!pip install git+https://github.com/lucasb-eyer/pydensecrf.git
from pydensecrf import densecrf
from pydensecrf.utils import unary_from_labels
import numpy as np<jupyter_output><empty_output><jupyter_text>The function below will refine a label image based on an undelying image, using a dense CRF. The model's hyperparameters (i.e. specified by you, not by training) are: ```compat_spat``` is a non-dimensional parameter that penalizes small pieces of segmentation that are spatially isolated -- it enforces more spatially consistent segmentations. Larger values means larger pieces of segmentation are allowed.```compat_col``` is a non-dimensional parameter that penalizes pieces of segmentation that are less uniform in color -- it enforces more consistent segmentations in colorspace. Larger values means pieces of segmentation with less similar image intesity are allowed.`theta_spat` and `theta_col` are tolerances in location and intensity, respectively. Larger values means pixel pairs can be considered to be the same class label with less similar location or intensity.`num_iter` is the number of iterations to run the CRF model inference<jupyter_code>def crf_labelrefine(input_image, predicted_labels):
    
    compat_spat=10
    compat_col=30
    theta_spat = 20
    theta_col = 80
    num_iter = 7
    
    h, w = input_image.shape[:2] #get image dimensions
    
    d = densecrf.DenseCRF2D(w, h, 2) #create a CRF object

    # For the predictions, densecrf needs 'unary potentials' which are labels (water or no water)
    predicted_unary = unary_from_labels(predicted_labels.astype('int')+1, num_classes, gt_prob= 0.51)
    
    # set the unary potentials to CRF object
    d.setUnaryEnergy(predicted_unary)

    # to add the color-independent term, where features are the locations only:
    d.addPairwiseGaussian(sxy=(theta_spat, theta_spat), compat=compat_spat, kernel=densecrf.DIAG_KERNEL,
                          normalization=densecrf.NORMALIZE_SYMMETRIC)

    input_image_uint = (input_image * 255).astype(np.uint8) #enfore unsigned 8-bit
    # to add the color-dependent term, i.e. 5-dimensional features are (x,y,r,g,b) based on the input image:    
    d.addPairwiseBilateral(sxy=(theta_col, theta_col), srgb=(5, 5, 5), rgbim=input_image_uint,
                           compat=compat_col, kernel=densecrf.DIAG_KERNEL, 
                           normalization=densecrf.NORMALIZE_SYMMETRIC)

    # Finally, we run inference to obtain the refined predictions:
    refined_predictions = np.array(d.inference(num_iter)).reshape(num_classes, h, w)
    
    # since refined_predictions will be a 2 x width x height array, 
    # each slice respresenting probability of each class (water and no water)
    # therefore we return the argmax over the zeroth dimension to return a mask
    return np.argmax(refined_predictions,axis=0)<jupyter_output><empty_output><jupyter_text>Create a model and, like before, compile it.<jupyter_code>model = res_unet((512, 512, 3), 32)
model.compile(optimizer = 'rmsprop', loss = dice_coef_loss, metrics = [dice_coef])
model.summary()<jupyter_output><empty_output><jupyter_text>Generate a random batch of images and associated label images, use the model to predict the label for the first image. Then refine that label using the CRF post-processing function. Run the cell repeatedly to see a variety of outputs and CRF refinements. How does the CRF model do? Try varying the model hyperparameters to see if you can make improvements<jupyter_code>x, y = next(train_generator)
ypred = model.predict(tf.expand_dims(x[0], axis=0))
ypred_crf = crf_labelrefine(x[0], ypred.squeeze())<jupyter_output><empty_output><jupyter_text>Use contours to visually compare the lake shorelines from ground truth (blue), model prediction (red) and CRF post-processing of the model prediction (white), then delete the model<jupyter_code>plt.figure(figsize=(10,10))
plt.imshow(x[0])
plt.contour(y[0].squeeze(), alpha=0.5, colors='b', label='Ground truth')
plt.contour(ypred.squeeze(), alpha=0.5, colors='r', label='Prediction')
plt.contour(ypred_crf.squeeze(), alpha=0.5, colors='w', label = 'CRF prediction')
plt.axis('off')
del model<jupyter_output><empty_output><jupyter_text>We've seen this function before, but this time we'll need to add a line that carries out the CRF post-processing on the input (training) mask<jupyter_code>def image_batch_generator(files, batch_size = 32, sz = (512, 512)):
  
  while True: # this is here because it will be called repeatedly by the training function
    
    #extract a random subset of files of length "batch_size"
    batch = np.random.choice(files, size = batch_size)    
    
    #variables for collecting batches of inputs (x) and outputs (y)
    batch_x = []
    batch_y = []
    
    #cycle through each image in the batch
    for f in batch:

        #preprocess the raw images 
        rawfile = f'nwpu_images/data/{f}'
        raw = Image.open(rawfile)
        raw = raw.resize(sz)
        raw = np.array(raw)

        #check the number of channels because some of the images are RGBA or GRAY
        if len(raw.shape) == 2:
            raw = np.stack((raw,)*3, axis=-1)

        else:
            raw = raw[:,:,0:3]
            
        #get the image dimensions, find the min dimension, then square the image off    
        nx, ny, nz = np.shape(raw)
        n = np.minimum(nx,ny)
        raw = raw[:n,:n,:] 
            
        batch_x.append(raw)
        
        #get the masks. 
        maskfile = rawfile.replace('nwpu_images','nwpu_label_images')+'_mask.jpg'
        mask = Image.open(maskfile)
        # the mask is 3-dimensional so get the max in each channel to flatten to 2D
        mask = np.max(np.array(mask.resize(sz)),axis=2)
        # water pixels are always greater than 100
        mask = (mask>200).astype('int')
        
        mask = mask[:n,:n]
        
        # use CRF to refine mask before it is used as a label
        mask = crf_labelrefine(raw, mask).squeeze()

        batch_y.append(mask)

    #preprocess a batch of images and masks 
    batch_x = np.array(batch_x)/255. #divide image by 255 to normalize
    batch_y = np.array(batch_y)
    batch_y = np.expand_dims(batch_y,3) #add singleton dimension to batch_y

    yield (batch_x, batch_y) #yield both the image and the label together<jupyter_output><empty_output><jupyter_text>We've also seen this function many times before, but again we'll need to add a line that carries out the CRF post-processing on the output mask (estimated by the model)<jupyter_code># inheritance for training process plot 
class PlotLearning(tf.keras.callbacks.Callback):

    def on_train_begin(self, logs={}):
        self.i = 0
        self.x = []
        self.losses = []
        self.val_losses = []
        self.acc = []
        self.val_acc = []
        #self.fig = plt.figure()
        self.logs = []
    def on_epoch_end(self, epoch, logs={}):
        self.logs.append(logs)
        self.x.append(self.i)
        self.losses.append(logs.get('loss'))
        self.val_losses.append(logs.get('val_loss'))
        self.acc.append(logs.get('dice_coef'))
        self.val_acc.append(logs.get('val_dice_coef'))
        self.i += 1
        print('i=',self.i,'loss=',logs.get('loss'),'val_loss=',logs.get('val_loss'),'dice_coef=',logs.get('dice_coef'),'val_dice_coef=',logs.get('val_dice_coef'))
        
        #choose a random test image and preprocess
        path = np.random.choice(test_files)
        infile = f'nwpu_images/data/{path}'
        raw = Image.open(infile)
        raw = np.array(raw.resize((512, 512)))/255.
        raw = raw[:,:,0:3]
        
        #predict the mask 
        pred = 255*model.predict(np.expand_dims(raw, 0)).squeeze()
        print(np.max(pred))
                
        #mask post-processing 
        msk  = (pred>60).astype('int') #100    
        
        # use CRF to refine mask before it is used as a label
        msk = crf_labelrefine(raw, msk).squeeze()
      
        msk = np.stack((msk,)*3, axis=-1)
          
        #show the mask and the segmented image 
        combined = np.concatenate([raw, msk, raw* msk], axis = 1)
        plt.axis('off')
        plt.imshow(combined)
        plt.show()<jupyter_output><empty_output><jupyter_text>Train the model like previously. This time both input and output masks of an image will be hopefully improved using the CRF algorithm<jupyter_code>batch_size = 2

filepath = 'res_dice_crf_unet'+str(batch_size)+'.h5'

train_generator = image_batch_generator(train_files, batch_size = batch_size)
test_generator  = image_batch_generator(test_files, batch_size = batch_size)
train_steps = len(train_files) //batch_size
test_steps = len(test_files) //batch_size
print(train_steps)
print(test_steps)

res_dice_history = model.fit_generator(train_generator, 
                    epochs = 100, steps_per_epoch = train_steps,
                    validation_data = test_generator, validation_steps = test_steps,
                    callbacks = build_callbacks(filepath, min_delta, min_lr, factor), verbose = 0,
                    use_multiprocessing=True)<jupyter_output><empty_output><jupyter_text>Plot the loss and validation loss, and dice coefficients to check the model training outputs<jupyter_code># summarize history for iou
plt.figure(figsize=(10,10))
plt.subplot(121)
plt.plot(res_history.history['loss'],'k',lw=1)
plt.plot(res_history.history['val_loss'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.ylabel('Loss')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.subplot(122)
plt.plot(res_dice_history.history['dice_coef'],'k',lw=1)
plt.plot(res_dice_history.history['val_dice_coef'],'r',lw=1)
plt.ylim(0,1)
plt.axhline(y=0.85)
plt.ylabel('Dice Coefficient')
plt.xlabel('Epoch number')
plt.legend(['train', 'test'], loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>Optimization strategy:  using ensemble predictions    It is common to train several models and average their predictions. This is called ensemble modeling because you are using the group of models to make a single prediction.     Tip.     We have several trained models. Here we will apply them all to the same imagery and average the masks, to see if that produces a more stable estimate                 Define two models, compile them, and assign them the weights from two different 'h5' files saved during model training using the `load_weights` utility<jupyter_code>model1 = res_unet((512, 512, 3), 32)
model1.compile(optimizer = 'rmsprop', loss = dice_coef_loss, metrics = [dice_coef])
model1.load_weights('weights/res_dice_crf_unet2.h5')

model2 = res_unet((512, 512, 3), 32)
model2.compile(optimizer = 'rmsprop', loss = dice_coef_loss, metrics = [dice_coef])
model2.load_weights('weights/res_dice_unet2.h5')<jupyter_output><empty_output><jupyter_text>Define a function for 1) getting estimates from models 1 and 2; 2) take the pixelwise maximum of the two; 3) obtain a CRF-refined estimated water mask, and 4) compute an IOU score evaluated against the real mask<jupyter_code>def get_pred(x, y, model1, model2):
    #predict the mask 
    pred1 = model1.predict(np.expand_dims(x, 0))
    pred2 = model2.predict(np.expand_dims(x, 0))
    
    #mask post-processing 
    msk  = np.maximum(pred1.squeeze(), pred2.squeeze())
    # binarize
    msk[msk >= 0.5] = 1 
    msk[msk < 0.5] = 0
    
    # use CRF to refine mask before it is used as a label
    msk = crf_labelrefine(x.squeeze(), msk).squeeze()
        
    # return the prediction and the IOU score of the prediction
    return msk, mean_iou(y, msk)<jupyter_output><empty_output><jupyter_text>Get a batch of images and labels and run the `get_pred` command on the first pair<jupyter_code>x, y = next(image_batch_generator)
ypred = get_pred(x[0], y[0], model1, model2)<jupyter_output><empty_output><jupyter_text>Show the label mask and the ensemble estimate side-by-side<jupyter_code>plt.subplot(121)
plt.imshow(x[0])
plt.imshow(y[0].squeeze(), cmap='gray', alpha=0.5)

plt.subplot(122)
plt.imshow(x[0])
plt.imshow(ypred.squeeze(), cmap='gray', alpha=0.5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11239317208079251, 'fraction_numerical': 0.018860653322149933, 'mean_word_length': 3.816990609303341, 'pattern_counts': {'":': 0, '<': 140, '<?xml version=': 0, '>': 147, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '34', 'dir': 'jupyter-structured-clean-dedup', 'id': '183853', 'license': 'MIT', 'n_tokens_mistral': 15276, 'n_tokens_neox': 14303, 'n_words': 3941, 'path': 'dl-satellite-image/p5_optimization/project_data/Part5_Optimizing.ipynb', 'repo_name': 'catalin-h/liveproject', 'size': '44122'}