{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# import os\n# i = 0\n# for dirname, _, filenames in os.walk('/kaggle/input'):  \n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n#         i += 1\n#         if i > 10:\n#             break","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import pydicom\nimport os\nfrom os import listdir\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"im_path = \"../input/osic-pulmonary-fibrosis-progressiont/\"\ntrain_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\nprint('Training data shape: ', train_df.shape)\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Get tabular data as features and observations","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.get_dummies(train_df, columns=['Sex'])\ntrain_df = pd.get_dummies(train_df, columns=['SmokingStatus'])\ntrain_df = train_df.rename(columns={\"Sex_Female\": \"Female\", \n                                    \"Sex_Male\": \"Male\",\n                                    \"SmokingStatus_Currently smokes\": \"CurrentlySmokes\",\n                                    \"SmokingStatus_Ex-smoker\": \"ExSmoker\",\n                                    \"SmokingStatus_Never smoked\": \"NeverSmoked\"})\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train_df.drop(['Patient','FVC'], axis=1)\ny = train_df['FVC']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test_df = pd.get_dummies(test_df, columns=['Sex'])\n# test_df = pd.get_dummies(test_df, columns=['SmokingStatus'])\n# test_df = test_df.rename(columns={  \"Sex_Male\": \"Male\",\n#                                     \"SmokingStatus_Ex-smoker\": \"ExSmoker\",\n#                                     \"SmokingStatus_Never smoked\": \"NeverSmoked\"})\n# test_df.insert(5, 'Female', np.zeros(5))\n# test_df.insert(7,'CurrentlySmokes',np.zeros(5))\n# test_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# XGBoost Regression","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Train-Test Splitting","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Splite data into training and testing\nfrom sklearn import model_selection\n\n# Reserve 20% for testing\nX_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.2)\n\nprint('training data has ' + str(X_train.shape[0]) + \n      ' observation with ' + str(X_train.shape[1]) + ' features')\nprint('test data has ' + str(X_test.shape[0]) + \n      ' observation with ' + str(X_test.shape[1]) + ' features')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Data Scaling by Standardization","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# standardization (x-mean)/std\n# normalization (x-x_min)/(x_max-x_min) ->[0,1]\n\nfrom sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler()\nscaler.fit(X_train)\nX_train = scaler.transform(X_train)\nX_test = scaler.transform(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import xgboost as xgb\nfrom xgboost import XGBRegressor\nregr_XGB = XGBRegressor()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"regr_XGB.fit(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import model_selection\nfrom sklearn.model_selection import GridSearchCV\ncv_score = model_selection.cross_val_score(regr_XGB, X_train, y_train, cv=5)\nprint(cv_score)\n\n# Possible hyperparamter options for XGBoost\n# Choose the number of trees, max depth and other\nparameters = {'max_depth': np.arange(5, 10),\n              'colsample_bytree': np.arange(0.6,1,0.1),\n              'subsample': np.arange(0.6,1,0.1),\n              'eta': np.logspace(-2, 0, 10)\n              }\nGrid_XGB = GridSearchCV(regr_XGB,parameters, cv=5)\nGrid_XGB.fit(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"best_XGB_model = Grid_XGB.best_estimator_\nbest_XGB_model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"regr_XGB_opt = XGBRegressor(base_score=0.5, booster='gbtree', colsample_bylevel=1,\n             colsample_bynode=1, colsample_bytree=0.8999999999999999, eta=0.01,\n             gamma=0, gpu_id=-1, importance_type='gain',\n             interaction_constraints='', learning_rate=0.300000012,\n             max_delta_step=0, max_depth=6, min_child_weight=1, missing=None,\n             monotone_constraints='()', n_estimators=100, n_jobs=0,\n             num_parallel_tree=1, random_state=0, reg_alpha=0, reg_lambda=1,\n             scale_pos_weight=1, subsample=0.7999999999999999,\n             tree_method='exact', validate_parameters=1, verbosity=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"regr_XGB_opt.fit(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(5,5))\nplt.scatter(y_test, y_pred,color = 'r', alpha = 0.3)\nplt.plot([min(y_test),max(y_test)],[min(y_test),max(y_test)], color = 'k')\nplt.xlabel('FVC$_{\\mathrm{test}}$')\nplt.ylabel('FVC$_{\\mathrm{pred}}$')\nplt.rcParams.update({'font.size': 22})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_squared_error\nmse = np.sqrt(mean_squared_error(y_test, y_pred))\nprint(\"RMSE: %f\" % (mse**0.5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dmatrix = xgb.DMatrix(data=X,label=y)\nparams = {\"objective\":\"reg:squarederror\",'colsample_bytree': 0.3,'learning_rate': 0.1,\n                'max_depth': 5, 'alpha': 10}\ncv_results = xgb.cv(dtrain=data_dmatrix, params=params, nfold=10,\n                    num_boost_round=200,early_stopping_rounds=10,metrics=\"rmse\", as_pandas=True, seed=123)\nprint((cv_results[\"test-rmse-mean\"]).tail(1))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"importances = regr_XGB.feature_importances_\nindices = np.argsort(importances)[::-1]\n# Print the feature ranking\nprint(\"Feature importance ranking by XGBoost Model:\")\nfor ind in range(X.shape[1]):\n    print (\"%s : %.4f\" %(X.columns[indices[ind]],importances[indices[ind]]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler()\nscaler.fit(X)\nX = scaler.transform(X)","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":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\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 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\ndef make_model():\n    z = L.Input((8,), name=\"Patient\")\n    x = L.Dense(128, activation=\"relu\", name=\"d1\")(z)\n    x = L.Dense(512, activation=\"relu\", name=\"d2\")(x)\n    #x = L.Dense(100, activation=\"relu\", name=\"d3\")(x)\n    p1 = L.Dense(3, activation=\"linear\", 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(0.8), optimizer=\"adam\", metrics=[score])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"net = make_model()\nprint(net.summary())\nprint(net.count_params())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\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.drop('Patient_Week', axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\nNFOLD = 5\nkf = KFold(n_splits=NFOLD)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ncnt = 0\npred = np.zeros((X.shape[0], 3))\nfor tr_idx, val_idx in kf.split(X):\n    cnt += 1\n    print(f\"FOLD {cnt}\")\n    net = make_model()\n    \n    net.fit(X[tr_idx], y[tr_idx], batch_size=50, epochs=500, \n            validation_data=(X[val_idx], y[val_idx]), verbose=0) \n    print(\"train\", net.evaluate(X[tr_idx], y[tr_idx], verbose=0, batch_size=50))\n    print(\"val\", net.evaluate(X[val_idx], y[val_idx], verbose=0, batch_size=50))\n    print(\"predict val...\")\n    pred[val_idx] = net.predict(X[val_idx], batch_size=50, verbose=0)\n    print(\"predict test...\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_absolute_error\nsigma_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":"plt.figure(figsize=(5,5))\nplt.scatter(y,pred[:,1], color = 'r', alpha = 0.3)\nplt.plot([1000,6000],[1000,6000], color = 'k')\nplt.xlabel('FVC$_{\\mathrm{test}}$')\nplt.ylabel('FVC$_{\\mathrm{pred}}$')\nplt.rcParams.update({'font.size': 22})\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mse = np.sqrt(mean_squared_error(y, pred[:,1]))\nprint(\"RMSE: %f\" % (mse**0.5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nidxs = np.random.randint(0, y.shape[0], 50)\nplt.figure(figsize=(10,10))\nplt.plot(y.values[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":{},"cell_type":"markdown","source":"## Fitting Model Selection","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def lin_decay(w, a, b):\n    return a * w + b\n\ndef exp_decay(w, r, b, c):\n    return c * np.exp(-r * w) + b\n\ndef log_decay(w, r, b, c):\n    return b / (1 + c * np.exp(-r * w))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\npatient_ids = os.listdir('/kaggle/input/osic-pulmonary-fibrosis-progression/train')\n\nFVC_list = []\nweek_list = []\npct_list = []\n\npars_lin, pcov_lin = [], []\npars_exp, pcov_exp = [], []\npars_log, pcov_log = [], []\nRSS_lin = RSS_exp = RSS_log = 0\n\n#fit_values = []\n#columns = ['Patient', 'r', 'b']\n    \nfor i in range(len(patient_ids)):\n    pid = patient_ids[i]\n    week = train_df.loc[train_df['Patient'] == pid]['Weeks']\n    week = week - min(week)\n    week_list.append(week)\n    FVC = train_df.loc[train_df['Patient'] == pid]['FVC']\n    FVC_list.append(FVC)\n    pct = train_df.loc[train_df['Patient'] == pid]['Percent']\n    pct_list.append(pct)\n    \n    try:\n        pars, pcov = curve_fit(lin_decay, xdata=week, ydata=FVC, p0=[-0.5, 3000])\n        pars_lin.append(pars)\n        pcov_lin.append(pcov)\n        RSS_lin += np.sum((FVC-lin_decay(week,*pars))**2)\n        \n        pars, pcov = curve_fit(exp_decay, xdata=week, ydata=FVC, p0=[-0.5, 500, 2000])\n        pars_exp.append(pars)\n        pcov_exp.append(pcov)\n        RSS_exp += np.sum((FVC-exp_decay(week,*pars))**2)\n        \n        pars, pcov = curve_fit(log_decay, xdata=week, ydata=FVC, p0=[0.5, -0.5, 5000])\n        pars_log.append(pars)\n        pcov_log.append(pcov)\n        RSS_log += np.sum((FVC-log_decay(week,*pars))**2)\n\n        #fit_values.append([pid] + pars.tolist())\n    except RuntimeError:\n        pass\n    \nprint(\"RSS for linear decay model: %.3f \\n\" % RSS_lin)   \nprint(\"RSS for exponential decay model: %.3f \\n\" % RSS_exp) \nprint(\"RSS for logistic decay model: %.3f \\n\" % RSS_log)    \n#     plt.plot(week,pct)\n#     plt.xlabel('Weeks')\n#     plt.ylabel('Percent')\n#fit_raw = pd.DataFrame(fit_values, columns=columns).sort_values(by='r')\n#fit_raw = fit_raw.reset_index(drop=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax = fit_raw['r'].hist()\nax.set_xlabel('r')\nax.set_ylabel('counts')\nax.set_title(\"Initial fitting with individual patient's meta data\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ss_res = np.dot((pct - pct_model(week, *pars)),(pct - pct_model(week, *pars)))\nss_res","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"z=train_df.groupby(['SmokingStatus','Sex'])['FVC'].count().to_frame().reset_index()\nz.columns = ['SmokingStatus', 'Sex', 'Count']\nz.style.background_gradient(cmap='YlOrRd') ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nimdir = '/kaggle/input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/'\n\n\nfig=plt.figure(figsize=(13, 10))\ncolumns = 5\nrows = 3\ninterval = 2\nimlist = os.listdir(imdir)\nfor i in range(1,columns*rows+1):\n    loc = 1\n    filename = imdir + \"/\" + str(i*interval) + \".dcm\"\n    ds = pydicom.dcmread(filename)\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(ds.pixel_array, cmap='YlOrRd')\n    plt.yticks([])\n    plt.xticks([])\n    \n\n# plt.tight_layout()\n# plt.show()","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}