{"cells":[{"metadata":{},"cell_type":"markdown","source":"Reference to this excellent notebook:\n\nhttps://www.kaggle.com/jagadish13/osic-baseline-elasticnet-eda\n\nI highly recommand read it\n\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Load Libraries\n","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold\nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom plotly.subplots import make_subplots\nimport plotly.express as px\nimport plotly.graph_objects as go","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Useful functions","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=2020):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    \nseed_everything(42)\n\ndef score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    return K.mean(metric)\n#============================#\ndef qloss(y_true, y_pred):\n    # Pinball loss for multiple quantiles\n    qs = [0.2, 0.50, 0.8]\n    q = tf.constant(np.array([qs]), dtype=tf.float32)\n    e = y_true - y_pred\n    v = tf.maximum(q*e, (q-1)*e)\n    return K.mean(v)\n#=============================#\ndef mloss(_lambda):\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n    return loss\n#=================","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Read and load data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(f\"{ROOT}/train.csv\")\ntrain.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\ntest = pd.read_csv(f\"{ROOT}/test.csv\")\n\nprint(\"add infos\")\nsub = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(test.drop('Weeks', axis=1), on=\"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Add reference to the concatenated data\ntrain['WHERE'] = 'train'\ntest['WHERE'] = 'val'\nsub['WHERE'] = 'test'\ndata = train.append([test, sub])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Add Variables","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = data.loc[data.Weeks == data.min_week]\nbase = base[['Patient','FVC']].copy()\nbase.columns = ['Patient','min_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\ndel base","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import LabelEncoder\nlb =LabelEncoder()\ndata[\"Sex\"] = lb.fit_transform(data[\"Sex\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FE = []\n\n\ndata['age'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\ndata['BASE'] = (data['min_FVC'] - data['min_FVC'].min() ) / ( data['min_FVC'].max() - data['min_FVC'].min() )\ndata['week'] = (data['base_week'] - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\ndata['percent'] = (data['Percent'] - data['Percent'].min() ) / ( data['Percent'].max() - data['Percent'].min() )\nFE += ['Sex','age','percent','week','BASE']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Use k-means Segmentation to create a cluster variable","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\nfrom sklearn.preprocessing import StandardScaler\nstd = StandardScaler()\n\nfeatures_scaled  =data[[\"age\",\"percent\",'Sex',\"BASE\"]]\n\nfrom sklearn.cluster import DBSCAN, KMeans\n## all data normalized\ninertia = []\nfor n in range(1 , 10):\n    algorithm = (KMeans(n_clusters = n ,init='k-means++', n_init = 100 ,max_iter=3000, \n                        tol=0.0001,  random_state= 111  , algorithm='elkan') )\n    algorithm.fit(features_scaled)\n    inertia.append(algorithm.inertia_)\nsns.set()\nplt.figure(1 , figsize = (15 ,6))\nplt.plot(np.arange(1 , 10) , inertia , 'o')\nplt.plot(np.arange(1 , 10) , inertia , '-' , alpha = 0.5)\nplt.xlabel('Number of Clusters') , plt.ylabel('Inertia')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Silhouette Score","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"range_n_clusters = [2,3,4,5,6,7,8,9]\nfrom sklearn.metrics import silhouette_score\n\nfor num_clusters in range_n_clusters:\n    \n    # intialise kmeans\n    kmeans = KMeans(n_clusters=num_clusters, max_iter=90)\n    kmeans.fit(features_scaled)\n    \n    cluster_labels = kmeans.labels_\n    \n    # silhouette score\n    silhouette_avg = silhouette_score(features_scaled, cluster_labels)\n    print(\"For n_clusters={0}, the silhouette score is {1}\".format(num_clusters, silhouette_avg))\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kmeans = KMeans(n_clusters=2, max_iter=3000)\nkmeans.fit(features_scaled)\n    \n\nfeatures_scaled[\"Cluster_Id\"] =kmeans.labels_\nfeatures_scaled[\"FVC\"] =data[\"FVC\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3-D Visualisation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import plotly.express as px\nfig = px.scatter_3d(features_scaled, x='percent', y='FVC', z='BASE',\n              color='Cluster_Id')\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(1 , figsize = (15 , 10))\nfor cluster in [0,1]:\n    plt.scatter(x = 'percent' , y = 'FVC' , data = features_scaled[features_scaled['Cluster_Id'] == cluster] ,\n                s = 10 , alpha = 1 , label = cluster)\nplt.xlabel('percent'), plt.ylabel('FVC') \nplt.title('FVC over percent')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"algorithm = (KMeans(n_clusters =2,init='k-means++', n_init = 10 ,max_iter=10000, \n                        tol=0.0001,  random_state= 111  , algorithm='elkan') )\nalgorithm.fit(features_scaled)\ndata['Cluster_Id'] = algorithm.labels_\n\ndel features_scaled","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Dummy SmokingStatus\n\nCOLS = ['SmokingStatus']\nfor col in COLS:\n    for mod in data[col].unique():\n        FE.append(mod)\n        data[mod] = (data[col] == mod).astype(int)\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = pd.get_dummies(data,columns=[\"Cluster_Id\"])\nfor col in (\"Cluster_Id_0\",\"Cluster_Id_1\"):\n    data[col] = data[col].astype(int)\n    \nFE+=[\"Cluster_Id_0\",\"Cluster_Id_1\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## DBSCAN clustering","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import metrics\n#varaible1 = std.fit_transform(data[\"FVC\"])\nfeatures_scaled  =data[[\"min_FVC\",\"percent\",\"BASE\",\"age\"]]\n\n\n## all data normalized\ninertia = []\n\nalgorithm = (DBSCAN(eps=0.6, metric='euclidean',\n                        metric_params=None, algorithm=\"auto\", leaf_size=20, p=None, n_jobs=None))\nalgorithm.fit(features_scaled)\nlabels = algorithm.labels_\ncore_samples_mask = np.zeros_like(labels, dtype=bool)\ncore_samples_mask[algorithm.core_sample_indices_] = True\nn_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)\nn_noise_ = list(labels).count(-1)# Printing the number of clusters and number of noise points (outliers)\nprint(\"Estimated number of clusters: %d\" % n_clusters_)\nprint(\"Estimated number of noise points: %d\" % n_noise_)# Evaluating the quality of clusters\ns = metrics.silhouette_score(features_scaled, algorithm.labels_)\nprint(f\"Silhouette Coefficient for the  Dataset Clusters: {s:.2f}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"algorithm = (DBSCAN(eps=0.5, min_samples=1, metric='euclidean',\n                        metric_params=None, algorithm=\"auto\", leaf_size=20, p=None, n_jobs=None))\nalgorithm.fit(features_scaled)\nlabels = algorithm.labels_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['DBSCAN_CLUSTER'] = algorithm.labels_\n\nfor mod in data[\"DBSCAN_CLUSTER\"].unique():\n    FE.append(mod)\n    data[mod] = (data[\"DBSCAN_CLUSTER\"] == mod).astype(int)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Resplit the data\ntrain = data.loc[data.WHERE=='train']\ntest = data.loc[data.WHERE=='val']\nsub = data.loc[data.WHERE=='test']\ndel data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y = train['FVC'].values\nz = train[FE].values\nze = sub[FE].values\npe = np.zeros((ze.shape[0], 3))\npred = np.zeros((z.shape[0], 3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Modeling","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras import layers\nfrom keras import regularizers\n\nimport keras\nC1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n#=============================#\n\ndef make_model(learning_rate):\n    z = L.Input((178,), name=\"Patient\")\n\n\n    x = L.Dense(200, activation=\"relu\",use_bias=True,kernel_regularizer=regularizers.l2(1e-7),\n               bias_regularizer=regularizers.l2(1e-7),\n               activity_regularizer=regularizers.l1(1e-7))(z)   \n    keras.layers.Dropout(0.1)\n\n\n    x = L.Dense(100, activation=\"relu\",use_bias=True,kernel_regularizer=regularizers.l2(1e-7),\n               bias_regularizer=regularizers.l2(1e-7),\n               activity_regularizer=regularizers.l1(1e-7))(x)   \n    keras.layers.Dropout(0.1)\n\n  \n    x = L.Dense(10, activation=\"relu\",use_bias=True)(x)   \n    keras.layers.Dropout(0.1)\n\n\n    p1 = L.Dense(3, activation=\"linear\", name=\"p1\",use_bias=True)(x)\n    p2 = L.Dense(3, activation=\"relu\", name=\"p2\",use_bias=True)(x)\n\n\n\n\n    \n    \n    preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\n    \n    model = M.Model(z, preds, name=\"CNN\")\n    #model.compile(loss=qloss, optimizer=\"adam\", metrics=[score])\n    model.compile(loss=mloss(0.8), optimizer=tf.keras.optimizers.Adam(lr=learning_rate, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False), metrics=[score])\n    return model\n\n''' x = L.Dense(90, activation=\"tanh\", name=\"d2v\",use_bias=True)(x)\n    x = L.Dense(45, activation=\"tanh\", name=\"d2f2\",use_bias=True)(x)\n    x = L.Dense(22, activation=\"tanh\", name=\"d2dfx\",use_bias=True)(x)\n    x = L.Dense(10, activation=\"tanh\", name=\"d4df5\",use_bias=True)(x) ''' \n\n'''kernel_regularizer=regularizers.l2(1e-7),\n               bias_regularizer=regularizers.l2(1e-7),\n               activity_regularizer=regularizers.l1(1e-7)'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Callbacks","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## EarlyStopping\nes_callback = keras.callbacks.EarlyStopping(monitor='val_loss', patience=15) \n\n## Keras API : Models often benefit from reducing the learning rate by a factor of 2-10 once learning stagnates.\n               #This callback monitors a quantity and if no improvement is seen for a 'patience' number of epochs,the learning rate is reduced.\nreduce_lr = tf.keras.callbacks.ReduceLROnPlateau(\n    monitor='val_loss', factor=0.5, patience=5, verbose=0, mode='auto',\n    min_delta=0.0001, cooldown=0, min_lr=0.0001\n)\n\n## Tensorboard\n#tensorboard_callback = tf.keras.callbacks.TensorBoard(log_dir=\"./logs\")\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NFOLD = 5\nkf = KFold(n_splits=NFOLD)\ncnt = 0\nBATCH_SIZE = 16\nfor tr_idx, val_idx in kf.split(z):\n    cnt += 1\n    print(f\"FOLD {cnt}\")\n    net = make_model(learning_rate=0.04)\n    net.fit(z[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=1000, \n            validation_data=(z[val_idx], y[val_idx]),callbacks=[es_callback,reduce_lr], verbose=0) #\n    print(\"train\", net.evaluate(z[tr_idx], y[tr_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"val\", net.evaluate(z[val_idx], y[val_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"predict val...\")\n    pred[val_idx] = net.predict(z[val_idx], batch_size=BATCH_SIZE, verbose=0)\n    print(\"predict test...\")\n    pe += net.predict(ze, batch_size=BATCH_SIZE, verbose=0) / NFOLD","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt = mean_absolute_error(y, pred[:, 1])\nunc = pred[:,2] - pred[:, 0]\nsigma_mean = np.mean(unc)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(pred[idxs, 0], label=\"q25\")\nplt.plot(pred[idxs, 1], label=\"q50\")\nplt.plot(pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(unc)\nplt.title(\"uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC1'] = 0.996  * pe[:, 1]\nsub['Confidence1'] = pe[:, 2] - pe[:, 0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.loc[~subm.FVC1.isnull(),'FVC'] = subm.loc[~subm.FVC1.isnull(),'FVC1']\nif sigma_mean<70:\n    subm['Confidence'] = sigma_opt\nelse:\n    subm.loc[~subm.FVC1.isnull(),'Confidence'] = subm.loc[~subm.FVC1.isnull(),'Confidence1']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"otest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nfor i in range(len(otest)):\n    subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'FVC'] = otest.FVC[i]\n    subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'Confidence'] = 0.1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}