!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index
!pip install ../input/efficientnet/efficientnet-1.1.0/ -f ./ --no-index

import efficientnet.tfkeras as efn

import os
import random
import pandas as pd
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns

from sklearn.metrics import mean_absolute_error
from sklearn.model_selection import KFold, train_test_split 
import tensorflow as tf
from tensorflow.keras import Model, backend
import tensorflow.keras.layers as L
import tensorflow.keras.models as M
from tensorflow.keras.utils import Sequence
from keras.utils.vis_utils import plot_model

import pydicom
import cv2
import time

def seed_everything(seed=2020):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    tf.random.set_seed(seed)
    
seed_everything(42)

config = tf.compat.v1.ConfigProto()
config.gpu_options.allow_growth = True
session = tf.compat.v1.Session(config=config)

palette_ro = ["#ee2f35", "#fa7211", "#fbd600", "#75c731", "#1fb86e", "#0488cf", "#7b44ab"]

ROOT = "../input/osic-pulmonary-fibrosis-progression/"

time.sleep(15)

train = pd.read_csv(ROOT + "train.csv")
test = pd.read_csv(ROOT + "test.csv")
sub = pd.read_csv(ROOT + "sample_submission.csv")

print("Training data shape: ", train.shape)
print("Test data shape: ", test.shape)

train.head(10)

time.sleep(15)

# %% [markdown]
# * `Patient`
# * `Weeks`
# * `FVC`
# * `Percent`
# * `Age`
# * `Sex`（`Male` / `Female`）
# * `SmokingStatus`
# # Explore CSV data 📊

train.isnull().sum()

test.isnull().sum()

dupRows_train = train[train.duplicated(subset=['Patient', 'Weeks'], keep=False)]

print("There are {} duplicate rows here ({:.2f} percent of the total).".format(len(dupRows_train), len(dupRows_train)/len(train)*100))
dupRows_train

time.sleep(20)

train.drop_duplicates(subset=['Patient', 'Weeks'], keep=False, inplace=True)

stats = []
for col in train.columns:
    stats.append((col,
                  train[col].nunique(),
                  train[col].value_counts().index[0],
                  train[col].value_counts().values[0],
                  train[col].isnull().sum() * 100 / train.shape[0],
                  train[col].value_counts(normalize=True, dropna=False).values[0] * 100,
                  train[col].dtype))
stats_df = pd.DataFrame(stats, columns=['Feature', 'Unique values', 'Most frequent item', 'Freuquence of most frequent item', 'Percentage of missing values', 'Percentage of values in the biggest category', 'Type'])
stats_df.sort_values('Percentage of missing values', ascending=False)

data = train.groupby("Patient").first().reset_index(drop=True)
data.head()

time.sleep(15)

#Multiple Quantile Regression
sub = pd.read_csv(ROOT + "sample_submission.csv")
sub.head()

sub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])
sub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))
sub =  sub[['Patient', 'Weeks', 'Confidence', 'Patient_Week']]
sub = sub.merge(test.drop('Weeks', axis=1), on="Patient")
sub.head()

train['WHERE'] = 'train'
test['WHERE'] = 'val'
sub['WHERE'] = 'test'
data = train.append([test, sub])

print(train.shape, test.shape, sub.shape, data.shape)
print(train.Patient.nunique(), test.Patient.nunique(), sub.Patient.nunique(), data.Patient.nunique())

data.head(10)

data['min_week'] = data['Weeks']
data.loc[data.WHERE=='test', 'min_week'] = np.nan
data['min_week'] = data.groupby('Patient')['min_week'].transform('min')

data.head(10)

base = data.loc[data.Weeks == data.min_week]
base = base[['Patient', 'FVC']].copy()
base.columns = ['Patient', 'base_FVC']
base['nb'] = 1
base['nb'] = base.groupby('Patient')['nb'].transform('cumsum')
base = base[base.nb==1]
base.drop('nb', axis=1, inplace=True)

base.head()

data = data.merge(base, on='Patient', how='left')
data['base_week'] = data['Weeks'] - data['min_week']
del base

data.head(10)

categorical_features = ['Sex', 'SmokingStatus']
features_nn = []
for col in categorical_features:
    for mod in data[col].unique():
        features_nn.append(mod)
        data[mod] = (data[col] == mod).astype(int)

data.head(10)

data['Percent_n'] = (data['Percent'] - data['Percent'].min() ) / ( data['Percent'].max() - data['Percent'].min() )
data['Age_n'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )
data['base_FVC_n'] = (data['base_FVC'] - data['base_FVC'].min() ) / ( data['base_FVC'].max() - data['base_FVC'].min() )
data['base_week_n'] = (data['base_week'] - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )
features_nn += ['Age_n', 'Percent_n', 'base_week_n', 'base_FVC_n']

print(features_nn)
data.head(10)

time.sleep(25)

train = data.loc[data.WHERE=='train']
test = data.loc[data.WHERE=='val']
sub = data.loc[data.WHERE=='test']
del data

train.shape, test.shape, sub.shape

C1, C2 = tf.constant(70, dtype="float32"), tf.constant(1000, dtype="float32")
def score(y_true, y_pred):
    tf.dtypes.cast(y_true, tf.float32)
    tf.dtypes.cast(y_pred, tf.float32)
    sigma = y_pred[:, 2] - y_pred[:, 0]
    fvc_pred = y_pred[:, 1]

    sigma_clip = tf.maximum(sigma, C1)
    delta = tf.abs(y_true[:, 0] - fvc_pred)
    delta = tf.minimum(delta, C2)
    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )
    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip * sq2)
    return backend.mean(metric)
#============================#
def qloss(y_true, y_pred):
    # Pinball loss for multiple quantiles
    qs = [0.2, 0.5, 0.8]
    q = tf.constant(np.array([qs]), dtype=tf.float32)
    e = y_true - y_pred
    v = tf.maximum(q*e, (q-1)*e)
    return backend.mean(v)
#=============================#
def mloss(_lambda):
    def loss(y_true, y_pred):
        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)
    return loss

def make_model():
    inp = L.Input(len(features_nn), name="Patient")
    x = L.Dense(100, activation="relu", name="d1")(inp)
    x = L.Dense(100, activation="relu", name="d2")(x)
    p1 = L.Dense(3, activation="linear", name="p1")(x)
    p2 = L.Dense(3, activation="relu", name="p2")(x)
    preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), 
                     name="preds")([p1, p2])
    
    model = M.Model(inp, preds, name="NeuralNet")
    model.compile(loss=mloss(0.64),    # changed from 0.8
                  optimizer=tf.keras.optimizers.Adam(lr=0.1, decay=0.01),
                  metrics=[score])
    return model

model = make_model()
model.summary()

plot_model(model)

time.sleep(15)

# ## Cross validation 
X_train = train[features_nn].values
X_test = sub[features_nn].values

y_train = train['FVC'].values

oof_train = np.zeros((X_train.shape[0], 3))
y_preds = np.zeros((X_test.shape[0], 3))

BATCH_SIZE = 128
EPOCHS = 804    # changed from 800
NFOLD = 5

kf = KFold(n_splits=NFOLD)

%time
for fold_id, (tr_idx, va_idx) in enumerate(kf.split(X_train)):
    print(f"FOLD {fold_id+1}")
    model = make_model()
    model.fit(X_train[tr_idx], y_train[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, 
              validation_data=(X_train[va_idx], y_train[va_idx]), verbose=0)
    print("train", model.evaluate(X_train[tr_idx], y_train[tr_idx], verbose=0, batch_size=BATCH_SIZE))
    print("val", model.evaluate(X_train[va_idx], y_train[va_idx], verbose=0, batch_size=BATCH_SIZE))
    oof_train[va_idx] = model.predict(X_train[va_idx], batch_size=BATCH_SIZE, verbose=0)
    y_preds += model.predict(X_test, batch_size=BATCH_SIZE, verbose=0) / NFOLD

fig, ax = plt.subplots(figsize=(12, 12))

idxs = np.random.randint(0, y_train.shape[0], 100)
ax.plot(y_train[idxs], label="ground truth", color=palette_ro[0])
ax.plot(oof_train[idxs, 0], label="q20", color=palette_ro[3], ls=':', alpha=0.5)
ax.plot(oof_train[idxs, 1], label="q50", color=palette_ro[4], ls=':', alpha=0.5)
ax.plot(oof_train[idxs, 2], label="q80", color=palette_ro[5], ls=':', alpha=0.5)
ax.legend(loc="best");

# We calculate the optimized 𝜎 (standard deviation) from the `oof_train`. `sigma_opt` is the mean absolute error between the correct value of each fold and the prediction (median), `sigma_unc` is the difference between the prediction (0.2 quantile) and the prediction (0.8 quantile), and `sigma_mean` is the mean value of the difference.<br> 

sigma_opt = mean_absolute_error(y_train, oof_train[:, 1])
sigma_unc = oof_train[:, 2] - oof_train[:, 0]
sigma_mean = np.mean(sigma_unc)
print(sigma_opt, sigma_mean)

print(sigma_unc.min(), sigma_unc.mean(), sigma_unc.max(), (sigma_unc>=0).mean())

print(np.mean(y_train / oof_train[:, 1]))

time.sleep(15)

fig, ax = plt.subplots(figsize=(16, 6))

sns.distplot(sigma_unc, ax=ax, color=palette_ro[1])
ax.set_title("uncertainty in prediction", fontsize=18);