{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import PIL\nimport os, re, time, tqdm\nimport numpy as np\nimport pandas as pd\nfrom tqdm.notebook import tqdm\n\nfrom scipy import stats\nimport shap\n\nfrom sklearn import utils, model_selection\nimport lightgbm as lgb\n\nimport matplotlib.pyplot as plt\nimport cv2\nimport pydicom\nimport tensorflow as tf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_dim = 128\nmin_layers = 200\nquantiles_thrs = [0.1, 0.3, 0.5, 0.7, 0.9]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_tomography(folder, dim=128, norm=False):\n    files = np.array(os.listdir(folder))\n    files = files[np.argsort([int(f.split('.')[0]) for f in files])]\n    images = []\n    for file in files:\n        dicom = pydicom.dcmread(os.path.join(folder, file))\n        try:\n            images += [cv2.resize(dicom.pixel_array, (dim, dim)).reshape(dim, dim)]\n        except RuntimeError:\n            pass\n    images = np.array(images)\n    if np.ndim(images) <= 2:\n        images = np.nan * np.ones((1, dim, dim))\n    if norm:\n        images = (images - images.min()) / (images.max() - images.min())\n    return np.array(images)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"df_train = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\ndf_test = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/test.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def movement(x, axis=None):\n    return np.mean(np.abs(np.diff(x, axis=axis)), axis=axis)\n\n\ndef compl_axis(ax):\n    assert isinstance(ax, list) or isinstance(ax, int)\n    axs = list()\n    if isinstance(ax, list):\n        for a in ax:\n            if isinstance(a, int):\n                axs.append(a)\n            elif isinstance(a, tuple) or isinstance(a, list):\n                axs.extend(list(a))\n    else:\n        axs = [ax]\n    return tuple([i for i in range(3) if i not in axs])\n\n\ndef build_fn2(fns, ax):\n    if isinstance(ax, int) or isinstance(ax, tuple):\n        ax = [ax]\n    assert len(fns) == 2\n    assert len(ax) == 1\n    def fn(x):\n        x = np.expand_dims(fns[0](x, axis=ax[0]), ax[0])\n        c_ax = compl_axis(ax)\n        # x = fns[1](x, axis=c_ax)\n        x = fns[1](np.ravel(x), axis=0)\n        x = float(np.squeeze(x))\n        return x\n    return fn\n\n\ndef build_fn3(fns, axs):\n    assert len(fns) == 3\n    assert len(axs) == 2\n    def fn(x):\n        x = np.expand_dims(fns[0](x, axis=axs[0]), axs[0])\n        x = np.expand_dims(fns[1](x, axis=axs[1]), axs[1])\n        x = fns[2](np.ravel(x), axis=0)\n        x = float(np.squeeze(x))\n        return x\n    return fn","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features_funcs = {\n    'mean': np.mean,\n    'std': np.std,\n}\nbase_aggrs = [('mean', np.mean), ('std', np.std)]\naggrs = base_aggrs + [('kurtosis', stats.kurtosis), ('skew', stats.skew)]\naggrs_mono = aggrs + [('move', movement)]\naxis = [('top_down', 0), ('left_right', 1), ('front_rear', 2)]\nplanes = [('frontal', (0, 1)), ('sagittal', (0, 2)), ('trasversal', (1, 2))]\n# plane-axis\nfor pl_name, pl in planes:\n    for fn0_name, fn0 in base_aggrs:\n        for fn1_name, fn1 in aggrs_mono:\n            if not (fn0_name == 'mean' and fn1_name == 'mean'):\n                features_funcs['%s__%s_%s' % (fn1_name, pl_name, fn0_name)] = build_fn2([fn0, fn1], [pl])\n# axis-plane\nfor ax_name, ax in axis:\n    for fn0_name, fn0 in base_aggrs:\n        for fn1_name, fn1 in aggrs:\n            if not (fn0_name == 'mean' and fn1_name == 'mean'):\n                features_funcs['%s__%s_%s' % (fn1_name, pl_name, fn0_name)] = build_fn2([fn0, fn1], [ax])\n# axis-axis-axis\nfor ax0_name, ax0 in axis:\n    for ax1_name, ax1 in axis:\n        for fn0_name, fn0 in base_aggrs:\n            for fn1_name, fn1 in aggrs_mono:\n                for fn2_name, fn2 in aggrs_mono:\n                    if ax0 != ax1 and (not (fn0_name == 'mean' and fn1_name == 'mean')) and (not (fn1_name == 'mean' and fn2_name == 'mean')):\n                        features_funcs['%s__%s_%s__%s_%s' % (fn2_name, ax1_name, fn1_name, ax0_name, fn0_name)] = build_fn3([fn0, fn1, fn2], [ax0, ax1])\nprint(\"%d features generated\" % len(features_funcs))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features = pd.DataFrame()\nfor i, patient in enumerate(tqdm(df_train.Patient.drop_duplicates(), desc='Patient loop')):\n    folder = os.path.join('../input/osic-pulmonary-fibrosis-progression/train', patient)\n    images = read_tomography(folder, dim=img_dim, norm=True)\n    depth = images.shape[0]\n    while images.shape[0] > min_layers:\n        images = images[::2]\n    f = {'Patient': patient, 'imgs_num': len(os.listdir(folder))}\n    for key in features_funcs.keys():\n        f['img__' + key] = features_funcs[key](images)\n        for thr in quantiles_thrs:\n            f['img__' + key + \"__thr%2d\" % int(100 * thr)] = features_funcs[key](images >= np.quantile(images, thr))\n    features = features.append(pd.DataFrame(f, index=[i]))\n    pass","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(6, 12))\nplt.subplot(2, 1, 1)\nplt.imshow(images[0])\nplt.subplot(2, 1, 2)\nm = images[0].mean(0); plt.plot((m - m.mean()) / m.std(), color='tab:blue')\nm = images[0].std(0); plt.plot((m - m.mean()) / m.std(), color='tab:orange')\nm = movement(images[0], 0); plt.plot((m - m.mean()) / m.std(), color='tab:green')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features.to_csv(\"all_features.csv\", index=False)\nfeatures.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features_test = pd.DataFrame()\nfor i, patient in enumerate(tqdm(df_test.Patient.drop_duplicates(), desc='Patient loop')):\n    folder = os.path.join('../input/osic-pulmonary-fibrosis-progression/test', patient)\n    images = read_tomography(folder, dim=img_dim, norm=True)\n    depth = images.shape[0]\n    while images.shape[0] > min_layers:\n        images = images[::2]\n    f = {'Patient': patient, 'imgs_num': len(os.listdir(folder))}\n    for key in features_funcs.keys():\n        f['img__' + key] = features_funcs[key](images)\n        for thr in quantiles_thrs:\n            f['img__' + key + \"__thr%2d\" % int(100 * thr)] = features_funcs[key](images >= np.quantile(images, thr))\n    features_test = features_test.append(pd.DataFrame(f, index=[i]))\n    pass","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features_test.to_csv(\"all_features_test.csv\", index=False)\nfeatures_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Warning 💀⚡⚠️ Overfitting alert\n\nThe risk of overfitting is significant, use at your own risk\n\nStrong regularization and variable selection are recommended (example below)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"features = pd.read_csv(\"all_features.csv\")\nprint(\"Features has %d columns\" % features.shape[1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Removing duplicates","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"MAX_CORRELATION = 0.995\n\ncorr = features.iloc[:100].replace(np.nan, features.iloc[:1000].mean(0)).corr()\nto_check = []\nfor f0 in corr.index:\n    f1s = corr.loc[f0].index[corr.loc[f0] > MAX_CORRELATION]\n    for f1 in f1s:\n        if f0 != f1:\n            to_check += [f0, f1]\nto_check = list(set(to_check))\ncorr = features[to_check].replace(np.nan, features[to_check].mean(0)).corr()\nduplicated = []\nfor i0, f0 in enumerate(tqdm(corr.index)):\n    f1s = corr.loc[f0].index[(corr.loc[f0] > MAX_CORRELATION) & (np.array(range(corr.shape[1])) > i)] \n    for f1 in f1s:\n        if f0 != f1:\n            duplicated += [(f0, f1)]\nprint(\"%d (almost) duplicated features found\" % len(duplicated))\nfeatures = features.drop([d for _, d in duplicated], axis=1)\nprint(\"Features has %d columns\" % features.shape[1])\nfeatures.to_csv(\"features_no_duplicates.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Iterated permutation importance for variable selection","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"features = pd.read_csv(\"features_no_duplicates.csv\")\ntrain_base = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv').sort_values(['Patient', 'Weeks'])\ntrain = train_base.copy()\nref = train.copy()[['Weeks', 'Patient', 'FVC', 'Percent', 'Age', 'Sex', 'SmokingStatus']]\ntar = train.copy()[['Weeks', 'Patient', 'FVC']]\ntrain = pd.merge(ref, tar, on='Patient', suffixes=['_base', ''])\ntrain['Week_passed'] = train['Weeks'] - train['Weeks_base']\ntrain = train[train['Week_passed'] != 0]\nbase_FE = ['FVC_base', 'Week_passed', 'Age', 'Percent']\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"metric = 'l1'\nselected_features = list(base_FE)\nres_features = list(features.drop('Patient', 1).columns)\nparams = {'n_estimators': 250, 'learning_rate': 0.2, 'reg_alpha': 1.0, 'reg_lambda': 1.0, 'subsample': 0.9}\nX, y = train.merge(features, on='Patient', how='left'), train.FVC\nidxT, idxV = list(model_selection.GroupKFold(10).split(X, y, train.Patient))[0]\nX_train, y_train = X.iloc[idxT], y.iloc[idxT]\nX_valid, y_valid = X.iloc[idxV], y.iloc[idxV]\ninteresting_ones = []\nto_remove = []\nfor _ in range(100):\n    model = lgb.LGBMRegressor(**params)\n    model.fit(X_train[selected_features], y_train, eval_set=[(X_valid[selected_features], y_valid)], eval_metric=metric, verbose=0)\n    baseline = min(model.evals_result_['valid_0'][metric])\n    print(\"Baseline\\t\\t%.4f\" % baseline)\n    best_val = baseline\n    best_f = None\n    trange = tqdm(res_features)\n    for f in trange:\n        if f in selected_features or f in to_remove:\n            continue\n        model.fit(X_train[selected_features + [f]], y_train, eval_set=[(X_valid[selected_features + [f]], y_valid)], eval_metric=metric, verbose=0)\n        val = min(model.evals_result_['valid_0'][metric])\n        if val < baseline:\n            interesting_ones.append(f)\n        if val < best_val:\n            best_val = val\n            best_f = f\n            trange.set_description(\"Found %.4f\" % best_val)\n            \n    if best_f is not None:\n        print(\"+%s\\t\\t%.4f\" % (best_f, best_val))\n        selected_features.append(best_f)\n    else:\n        print(\"STOP\")\n        break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features[['Patient'] + [f for f in selected_features if f in features.columns]].to_csv(\"features_selected.csv\", index=False)\nfeatures[['Patient'] + [f for f in list(set(interesting_ones)) if f in features.columns]].to_csv(\"features_interesting.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'learning_rate': 0.02, 'reg_lambda': 50, 'reg_alpha': 50, 'n_estimators': 1000, 'subsample': 0.5}\nmodel = lgb.LGBMRegressor(**params)\nmodel.fit(X_train[base_FE], y_train, eval_set=[(X_valid[base_FE], y_valid)], eval_metric=metric, verbose=0)\nmodel_imgs = lgb.LGBMRegressor(**params)\nmodel_imgs.fit(X_train[selected_features], y_train, eval_set=[(X_valid[selected_features], y_valid)], eval_metric=metric, verbose=0)\nplt.plot(model.evals_result_['valid_0'][metric], color='tab:blue', label='base data')\nplt.plot(model_imgs.evals_result_['valid_0'][metric], color='tab:red', linestyle='--', label='base data + imgs')\nplt.yscale('log'); plt.legend();\nprint(\"%.6f ==> %.6f\" % (min(model.evals_result_['valid_0'][metric]), min(model_imgs.evals_result_['valid_0'][metric])))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}