{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install /kaggle/input/autograd/ -f ./ --no-index\n!pip install /kaggle/input/autogradgamma/ -f ./ --no-index\n!pip install /kaggle/input/lifelines/ -f ./ --no-index\n!pip install /kaggle/input/ngboost/ -f ./ --no-index","execution_count":null,"outputs":[]},{"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, GroupKFold\nfrom sklearn.metrics import make_scorer\nimport lightgbm as lgb\nimport typing as tp\nfrom scipy.optimize import lsq_linear\nfrom sklearn.linear_model import ElasticNet\nfrom tqdm.keras import TqdmCallback\n\nfrom ngboost.ngboost import NGBoost\nfrom ngboost.learners import default_tree_learner\nfrom ngboost.distns import Normal, LogNormal\nfrom ngboost.scores import MLE\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M","execution_count":null,"outputs":[]},{"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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = pd.read_csv(f\"{ROOT}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nchunk = 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(chunk.drop('Weeks', axis=1), on=\"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\ndata = tr.append([chunk, sub])","execution_count":null,"outputs":[]},{"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":"COLS = ['Sex','SmokingStatus'] #,'Age'\nFE = []\nfor col in COLS:\n    for mod in data[col].unique():\n        FE.append(mod)\n        data[mod] = (data[col] == mod).astype(int)\n#=================","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"percent__ = data[data['base_week'] == 0][['Patient', 'Percent']].values\npercent__ = dict(percent__)\n\npercent_list__ = []\nfor d in data.values:\n    percent_list__.append(percent__[d[0]])\n    \ndata['percent_base'] = percent_list__","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def preproc_sex_smoking_status(row):\n    if row['SmokingStatus'] == 'Currently smokes':\n        if row['Sex'] == 'Male':\n            return 1\n        if row['Sex'] == 'Female':\n            return 2\n    if row['SmokingStatus'] == 'Ex-smoker':\n        if row['Sex'] == 'Male':\n            return 3\n        if row['Sex'] == 'Female':\n            return 4\n    if row['SmokingStatus'] == 'Never smoked':\n        if row['Sex'] == 'Male':\n            return 5\n        if row['Sex'] == 'Female':\n            return 6","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['person_type'] = np.nan\n\nfor pid, row in enumerate(data.iloc):\n    data['person_type'][pid] = preproc_sex_smoking_status(row)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.join(pd.get_dummies(data['person_type'], prefix='person_type'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def to_normalize(cols):\n    normalized_cols = []\n    for idx, col in enumerate(cols):\n        new_col_name = col + '_normalized'\n        data[new_col_name] = (data[col] - data[col].min() ) / (data[col].max() - data[col].min())\n        normalized_cols.append(new_col_name)\n    return normalized_cols\n\ncols_to_normalize = ['Age', 'base_week', 'min_FVC', 'percent_base', 'Weeks']\n\nnormalized_cols = to_normalize(cols_to_normalize)\n\nFE2 = ['person_type_1.0', 'person_type_2.0', 'person_type_3.0',\n       'person_type_4.0', 'person_type_5.0', 'person_type_6.0'] + normalized_cols\n\nFE3 = ['person_type'] + cols_to_normalize","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = data.loc[data.WHERE=='train']\nchunk = data.loc[data.WHERE=='val']\nsub = data.loc[data.WHERE=='test']\n# del data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr.shape, chunk.shape, sub.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### BASELINE NN "},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_model(nh, delta):\n    \n    C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n    #=============================#\n    def 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    #============================#\n    def qloss(y_true, y_pred):\n        # Pinball loss for multiple quantiles\n        qs = [0.2, 0.5, 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    #=============================#\n    def 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    \n    z = L.Input((nh, 1), name=\"Patient\")\n    x = L.Dense(100, activation=\"relu\", name=\"d1\")(z)\n    conv = L.Conv1D(filters = 2, kernel_size = 2)(x)\n    x = L.Flatten()(conv)\n    x = L.Dense(100, activation=\"relu\", name=\"d2\")(x)\n    p1 = L.Dense(3, activation=\"relu\", name=\"p1\")(x)\n    p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\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(delta), optimizer=tf.keras.optimizers.Adam(lr=0.01, beta_1=0.9, beta_2=0.999, epsilon=1e-07, amsgrad=False), metrics=[score])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y = tr['FVC'].values\n# z = tr[FE].values\n# ze = sub[FE].values\n# nh = z.shape[1]\n\npe1 = np.zeros((sub[FE2].values.shape[0], 2))\npe2 = np.zeros((sub[FE2].values.shape[0], 2))\n\npred = np.zeros((tr[FE2].values.shape[0], 3))\n\nz2 = tr[FE2].values\n\nze2 = sub[FE2].values\nnh2 = ze2.shape[1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplace_log_likelihood(actual_fvc, predicted_fvc, confidence, return_values = False):\n    \"\"\"\n    Calculates the modified Laplace Log Likelihood score for this competition.\n    \"\"\"\n    sd_clipped = np.maximum(confidence, 70)\n    delta = np.minimum(np.abs(actual_fvc - predicted_fvc), 1000)\n    metric = - np.sqrt(2) * delta / sd_clipped - np.log(np.sqrt(2) * sd_clipped)\n\n    if return_values:\n        return metric\n    else:\n        return np.mean(metric)\n    \nscorer_210 = make_scorer(lambda y_true, y_pred: laplace_log_likelihood(y_true, y_pred, 210))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"net = make_model(nh2, 1)\nprint(net.summary())\nprint(net.count_params())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class OSICLossForLGBM:\n    \"\"\"\n    Custom Loss for LightGBM.\n    \n    * Objective: return grad & hess of NLL of gaussian\"\n    * Evaluation: return competition metric\n    \"\"\"\n    \n    def __init__(self, epsilon: float=1) -> None:\n        \"\"\"Initialize.\"\"\"\n        self.name = \"osic_loss\"\n        self.n_class = 2  # FVC & Confidence\n        self.epsilon = epsilon\n    \n    def __call__(self, preds: np.ndarray, labels: np.ndarray, weight: tp.Optional[np.ndarray]=None) -> float:\n        \"\"\"Calc loss.\"\"\"\n        sigma_clip = np.maximum(preds[:, 1], 70)\n        Delta = np.minimum(np.abs(preds[:, 0] - labels), 1000)\n        loss_by_sample = - np.sqrt(2) * Delta / sigma_clip - np.log(np.sqrt(2) * sigma_clip)\n        loss = np.average(loss_by_sample, weight)\n        \n        return loss\n    \n    def _calc_grad_and_hess(\n        self, preds: np.ndarray, labels: np.ndarray, weight: tp.Optional[np.ndarray]=None\n    ) -> tp.Tuple[np.ndarray]:\n        \"\"\"Calc Grad and Hess\"\"\"\n        mu = preds[:, 0]\n        sigma = preds[:, 1]\n        \n        sigma_t = np.log(1 + np.exp(sigma))\n        grad_sigma_t = 1 / (1 + np.exp(- sigma))\n        hess_sigma_t = grad_sigma_t * (1 - grad_sigma_t)\n        \n        grad = np.zeros_like(preds)\n        hess = np.zeros_like(preds)\n        grad[:, 0] = - (labels - mu) / sigma_t ** 2\n        hess[:, 0] = 1 / sigma_t ** 2\n        \n        tmp = ((labels - mu) / sigma_t) ** 2\n        grad[:, 1] = 1 / sigma_t * (1 - tmp) * grad_sigma_t\n        hess[:, 1] = (\n            - 1 / sigma_t ** 2 * (1 - 3 * tmp) * grad_sigma_t ** 2\n            + 1 / sigma_t * (1 - tmp) * hess_sigma_t\n        )\n        if weight is not None:\n            grad = grad * weight[:, None]\n            hess = hess * weight[:, None]\n        return grad, hess\n    \n    def return_loss(self, preds: np.ndarray, data: lgb.Dataset) -> tp.Tuple[str, float, bool]:\n        \"\"\"Return Loss for lightgbm\"\"\"\n        labels = data.get_label()\n        weight = data.get_weight()\n        n_example = len(labels)\n        \n        # # reshape preds: (n_class * n_example,) => (n_class, n_example) =>  (n_example, n_class)\n        preds = preds.reshape(self.n_class, n_example).T\n        # # calc loss\n        loss = self(preds, labels, weight)\n        \n        return self.name, loss, True\n    \n    def return_grad_and_hess(self, preds: np.ndarray, data: lgb.Dataset) -> tp.Tuple[np.ndarray]:\n        \"\"\"Return Grad and Hess for lightgbm\"\"\"\n        labels = data.get_label()\n        weight = data.get_weight()\n        n_example = len(labels)\n        \n        # # reshape preds: (n_class * n_example,) => (n_class, n_example) =>  (n_example, n_class)\n        preds = preds.reshape(self.n_class, n_example).T\n        # # calc grad and hess.\n        grad, hess =  self._calc_grad_and_hess(preds, labels, weight)\n\n        # # reshape grad, hess: (n_example, n_class) => (n_class, n_example) => (n_class * n_example,) \n        grad = grad.T.reshape(n_example * self.n_class)\n        hess = hess.T.reshape(n_example * self.n_class)\n        \n        return grad, hess\n    \ndef find_optimal_solution(preds, targets):\n    A = np.array(preds).T\n    res = lsq_linear(A, targets, lsq_solver='exact', method='trf', tol=1e-5, verbose=2)\n    return A.dot(res.x), res.x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cnt = 0\ndelta = 0.4\nEPOCHS = 2000\nBATCH_SIZE = 512\nMODELS = 1\nNFOLD = 5\ngkf = GroupKFold(n_splits = NFOLD)\nNFOLD_MODELS = NFOLD * MODELS\nval_scores = []\nlr_scheduler = tf.keras.callbacks.ReduceLROnPlateau(\n                monitor=\"val_loss\",\n                factor=0.05,\n                patience=100,\n                verbose=0,\n                mode=\"auto\",\n                min_delta=0.0001,\n                cooldown=0,\n                min_lr=0.0001)\npred = np.zeros((z2.shape[0], 3))\ncategorical_features = ['person_type']\nfold_n = 0\n\nfor tr_idx, val_idx in gkf.split(z2, y, data[data.WHERE=='train'].Patient.values):\n    fold_n += 1\n    seed_everything(42)\n\n    ngb = NGBoost(Base = default_tree_learner, Dist = Normal, Score=MLE, natural_gradient = True, verbose = False).fit(tr[FE3].iloc[tr_idx].values, y[tr_idx])\n    ngb_val_pred = ngb.predict(tr[FE3].iloc[val_idx].values)\n    ngb_val_dists = ngb.pred_dist(tr[FE3].iloc[val_idx].values, 1)\n    a, b = ngb_val_dists.dist.interval(0.1)\n    ngb_val_conf = b - a\n    \n    ngb_subm_pred = ngb.predict(sub[FE3].values)\n    ngb_subm_dists = ngb.pred_dist(sub[FE3].values, 1)\n    a, b = ngb_subm_dists.dist.interval(0.1)\n    ngb_subm_conf = b - a\n    \n    net = make_model(nh2, delta)\n    net.fit(z2[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, \n            validation_data=(z2[val_idx], y[val_idx]), verbose=0, callbacks = [TqdmCallback(verbose = 0), lr_scheduler])\n    pred[val_idx] = net.predict(z2[val_idx], batch_size=BATCH_SIZE, verbose=0)\n\n    nn_val_pred = pred[val_idx]\n    nn_val_pred_conf = nn_val_pred[:, 2] - nn_val_pred[:, 0]\n    nn_val_loss = laplace_log_likelihood(nn_val_pred[:, 1], y[val_idx], nn_val_pred_conf)\n        \n    train_data = lgb.Dataset(tr[FE3].iloc[tr_idx], label = y[tr_idx], categorical_feature = categorical_features)\n    test_data = lgb.Dataset(tr[FE3].iloc[val_idx], label = y[val_idx], categorical_feature = categorical_features)\n    \n    lgb_model_param = {\n    'num_class': 2,\n    'metric': 'None',\n    'boosting_type': 'gbdt',\n    'learning_rate': 0.05,\n    'seed': 42,\n    \"subsample\": 0.4,\n    \"subsample_freq\": 1,\n    'max_depth': 1,\n    'verbosity': 0    \n    }\n    \n    lgb_fit_param = {\n        \"num_boost_round\": 500000,\n        \"verbose_eval\": 5000,\n        \"early_stopping_rounds\": 500,\n    }\n    \n    loss = OSICLossForLGBM()\n    \n    model = lgb.train(lgb_model_param,\n                       train_data,\n                       valid_sets = test_data,\n                       **lgb_fit_param, \n                      fobj=loss.return_grad_and_hess,\n                      feval=loss.return_loss)\n    lgb_val_predict = model.predict(tr[FE3].iloc[val_idx].values)\n    lgb_val_loss = laplace_log_likelihood(lgb_val_predict[:, 0], y[val_idx], nn_val_pred_conf)\n    lgb_val_loss_conf = laplace_log_likelihood(lgb_val_predict[:, 0], y[val_idx], lgb_val_predict[:, 1])\n    predicted = model.predict(sub[FE3].values)\n    \n    el_fitted = ElasticNet(alpha=0.3, l1_ratio = 0.8).fit(z2[tr_idx], y[tr_idx])\n    el = el_fitted.predict(z2[val_idx])\n    el_subm = el_fitted.predict(ze2)\n    \n    print(f\"Loss Keras #{fold_n}: {nn_val_loss}\")\n    print(f\"Loss LGBM #{fold_n}: {lgb_val_loss}, {lgb_val_loss_conf}\")\n    \n    optimal = find_optimal_solution([nn_val_pred[:, 1], lgb_val_predict[:, 0], el, ngb_val_pred], y[val_idx])\n    optimal_loss = laplace_log_likelihood(optimal[0], y[val_idx], nn_val_pred_conf)\n    optimal_conf = laplace_log_likelihood(optimal[0], y[val_idx], lgb_val_predict[:, 1])\n    optimal_conf_updated = laplace_log_likelihood(optimal[0], y[val_idx], np.mean([lgb_val_predict[:, 1], nn_val_pred_conf], axis = 0))\n    optimal_conf_updated_ngb = laplace_log_likelihood(optimal[0], y[val_idx], np.mean([lgb_val_predict[:, 1], nn_val_pred_conf, ngb_val_conf], axis = 0))\n    print(\"Optimal loss:\", optimal_loss, optimal_conf, optimal_conf_updated, optimal_conf_updated_ngb)\n    print(\"Optimal coefficients:\", optimal[1])\n    \n    subm_predict = net.predict(ze2, batch_size=BATCH_SIZE, verbose=0)\n    subm_fvc = np.array([subm_predict[:, 1], predicted[:, 0], el_subm, ngb_subm_pred]).T.dot(optimal[1])\n    subm_conf = np.mean([subm_predict[:, 2] - subm_predict[:, 0], predicted[:, 1], ngb_subm_conf], axis = 0)\n    subm_predict = np.array([subm_fvc, subm_conf]).T\n  \n    pe1 += subm_predict / NFOLD_MODELS\n    \n    print()","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":{},"cell_type":"markdown","source":"### PREDICTION"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC1'] = pe1[:, 0]\nconf = pe1[:, 1]\nsub['Confidence1'] = conf # np.minimum(conf, conf.mean())","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()].head(10)","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":"subm.describe()","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":{"trusted":true},"cell_type":"code","source":"subm","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}