# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

from __future__ import print_function
import os
from keras.layers import merge
from keras.layers.advanced_activations import LeakyReLU
from keras.optimizers import Nadam
from keras.metrics import MSE
from keras.layers.convolutional import Convolution2D, MaxPooling2D, ZeroPadding2D, AveragePooling2D
from keras.layers.core import Dense, Activation, Flatten
from keras.layers.normalization import BatchNormalization
from keras.models import Model
from keras.callbacks import ProgbarLogger, ModelCheckpoint
from keras.layers import Input
from keras.preprocessing.image import load_img, img_to_array, ImageDataGenerator
from sklearn.metrics import roc_auc_score, confusion_matrix, classification_report
from scipy.stats import spearmanr
import keras.backend as K
import numpy as np
import tensorflow as tf
import datetime
import scandir
import cv2
from PIL import Image
from ml_metrics import quadratic_weighted_kappa
os.environ['KERAS_BACKEND'] = 'tensorflow'
os.environ['CUDA_HOME'] = '/usr/local/cuda-7.5'


def read_img(img_path):
    '''This function returns a preprocessed image
    '''
    dim_ordering = K.image_dim_ordering()
    img = load_img(img_path, target_size=(512, 512))
    img = img_to_array(img, dim_ordering=dim_ordering)

    if dim_ordering == 'th':
        img = img[::-1, :, :]
    else:
        img = img[:, :, ::-1]
    return img


def get_patchblock(inp, idx):
    if K.image_dim_ordering() == 'tf':
        inp_shape=(512, 512, 3)
        bn_axis = 3
    else:
        inp_shape=shape=(3, 512, 512)
        bn_axis = 1
    dim_ordering = K.image_dim_ordering()
        
    with tf.device('/gpu:0'):
        out = Convolution2D(32, 7, 7, subsample=(2, 2),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv1_{0}'.format(idx), input_shape=(inp_shape))(inp)
        out = BatchNormalization(axis=bn_axis, name='bn_conv1_{0}'.format(idx))(out)
        out = Activation('relu')(out) LeakyReLU(alpha=0.5)
        out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)
        
        out = Convolution2D(32, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv2_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv2_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
        out = Convolution2D(32, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv3_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv3_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)
        
        out = Convolution2D(64, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv4_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv4_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
    with tf.device('/gpu:1'):
        out = Convolution2D(64, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv5_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv5_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)

        out = Convolution2D(128, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv6_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv6_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
        out = Convolution2D(128, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv7_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv7_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
        
    with tf.device('/gpu:2'):
        out = Convolution2D(128, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv8_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv8_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
        out = Convolution2D(128, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv9_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv9_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)

        out = Convolution2D(256, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv10_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv10_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
    with tf.device('/gpu:3'):
        out = Convolution2D(256, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv11_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv11_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
        out = Convolution2D(256, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv12_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv12_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        
        out = Convolution2D(256, 3, 3, subsample=(1, 1),
                            init='he_normal', border_mode='same', dim_ordering=dim_ordering,
                            name='conv13_{0}'.format(idx), input_shape=(inp_shape))(out)
        out = BatchNormalization(axis=bn_axis, name='bn_conv13_{0}'.format(idx))(out)
        out = Activation(LeakyReLU(alpha=0.5))(out)
        out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)
        out = MaxPooling2D((7, 7), dim_ordering=dim_ordering)(out)

    return out

patches_per_slice = 5
def get_mergenet():    
    if K.image_dim_ordering() == 'tf':
        inp_shape=(512, 512, 3)
        concat_axis = 3
    else:
        inp_shape=(3, 512, 512)
        concat_axis = 1
    
    patch_nets = []
    input_list = []
    for i in range(patches_per_slice):
        inp = Input(inp_shape)
        model = get_patchblock(inp, i)
        patch_nets.append(model)
        input_list.append(inp)
        
    out = merge(patch_nets, mode='concat', concat_axis=concat_axis)
    out = Flatten()(out)
    out = Dense(2, activation='softmax', name='fc')(out)
    return Model(input_list, out)


print('Started to build model at {0}'.format(datetime.datetime.utcnow()))
model = get_mergenet()
model.summary()
print('Finished building model at {0}'.format(datetime.datetime.utcnow()))

#TODO: insert my own loss and metrics functions (one per output), equal to the ones specified by TUPAC
print('Started to compile model at {0}'.format(datetime.datetime.utcnow()))
model.compile(optimizer='adam',
              loss='mse',
              metrics=['mse'])
print('Finished compiling model at {0}'.format(datetime.datetime.utcnow()))


def get_data(dir_path, res = 512):
    annotations = open('../annotations/training_ground_truth.csv', 'r')
    lines = annotations.readlines()
    images = []
    labels = []
    for subdir, _, files in scandir.walk(dir_path):
        subdir = subdir.replace('\\', '/')  # windows path fix
        subdir_split = subdir.split('/')
        if len(subdir_split[3])>0: study_id = int(subdir_split[3].lstrip("0"))
        else: continue
        label = [float(lines[study_id].split(',')[0]), float(lines[study_id].split(',')[1])]
        imgs = []
        for f in files:
            tiff_path = os.path.join(subdir, f)
            if not tiff_path.endswith('.tiff'):
                continue
            img = read_img(tiff_path)
            imgs.append(img)
        images.append(imgs)
        labels.append(label)
    annotations.close()
    images_list = []
    print(np.asarray(images))
    for x in np.split(np.swapaxes(images, 0, 1), 5):
        images_list.append(np.squeeze(x))
    return images_list, np.asarray(labels)


# Load Data
images_train, labels_train = get_data('../example_images/train/')
images_val, labels_val = get_data('../example_images/val/')
print(images_train[0].shape[3]!=512, len(labels_train))


# +
# featurewise standardization, normalization and augmentation (horizontal and vertical flips)
def augment_data(X):
    X = np.asarray(X).swapaxes(1, 0)
    mean = np.mean(X, axis=0)
    X -= mean
    std = np.std(X, axis=0)
    X /= (std + 1e-7)
    X = X.swapaxes(0, 1)
    if np.random.random() < 0.5:
        X = flip_axis(X, 2)
    if np.random.random() < 0.5:
        X = flip_axis(X, 3)
    return X
        
def flip_axis(X, axis):
    X = np.asarray(X).swapaxes(axis, 0)
    X = X[::-1, ...]
    X = X.swapaxes(0, axis)
    return X

print('Started data augmentation at {0}'.format(datetime.datetime.utcnow()))
images_train = augment_data(images_train)
images_val = augment_data(images_val)
print('Finished data augmentation at {0}'.format(datetime.datetime.utcnow()))
#print(np.mean(images_train[0]), np.std(images_train[0]))

# +
batch_size = 3
nb_epoch = 10
patches_per_slice = 5
#class_weights = {}
callbacks = [ProgbarLogger(),
             ModelCheckpoint('mergenet_weights.hdf5', monitor='val_loss', verbose=1,
                             save_best_only=True, mode='auto')]

print('Started fitting at {0}'.format(datetime.datetime.utcnow()))
history = model.fit(images_train, labels_train, batch_size,
                              #samples_per_epoch=images_train.shape[0],
                              #class_weights = class_weights,
                              nb_epoch=nb_epoch, verbose=1,
                              validation_data=(images_val, labels_val),
                              #nb_val_samples = 1,
                              callbacks=callbacks)
print('Finished fitting at {0}'.format(datetime.datetime.utcnow()))

# +
# Decode predictions to one hot encoding
# Use ranking approach from https://github.com/ChenglongChen/Kaggle_CrowdFlower/blob/master/BlogPost/BlogPost.md
# The predicted values are mapped to classes based on the CDF of the ref values.
hist = np.bincount(labels_train[:,0].astype(int))
cdf = np.cumsum(hist) / float(sum(hist))
np.savetxt("../example_images/train/cdf.txt", cdf)

def getScore(pred, cdf, valid=False):
    num = pred.shape[0]
    output = np.asarray([4]*num, dtype=int)
    rank = pred.argsort()
    output[rank[:int(num*cdf[0]-1)]] = 1
    output[rank[int(num*cdf[0]):int(num*cdf[1]-1)]] = 2
    output[rank[int(num*cdf[1]):int(num*cdf[2]-1)]] = 3
    if valid:
        cutoff = [ pred[rank[int(num*cdf[i]-1)]] for i in range(3) ]
        return output, cutoff
    return output


# -

print('Metrics on validation set:\n------------------------------------------------------------\n')
pred = model.predict(images_val, verbose=0)
pred_class = getScore(pred[:,0], cdf)
print('Confusion Matrix:\n', confusion_matrix(labels_val[:,0], pred_class), '\n')
print(classification_report(labels_val[:,0], pred_class), '\n')
print("Quadratic Weighted Cohen's Kappa: ", quadratic_weighted_kappa(pred[:,0], labels_val[:,0]))
print("Spearman's Correlation Coëfficient: ", spearmanr(pred[:,1], labels_val[:,1])[0])

print(history.history.keys())



StarcoderdataJupyterScriptsDedupFiltered