{"cells":[{"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 tensorflow as tf\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport tensorflow_probability as tfp\nimport seaborn as sns\nimport pydicom\nimport os\nimport re\nimport time\nimport gzip\nfrom tqdm import tqdm\nimport seaborn as sns\nimport sklearn.preprocessing\ntfd = tfp.distributions\ntfb = tfp. bijectors\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn.model_selection import train_test_split\nimport ct_scan_processing2_module as ct\nimport osic_utils\nfrom tqdm import tqdm_notebook as tqdm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Determine patient ids from train "},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_stats = pd.read_csv('../usr/lib/ct_scan_processing2/lung_statistics.csv')\ncols = lung_stats.columns[~pd.Series(lung_stats.columns).isin(['PatientId'])]\nrow, col_ = (3,2)\nfig, ax = plt.subplots(row, col_, figsize = (15,9))\nplt.subplots_adjust(bottom = -0.4)\nfor i, col in enumerate(cols):\n    ridx = int(i/col_)\n    cidx = i - col_ * ridx      \n    sns.distplot(lung_stats[col], ax = ax[ridx][cidx])\n    ax[ridx][cidx].set_xlabel('')\n    ax[ridx][cidx].set_title(col)\nfig.delaxes(ax[row - 1][col_ - 1])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_stats.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_dcm_dict = {}\nroot_dir = \"/kaggle/input/osic-pulmonary-fibrosis-progression/train/\"\n\nremove_ids = []\nfor dirname, _, filenames in os.walk(root_dir):\n    if 'ID' in dirname:\n        dirname_ = dirname.replace(root_dir, \"\")\n        # check if the id is among the ids to be removed\n        s = [i for i in remove_ids if dirname_ == i]\n        if len(s) == 0:            \n            patient_dcm_dict[dirname_] = filenames\n        \npatient_dcm_dict = {k: i for i, k in enumerate(sorted(patient_dcm_dict.keys()))}\nprint(\"A total of\", len(patient_dcm_dict.keys()), \"patients\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_dcm_dict = {k: i for i, k in enumerate(sorted(patient_dcm_dict.keys()))}\n{k: v for i, (k, v) in enumerate(patient_dcm_dict.items()) if i < 5}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sort the lung statistics dataset by the patient dictionary"},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_stats = lung_stats.set_index('PatientId') \\\n.loc[list(patient_dcm_dict.keys())] \\\n.reset_index()\n\nlung_stats['lung_vol'] = lung_stats['lung_vol']/10000\nlung_stats.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We now create a tf dataset for easy batching and shuffling."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"markdown","source":"## Preprocess structured data"},{"metadata":{},"cell_type":"markdown","source":"Outlined below is a pipeline for processing the data ranging from structured to CT scan images."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = pd.read_csv(\"/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv\")\ntest_data = pd.read_csv(\"/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv\")\n\n### fit encoders before running preprocessing data ###\n# one-hot encoder\noh_encoder = OneHotEncoder(sparse = False, dtype = 'int32')\noh_encoder.fit(train_data['SmokingStatus'].values.reshape(-1, 1))\n\n### Standard scalers ####\n# AGE\nage_scaler = StandardScaler()\nage_scaler.fit(train_data[['Patient', 'Age']].drop_duplicates()['Age'] \\\n               .values.reshape(-1, 1))\n# WEEKS\nweek_scaler = StandardScaler()\nweek_scaler.fit(train_data[['Patient', 'Weeks']].drop_duplicates()['Weeks'] \\\n               .values.reshape(-1,1))\n# PERCENTAGE\npct_scaler = StandardScaler()\npct_scaler.fit(train_data['Percent'].values.reshape(-1,1))\n# BASELINE FVC\nbaseline_scaler = StandardScaler()\nbaseline_scaler.fit(train_data['FVC'].values.reshape(-1,1))\n\nencoder_dict = {\n    'Age': age_scaler,\n    'Weeks': week_scaler,\n    'pct_bsl': pct_scaler,\n    'baseline': baseline_scaler,\n    'SmokingStatus': oh_encoder\n}\nX, pid, lbl = osic_utils.preprocess_structured_data(\n    train_data, \n    patient_dcm_dict, \n    encoder_dict)\nprint(X.shape)\nprint(pid.shape)\nprint(lbl.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Concatenate the lung statistics with the structured data features"},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_stats_train = tf.gather(lung_stats.iloc[:,1:].values, \n                             tf.constant(pid, dtype = tf.int32), axis = 0)\nX = tf.concat([X,lung_stats_train], axis = 1)\nX","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"holdout = osic_utils.transform_holdout_data(test_data)\n\n# there are no actual labels. \n#lbl_holdout is just there as placeholder for consistency with train tf dataset generation\npatient_dcm_holdout = {k: i for i, k in enumerate(holdout.Patient.unique())}\nX_holdout, pid_holdout, lbl_holdout = osic_utils.preprocess_structured_data(holdout, \n                                                         patient_dcm_holdout,\n                                                         encoder_dict)\nprint(X_holdout.shape)\nprint(pid_holdout.shape)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create tf dataset for train\ndataset = osic_utils.DatasetGen(patient_dcm_dict, split = None, seed = 300, \n                         batch_size = 1, \n                         root_dir = '/kaggle/input/osic-pulmonary-fibrosis-progression/train/',\n                         shuffle = False)\n\n# create tf dataset for test\n\"\"\"\nholdout_ds = osic_utils.DatasetGen(patient_dcm_holdout, \n                                   split = None,\n                                   seed = 300, batch_size =1,\n                                   root_dir = '/kaggle/input/osic-pulmonary-fibrosis-progression/test/',\n                                   struct_data = X_holdout, label = lbl_holdout, \n                                   id_tensor = pid_holdout,\n                                   return_stats_only = True,\n                                  shuffle = False)\n\"\"\"\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_joint_dist_coroutine(num_ids, pids, X):\n    def model():\n        X_cst = tf.cast(X, tf.float32)\n        # random intercepts\n        Root = tfd.JointDistributionCoroutine.Root\n        patient_scale = yield Root(tfd.HalfCauchy(loc = 0, scale = 5.))\n        intercept = yield Root(tfd.Normal(loc = 0, scale = 10.))\n        patient_prior = yield tfd.MultivariateNormalDiag(loc = tf.zeros(num_ids), \n                                                       scale_identity_multiplier = patient_scale)\n        int_resp = tf.gather(patient_prior, pids, axis = -1) + intercept[...,tf.newaxis]\n        \n        # random slopes for week var\n        beta_week = yield Root(tfd.Normal(loc = 0, scale = 10.))\n        beta_week_scale = yield Root(tfd.HalfCauchy(loc = 0, scale = 5.))\n        beta_week_prior = yield tfd.MultivariateNormalDiag(loc = tf.zeros(num_ids),\n                                                          scale_identity_multiplier = beta_week_scale)\n        bw_resp = (tf.gather(beta_week_prior, pids, axis = -1)  + beta_week[..., tf.newaxis]) * X_cst[:,0]\n        \n        # other variables\n        betas = yield Root(tfd.MultivariateNormalDiag(loc = tf.zeros(tf.shape(X_cst)[1] - 1), \n                                                scale_identity_multiplier = 10.))\n        other_vars_resp = tf.tensordot(betas, tf.transpose(X_cst[:,1:]), axes = 1)    \n        total_response = int_resp + bw_resp + other_vars_resp\n        \n        # response scale\n        resp_scale = yield Root(tfd.HalfCauchy(loc =0, scale = 5.))\n        \n        yield tfd.Normal(loc = total_response,  scale = resp_scale[...,tf.newaxis])\n        \n    return tfd.JointDistributionCoroutineAutoBatched(model)\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Try sampling from the joint distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"dist = make_joint_dist_coroutine(176, tf.constant(pid, tf.int32), X)\n[i.shape for i in dist.sample(2)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Stochastic Variational Inference on hierarchical linear model"},{"metadata":{"trusted":true,"_kg_hide-input":false,"_kg_hide-output":false},"cell_type":"code","source":"\"\"\"\n# define parameters to be optimized\n_init_loc = lambda name, shape = (): tf.Variable(\n    tf.random.uniform(shape, name = name,  minval=-2., maxval=2.))\n\n_init_scale =lambda name, shape = (): tfp.util.TransformedVariable(\n    initial_value=tf.random.uniform(shape, minval=0.01, maxval=1.),\n    bijector=tfb.Softplus(), name = name)\n\n### SET PARAMETERS HERE ####\nnum_epochs = 1\nprint_every = 5\nnum_ids = len(patient_dcm_dict.keys())\nnum_draws = 2\n\n# intercept & random intercepts\na0 = _init_loc(name = \"alpha0\")\na0_sigma = _init_scale(name = \"alpha0_sigma\")\npatient_scale = _init_loc(name = \"patient_scale\")\npatient_scale_sigma = _init_scale(name = \"patient_scale_sigma\")\na = _init_loc(name = \"alphas\", shape = [num_ids])\na_sigma = _init_scale(name = \"alphas_sigma\", shape = [num_ids])\n\nint_vars = [patient_scale, patient_scale_sigma.trainable_variables[0],\n           a0, a0_sigma.trainable_variables[0],\n           a, a_sigma.trainable_variables[0]]\n\n# random slope (week variable)\nbw = _init_loc(name = \"bw\")\nbw_sigma = _init_scale(name = \"bw_sigma\")\nbws_scale = _init_loc(name = \"bw_scale\")\nbws_scale_sigma = _init_scale(name = \"bws_scale_sigma\")\nbws = _init_loc(name = \"bws\", shape = [num_ids])\nbws_sigma = _init_scale(name = \"bws_sigma\", shape = [num_ids])\n\nrnd_slope_vars  = [bw, bw_sigma.trainable_variables[0], \n                   bws_scale, bws_scale_sigma.trainable_variables[0],\n                   bws, bws_sigma.trainable_variables[0]]\n\n# other vars\nb = _init_loc(name = \"betas\", shape = [11])\nb_sigma = _init_scale(name = \"betas_sigma\", shape = [11])\n\nother_vars = [b, b_sigma.trainable_variables[0]]\n\n# response scale\nresp_scale = _init_loc(name = \"resp_scale\")\nresp_scale_sigma = _init_scale(name = \"resp_scale_sigma\")\nresponse_scale_vars = [resp_scale, resp_scale_sigma.trainable_variables[0]]\n\ntrainable_vars = int_vars + rnd_slope_vars + other_vars + response_scale_vars \n\noptimizer = tf.optimizers.Adam(learning_rate = 0.001)\nepoch_loss = []\nstart = time.time()\nfor e in range(num_epochs):\n    batch_loss = []\n    i = 0\n    for ids, X, lbl, idx, img_stats in dataset.train:\n        with tf.GradientTape() as tape:      \n            lbl = tf.cast(lbl/1000, tf.float32) #reduce the magnitude of the labels\n            shp = tf.shape(ids)\n            # concatenate structured data and extracted image features\n            X = tf.concat([tf.cast(X, tf.float32), tf.cast(img_stats, tf.float32)], axis = 1)\n            jd = make_joint_dist_coroutine(shp.numpy()[0], idx, X)\n        \n            # create surrogate posterior dynamically\n            def variational_model_fn():\n                return tfd.JointDistributionSequentialAutoBatched([\n                  tfb.Softplus()(tfd.Normal(patient_scale, patient_scale_sigma)),  # scale_prior\n                  tfd.Normal(a0, a0_sigma),                                        # intercept\n                  tfd.Normal(tf.gather(a, ids), tf.gather(a_sigma, ids)),          # patient prior\n                  tfd.Normal(bw, bw_sigma),                                        # week random slope\n                  tfb.Softplus()(tfd.Normal(bws_scale, bws_scale_sigma)),          # bw scale prior\n                  tfd.Normal(tf.gather(bws, ids), tf.gather(bws_sigma, ids)),      # patient by week prior\n                  tfd.Normal(b, b_sigma),                                          # other vars prior\n                  tfb.Softplus()(tfd.Normal(resp_scale, resp_scale_sigma)),        # response scale\n                    \n                ])\n\n            # create the surrogate posterior\n            surrogate_pos = variational_model_fn()\n            # calculate losses\n            loss = tfp.vi.monte_carlo_variational_loss(\n                lambda *args: jd.log_prob(*args, lbl),\n                surrogate_pos,\n                sample_size = num_draws,\n                use_reparameterization = True\n            )\n            \n        # compute gradients \n        grad = tape.gradient(loss, trainable_vars)        \n        optimizer.apply_gradients(zip(grad, trainable_vars))\n        print(\"epoch {}, batch {} loss: {}\".format(e, i, loss))\n        i += 1\n        batch_loss.append(loss)\n        \n    epoch_loss.append(np.mean(batch_loss))\n    if e % print_every == 0:\n        print(\"epoch {} loss: {}\".format(e, np.mean(batch_loss)))\nend = time.time()\n\nprint(\"Total processing time:\", str((end - start)/60), \"mins\")\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_ids = len(patient_dcm_dict.keys())\njd = make_joint_dist_coroutine(num_ids, pid, X)\ndef target_log_prob_fn(*args):\n    return jd.log_prob(*args, lbl/1000)\n    \n_init_loc = lambda shape=(): tf.Variable(\n    tf.random.uniform(shape, minval=-2., maxval=2.))\n_init_scale = lambda shape=(): tfp.util.TransformedVariable(\n    initial_value=tf.random.uniform(shape, minval=0.01, maxval=1.),\n    bijector=tfb.Softplus())\n\nnum_other_vars = X.shape[1] - 1\n\nsurrogate_posterior = tfd.JointDistributionSequentialAutoBatched([\n                  tfb.Softplus()(tfd.Normal(_init_loc(), _init_scale())),                       # scale_prior\n                  tfd.Normal(_init_loc(), _init_scale()),                                       # intercept\n                  tfd.Normal(_init_loc(shape = [num_ids]), _init_scale(shape = [num_ids])),      # patient prior\n                  tfd.Normal(_init_loc(), _init_scale()),                                       # week random slope\n                  tfb.Softplus()(tfd.Normal(_init_loc(), _init_scale())),                       # week scale prior\n                  tfd.Normal(_init_loc(shape = [num_ids]), _init_scale(shape = [num_ids])),      # patient by week prior\n                  tfd.Normal(_init_loc(shape = [num_other_vars]), \n                             _init_scale(shape = [num_other_vars])),                            # other vars prior\n                  tfb.Softplus()(tfd.Normal(_init_loc(), _init_scale())),                       # response scale\n                    \n                ])\n\n\noptimizer = tf.optimizers.Adam(learning_rate=1e-4)\n\nstart = time.time()\nlosses = tfp.vi.fit_surrogate_posterior(\n    target_log_prob_fn,  \n    surrogate_posterior,\n    optimizer=optimizer,\n    num_steps= 100000, \n    seed=42,\n    sample_size= 500)\nend = time.time()\n\nprint(\"processing time:\", str((end - start)/60), \"min\")\n\n(scale_prior_, \n intercept_,  \n patient_weights,\n week_slope,\n week_scale_prior,\n week_patient_weights,\n other_vars_weights,\n response_scale), _ = surrogate_posterior.sample_distributions()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Model Checking"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('intercept:', intercept_.mean())\nprint('week weight:', week_slope.mean())\nprint('other var weights', other_vars_weights.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize = (10, 6))\nplt.plot(losses)\nplt.title('ELBO loss across epochs', size = 15)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"def predict_from_params(X):\n    other_vars_weights_ = other_vars_weights.mean().numpy().reshape(-1, 1)\n    intercept = intercept_.mean().numpy().reshape(-1, 1)\n    week_slope_ = week_slope.mean().numpy().reshape(-1, 1)\n    expected_val = (np.matmul(X[:, 1:], other_vars_weights_)\n                    + week_slope_ * X[:, :1]\n                    + intercept)\n    return expected_val","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Check train fit"},{"metadata":{"trusted":true},"cell_type":"code","source":"min_FVC = np.min(lbl/1000)\nmax_FVC = np.max(lbl/1000)\nplt.subplots(1,1, figsize = (8,6))\nplt.scatter(np.squeeze(predict_from_params(X)), lbl/1000, alpha =0.3)\nplt.xlabel('Predicted', size = 10)\nplt.ylabel('Actual', size = 10)\nplt.plot([min_FVC, max_FVC], [min_FVC, max_FVC], \"--\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds_holdout = []\nfor i, (ids, idx) in enumerate(patient_dcm_holdout.items()):\n    path = \"/kaggle/input/osic-pulmonary-fibrosis-progression/test/\" + ids\n    stats = ct.image_processing_pipeline(path, return_stats_only = True)\n    stats[0,0] = stats[0,0]/1000\n    # get corresponding feature set\n    idxs = np.where(pid_holdout == idx)[0]\n    X_h = np.take(X_holdout, idxs, axis = 0)   \n    pred = predict_from_params(np.concatenate([X_h, np.repeat(stats, idxs.shape[0], axis=0)], axis=1))\n    preds_holdout.append(pred)\n    \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"holdout['FVC'] = np.concatenate(preds_holdout, axis = 0) * 1000\nholdout['Confidence'] = np.mean(response_scale.sample(50000).numpy()) * 1000\nholdout['Patient_Week'] = holdout[['Patient', 'Weeks', 'Weeks_base']] \\\n.apply(lambda x: '{}_{}'.format(x['Patient'], x['Weeks'] + x['Weeks_base']), axis = 1)\nholdout\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = holdout[['Patient_Week',  'FVC', 'Confidence']]\nsubmission.to_csv(\"submission.csv\", index = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","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}