{"cells":[{"metadata":{},"cell_type":"markdown","source":"- First of all, congratulations to all kagglers and thanks to organizers. I have learned a lot from this competition, so I would like to share my solution here. I hope some pepole find new ideas from my solution.\n- This notebook is based on the following great notebooks.\n  - https://www.kaggle.com/khoongweihao/efficientnets-quantile-regression-inference\n  - https://www.kaggle.com/ttahara/osic-baseline-lgbm-with-custom-metric\n  - https://www.kaggle.com/allunia/pulmonary-dicom-preprocessing\n  - https://www.kaggle.com/unforgiven/osic-comprehensive-eda"},{"metadata":{},"cell_type":"markdown","source":"# Imports"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index\n!pip install ../input/efficientnet/efficientnet-1.1.0/ -f ./ --no-index\n#!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import os\nimport gc\nimport cv2\nimport random\nimport pydicom\nimport operator\nimport typing as tp\nimport pandas as pd\nimport numpy as np \nimport pymc3 as pm\nimport tensorflow as tf \nimport matplotlib.pyplot as plt \nimport seaborn as sns\nfrom PIL import Image\nfrom tqdm.notebook import tqdm\nfrom datetime import datetime\nfrom IPython.display import display\n\nfrom scipy.stats import describe\nfrom scipy.ndimage import binary_fill_holes\nfrom skimage import measure, morphology\nfrom skimage.filters import threshold_otsu, median\nfrom skimage.segmentation import clear_border\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler\nfrom sklearn.model_selection import train_test_split, StratifiedKFold, KFold, TimeSeriesSplit, GroupKFold, GroupShuffleSplit\nfrom optuna.visualization import plot_optimization_history\n\nimport optuna\nimport lightgbm as lgb\nimport efficientnet.tfkeras as efn\nimport tensorflow_addons as tfa\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.activations as activations\nfrom tensorflow.keras import Model\nfrom tensorflow.keras.utils import Sequence\nfrom tensorflow.keras.optimizers import Nadam\nfrom tensorflow_addons.optimizers import RectifiedAdam\nfrom tensorflow.keras.layers import (\n    Dense, Dropout, Activation, Flatten, Input, BatchNormalization, GlobalAveragePooling2D, Add, Conv2D, AveragePooling2D, \n    LeakyReLU, Concatenate, DepthwiseConv2D, Maximum, Layer)\n\ndef 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\npd.set_option(\"max_columns\", 100)\npd.set_option('display.max_rows', 150)\nseed_everything(42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"config  = tf.compat.v1.ConfigProto()\nconfig.gpu_options.allow_growth = True\nsession = tf.compat.v1.Session(config=config)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BAD_IDS = ['ID00011637202177653955184', 'ID00052637202186188008618']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Linear Decay (based on EfficientNets)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define the last 3 and last 6 in order to get solid cv.\ndef add_last(data, add3=True, add_other=0):\n    data[\"last3\"] = False\n    if add_other != 0:\n        data[\"last\"+str(add_other)] = False\n    for patient in data.Patient.unique():\n        weeks   = data.query(\"Patient == @patient & WHERE == 'train'\").Weeks\n        weeks_3 = sorted(weeks)[-3:]\n        weeks_other = sorted(weeks)[-add_other:]\n        data.loc[(data.Patient==patient)&(data.Weeks.isin(weeks_3)), \"last3\"] = True\n        if add_other != 0:\n            data.loc[(data.Patient==patient)&(data.Weeks.isin(weeks_other)), \"last\"+str(add_other)] = True\n    return data\n\ndef get_tab(df):\n    vector = [(df.Age.values[0] - 30) / 30] \n    if df.Sex.values[0] == 'male':\n        vector.append(0)\n    else:\n        vector.append(1)\n    if   df.SmokingStatus.values[0] == 'Never smoked':\n        vector.extend([0,0])\n    elif df.SmokingStatus.values[0] == 'Ex-smoker':\n        vector.extend([1,1])\n    elif df.SmokingStatus.values[0] == 'Currently smokes':\n        vector.extend([0,1])\n    else:\n        vector.extend([1,0])\n    return np.array(vector) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Funnel ReLu(https://arxiv.org/abs/2007.11824)\nclass FReLU(Layer):\n    def __init__(self, kernel_size=3):\n        super(FReLU, self).__init__()\n        self.conv = DepthwiseConv2D(kernel_size=kernel_size, padding=\"same\")\n        self.bn   = BatchNormalization()\n\n    def call(self, inputs):\n        x1 = self.conv(inputs)\n        x1 = self.bn(x1)\n        x  = K.maximum(inputs, x1)\n        return x\n\n    def compute_output_shape(self, input_shape):\n        return input_shape\n    \n    def get_config(self):\n        config = super().get_config().copy()\n        config.update({\n            'conv': self.conv,\n            'bn': self.bn})\n        return config\n\ndef get_efficientnet(model, shape):\n    models_dict = {\n        'b0': efn.EfficientNetB0(input_shape=shape,weights=None,include_top=False),\n        'b1': efn.EfficientNetB1(input_shape=shape,weights=None,include_top=False),\n        'b2': efn.EfficientNetB2(input_shape=shape,weights=None,include_top=False),\n        'b3': efn.EfficientNetB3(input_shape=shape,weights=None,include_top=False),\n        'b4': efn.EfficientNetB4(input_shape=shape,weights=None,include_top=False),\n        'b5': efn.EfficientNetB5(input_shape=shape,weights=None,include_top=False),\n        'b6': efn.EfficientNetB6(input_shape=shape,weights=None,include_top=False),\n        'b7': efn.EfficientNetB7(input_shape=shape,weights=None,include_top=False)\n    }\n    return models_dict[model]\n\ndef redifine_efn_model(model):\n    for i, layer in enumerate(model.layers):\n        if i==0:\n            input = layer.input\n            x     = input\n            seblock_on = False\n            block_prev = \"\"\n            continue\n            \n        block_curt = layer.name[:layer.name.find(\"_\")]\n        if block_prev != block_curt:\n            input_block = x\n            \n        # Change from current activation to FReLU\n        if   -1 < layer.name.find(\"activation\"):\n            x = FReLU()(x)\n        # SEBLOCKS\n        elif -1 < layer.name.find(\"se\"):\n            if -1 < layer.name.find(\"excite\"):\n                x = layer([x, x_seblock])\n                seblock_on = False\n            elif -1 < layer.name.find(\"reduce\"):\n                layer.activation = activations.linear\n                x_seblock = layer(x_seblock)\n                x_seblock = FReLU()(x_seblock)\n            elif not seblock_on:\n                seblock_on = True \n                x_seblock  = layer(x)\n            else:\n                x_seblock  = layer(x_seblock)\n        # Add Layer\n        elif -1 < layer.name.find(\"add\"):\n            x = layer([x, input_block])\n        # Normal Layer\n        else:\n            x = layer(x)      \n        block_prev = block_curt\n    return Model(input, x)\n\ndef build_model(shape=(512, 512, 1), model_class=None, act=None):\n    inp  = Input(shape=shape)\n    base = get_efficientnet(model_class, shape)\n    if act == \"FReLU\":\n        base = redifine_efn_model(base)\n    x = base(inp)\n    x = GlobalAveragePooling2D()(x)\n    inp2 = Input(shape=(4,))\n    x2 = tf.keras.layers.GaussianNoise(0.2)(inp2)\n    x  = Concatenate()([x, x2]) \n    x  = Dropout(0.5)(x) \n    x  = Dense(1)(x)\n    model = Model([inp, inp2] , x)\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score(fvc_true, fvc_pred, sigma):\n    sigma_clip = np.maximum(sigma, 70) # changed from 70, trie 66.7 too\n    delta  = np.abs(fvc_true - fvc_pred)\n    delta  = np.minimum(delta, 1000)\n    sq2    = np.sqrt(2)\n    metric = (delta / sigma_clip)*sq2 + np.log(sigma_clip* sq2)\n    return np.mean(metric)\n\ndef get_img(path, shapes):\n    def rescale(img, ri, rs):\n        return (img - ri) / (rs * 1000)\n    d    = pydicom.dcmread(path)\n    img  = d.pixel_array\n    img[img <= -1000] = 0\n    img  = img.astype(np.int32)\n    W, H = img.shape\n    if W == H:\n        img = rescale(img, d.RescaleIntercept, d.RescaleSlope)\n        if W == shapes[0]:\n            return img\n        else:\n            return cv2.resize(img, shapes)\n    iimg   = Image.fromarray(img, mode=\"I\")\n    left   = (img.shape[1]-512)/2\n    right  = (img.shape[1]+512)/2\n    top    = (img.shape[0]-512)/2\n    bottom = (img.shape[0]+512)/2\n    iimg   = iimg.crop((left, top, right, bottom))\n    iimg   = iimg.resize(shapes, resample=Image.LANCZOS)\n    return rescale(np.array(iimg), d.RescaleIntercept, d.RescaleSlope)\n\ndef get_result(models, shapes, vl_p):\n    m = [[],[],[],[],[],[],[],[],[]]\n    df_results = pd.DataFrame()\n    metric = []\n    for p in tqdm(vl_p):\n        if p in BAD_IDS: continue\n        x    = [] \n        tab  = [] \n        path = '../input/osic-pulmonary-fibrosis-progression/train/%s/' % p\n\n        ldir = os.listdir(path)\n        ldir = random.sample(ldir, 30) if 30 < len(ldir) else ldir\n        for i in ldir:\n            x.append(get_img(path + i, shapes))\n            tab.append(get_tab(train.loc[train.Patient == p, :])) \n        tab = np.array(tab) \n\n        percent_true = train.Percent.values[train.Patient == p]\n        fvc_true     = train.FVC    .values[train.Patient == p]\n        weeks_true   = train.Weeks  .values[train.Patient == p]\n        \n        x  = np.expand_dims(x, axis=-1)\n        for idx, model in enumerate(models):\n            _a = model.predict([x, tab])\n            for q in range(1, 10):\n                a          = np.quantile(_a, q / 10)\n                fvc        = a * (weeks_true - weeks_true[0]) + fvc_true[0]\n                percent    = percent_true[0] - a * abs(weeks_true - weeks_true[0])\n                lll_score  = score(fvc_true, fvc, percent)\n                df_results = df_results.append(pd.DataFrame({\"model_id\":idx, \"q\":q, \"Patient\":p, \"Weeks\":weeks_true, \n                                                             \"lll_score\":lll_score, \"fvc_true\":fvc_true, \"fvc\":fvc, \"confidence\":percent}))\n    metric = df_results.groupby(\"q\", as_index=False).lll_score.mean()\n    metric = metric.sort_values(by=\"q\").reset_index(drop=True)\n    print(metric)\n    if len(np.unique(metric.lll_score))==1:\n        q = 5\n    else:\n        q = (np.argmin(metric.lll_score) + 1)\n    df_results = df_results.query(\"q == @q\").reset_index(drop=True)\n    q /= 10\n    \n    sub  = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv') \n    test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv') \n    A_test, B_test, P_test, W, FVC = {},{},{},{},{} \n    STD, WEEK = {},{} \n    for p in test.Patient.unique():\n        x = []; tab = []; a_test = 0; b_test = 0\n        path = '../input/osic-pulmonary-fibrosis-progression/test/%s/' % p\n        \n        ldir = os.listdir(path)\n        ldir = random.sample(ldir, 30) if 30 < len(ldir) else ldir\n        for i in ldir:\n            x.append(get_img(path + i, shapes)) \n            tab.append(get_tab(test.loc[test.Patient == p, :])) \n            \n        tab = np.array(tab) \n        x   = np.expand_dims(x, axis=-1)\n        for model in models:\n            _a  = model.predict([x, tab]) \n            a   = np.quantile(_a, q)\n            a_test    += a \n            b_test    += test.FVC    .values[test.Patient == p] - a*test.Weeks.values[test.Patient == p]\n            P_test[p]  = test.Percent.values[test.Patient == p] \n            WEEK[p]    = test.Weeks  .values[test.Patient == p]\n        A_test[p] = a_test / len(models)\n        B_test[p] = b_test / len(models)\n\n    for k in sub.Patient_Week.values:\n        p, w = k.split('_')\n        w    = int(w) \n        fvc  = A_test[p] * w + B_test[p]\n        sub.loc[sub.Patient_Week == k, 'FVC']        = fvc\n        sub.loc[sub.Patient_Week == k, 'Confidence'] = (P_test[p] - A_test[p] * abs(WEEK[p] - w))\n        \n    return df_results, sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv') \ntrain = train.sort_values(by=[\"Patient\",\"Weeks\"]).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"USE_ALL = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nmodels = [[], ['b5']]\npaths  = ['../input/osic-model-weights/', '../input/osic-own-model-weights/']\nacts   = [\"\", \"FReLU\"]\nshape  = (512, 512)\n\nif USE_ALL:\n    vl_p    = train.Patient.unique()\nelse:\n    _, vl_p = train_test_split(train.Patient.unique(), shuffle=True, train_size=0.8)\n    \nloaded_models = []\nprint(\"Loading models\")\nfor ms, path, act in zip(models, paths, acts):\n    for m in ms:\n        if m == \"\": continue\n        files = sorted([f for f in os.listdir(path) if m in f])\n        for f in files:\n            print(path + f)\n            model = build_model((*shape,1), m, act)\n            model.load_weights(path + f)\n            loaded_models.append(model)\n            del model\n            gc.collect()\nprint(\"Calculating results\")\nres, sub  = get_result(loaded_models, shape, vl_p)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_result_img = res.groupby([\"Patient\",\"Weeks\"], as_index=False)[[\"fvc\",\"confidence\"]].mean()\ndf_result_img[\"WHERE\"] = \"train\"\ndf_result_img = add_last(df_result_img)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient   = df_result_img.query(\"last3\").Patient\nlast3_fvc = df_result_img.query(\"last3\").fvc\nlast3_con = df_result_img.query(\"last3\").confidence\nlast3_pred_img = np.array([last3_fvc-(last3_con/2),\n                           last3_fvc,\n                           last3_fvc+(last3_con/2)]).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(sub.shape)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission_img.csv\", index=False)\n\nimg_sub = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# DICOM preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_slice_thickness(s0, s1):\n    if \"ImagePositionPatient\" in s0.dir():\n        return np.abs(s0.ImagePositionPatient[2] - s1.ImagePositionPatient[2])\n    if \"SliceLocation\" in s0.dir():\n        return np.abs(s0.SliceLocation - s1.SliceLocation)\n    return s0.SliceThickness\n\ndef load_scans(dcm_path):\n    files     = os.listdir(dcm_path)\n    file_nums = [np.int(file.split(\".\")[0]) for file in files]\n    sorted_file_nums = np.sort(file_nums)[::-1]\n    slices = [pydicom.dcmread(dcm_path + \"/\" + str(file_num) + \".dcm\" ) for file_num in sorted_file_nums]\n    st     = get_slice_thickness(slices[0],slices[1])\n    return slices, st, slices[0].PixelSpacing\n\ndef transform_to_hu(slices):\n    images = np.stack([file.pixel_array for file in slices]).astype(np.int16)\n\n    # convert ouside pixel-values to air:\n    # I'm using <= -1000 to be sure that other defaults are captured as well\n    images[images <= -1000] = 0\n    \n    # convert to HU\n    intercept = slices[0].RescaleIntercept\n    slope     = slices[0].RescaleSlope\n    for n in range(len(images)):\n        \n        if slope != 1:\n            images[n] = slope * images[n].astype(np.float64)\n            images[n] = images[n].astype(np.int16)\n        images[n] += np.int16(intercept)\n    \n    return np.array(images, dtype=np.int16)\n\ndef lung_segment(img):\n    if len(np.unique(img))==1:\n        lungs  = np.zeros_like(img)\n    else:\n        thresh = threshold_otsu(img)\n        binary = img <= thresh\n        lungs  = median(clear_border(binary))\n        lungs  = morphology.binary_closing(lungs, selem=morphology.disk(7))\n        lungs  = binary_fill_holes(lungs)\n\n    final = lungs*img\n    final[final == 0] = np.min(img)\n\n    return final, lungs\n\ndef lung_segment_stack(imgs):\n    masks    = np.empty_like(imgs)\n    segments = np.empty_like(imgs)\n\n    for i, img in enumerate(imgs):\n        seg, mask = lung_segment(img)\n        segments[i,:,:] = seg\n        masks[i,:,:]    = mask\n        \n    return segments, masks\n\ndef lung_volume(st, ps, masks):\n    return np.round(np.sum(masks) * st * ps[0] * ps[1], 3)\n\ndef hist_analysis(segmented):\n    values = segmented.flatten()\n    values = values[values >= -1000]\n    stats  = describe(values)\n    return stats.mean, stats.variance, stats.skewness, stats.kurtosis\n\ndef chest_measurements(ps, masks):\n    middle_slice = masks[len(masks)//2]\n    lung_area    = np.round(np.sum(middle_slice.flatten()) * ps[0] * ps[1], 3)\n    conv_h = morphology.convex_hull_image(middle_slice)\n    props  = measure.regionprops(measure.label(conv_h))\n    chest_diameter = np.round(props[0].major_axis_length, 3)\n    chest_circ     = np.round(props[0].perimeter, 3)\n    return lung_area, chest_diameter, chest_circ\n\ndef resize_scan(scan, new_shape):\n    # read slice as 32 bit signed integers\n    img = Image.fromarray(scan, mode=\"I\")\n    # do the resizing\n    img = img.resize(new_shape, resample=Image.LANCZOS)\n    # convert back to 16 bit integers\n    resized_scan = np.array(img, dtype=np.int16)\n    return resized_scan\n\ndef crop_scan(scan):\n    img = Image.fromarray(scan, mode=\"I\")\n    \n    left   = (scan.shape[0]-512)/2\n    right  = (scan.shape[0]+512)/2\n    top    = (scan.shape[1]-512)/2\n    bottom = (scan.shape[1]+512)/2\n\n    img = img.crop((left, top, right, bottom))\n    # convert back to 16 bit integers\n    cropped_scan = np.array(img, dtype=np.int16)\n    return cropped_scan\n\ndef crop_and_resize(scan, new_shape):\n    img = Image.fromarray(scan, mode=\"I\")\n    \n    left   = (scan.shape[0]-512)/2\n    right  = (scan.shape[0]+512)/2\n    top    = (scan.shape[1]-512)/2\n    bottom = (scan.shape[1]+512)/2\n    \n    img = img.crop((left, top, right, bottom))\n    img = img.resize(new_shape, resample=Image.LANCZOS)\n    \n    cropped_resized_scan = np.array(img, dtype=np.int16)\n    return cropped_resized_scan\n\ndef resize_image(images, resized_shape):\n    images         = images.astype(np.int32)\n    resized_images = np.zeros((images.shape[0], resized_shape[0], resized_shape[1]),\n                              dtype=np.int16)\n    # if squared:\n    if images.shape[1] == images.shape[2]:\n        # if size is not as desired\n        if images.shape[1] == resized_shape[0]:\n            return images\n        # else resize:\n        else:\n            for s in range(images.shape[0]):\n                resized_images[s] = resize_scan(images[s,:,:], resized_shape)\n\n    # if non-squared - do a center crop to 512, 512 and then resize to desired shape\n    else:\n        for s in range(images.shape[0]):\n            # if desired shape is 512x512:\n            if resized_shape[0]==512:\n                resized_images[s] = crop_scan(images[s,:,:])\n            else:\n                resized_images[s] = crop_and_resize(images[s,:,:], resized_shape)\n    return resized_images\n\ndef calculate_lung_area(path, train_test, resize_shape=(512,512)):\n    dicom_data = []\n    for i, patient in enumerate(tqdm(os.listdir(path + train_test))):\n        if patient in ['ID00011637202177653955184', 'ID00052637202186188008618']:\n            continue\n        scans, st, ps = load_scans(path + train_test + \"/\" + patient + \"/\")\n        images = transform_to_hu(scans)\n        images = resize_image(images, resize_shape)\n        segmented, masks = lung_segment_stack(images)\n        del scans, images\n        gc.collect()\n        vol     = lung_volume(st, ps, masks)\n        m,v,s,k = hist_analysis(segmented)\n        lung_area, chest_diameter, chest_circ = chest_measurements(ps, masks)\n        del segmented, masks\n        gc.collect()\n        dicom_data.append([patient, train_test, vol, m, v, s, k, \n                           lung_area, chest_diameter, chest_circ])\n        \n    cols = [\"Patient\",\"WHERE\",\"volume\",\"mean\",\"variance\",\"skewness\",\"kurtosis\",\n            \"lung_area\",\"chest_diameter\",\"chest_circ\"]\n    return pd.DataFrame(dicom_data, columns=cols)\n\ndef get_dicom_Data(method=\"type1\", path=\"\", resize=None):\n    if method==\"type1\":\n        def get_window_value(feature):\n            if type(feature) == pydicom.multival.MultiValue: return np.int(feature[0])\n            return np.int(feature)\n\n        def load_dicom_data(train_test):\n            dicom_data = []\n            path     = ROOT + train_test + \"/\"\n            patients = os.listdir(path)\n            for patient in patients:\n                example_dcm = os.listdir(path + patient)[0]\n                dataset     = pydicom.dcmread(path + patient + \"/\" + example_dcm)\n                spacing     = dataset.PixelSpacing\n                dicom_data.append([patient, train_test, dataset.Rows, dataset.Columns,\n                                   get_window_value(dataset.WindowWidth),\n                                   get_window_value(dataset.WindowCenter),\n                                   dataset.SliceThickness, spacing[0], spacing[1]])\n            return dicom_data\n            \n        dicom_data_train = load_dicom_data(\"train\")\n        dicom_data_test  = load_dicom_data(\"test\")\n\n        cols = [\"Patient\",\"WHERE\",\"rows\",\"columns\",\"window_width\",\"window_level\",\n                \"slice_thickness\",\"ps_r\",\"ps_c\"]\n        df_dicom_train = pd.DataFrame(dicom_data_train, columns=cols)\n        df_dicom_test  = pd.DataFrame(dicom_data_test,  columns=cols)\n        df_dicom       = pd.concat([df_dicom_train, df_dicom_test]).reset_index(drop=True)\n        df_dicom[\"area\"]       = df_dicom.rows * df_dicom[\"columns\"]\n        df_dicom[\"ps_area\"]    = df_dicom.ps_r * df_dicom.ps_c\n        df_dicom[\"r_distance\"] = df_dicom.ps_r * df_dicom.rows\n        df_dicom[\"c_distance\"] = df_dicom.ps_c * df_dicom[\"columns\"]\n        df_dicom[\"area_cm2\"]         = 0.1 * df_dicom.r_distance * 0.1 * df_dicom.c_distance\n        df_dicom[\"slice_volume_cm3\"] = 0.1 * df_dicom.slice_thickness * df_dicom.area_cm2\n        df_dicom = df_dicom[[\"Patient\",\"WHERE\",\"ps_area\",\"area_cm2\",\"slice_volume_cm3\"]]\n        \n    else:\n        df_dicom_train = calculate_lung_area(path, \"train\", resize)\n        df_dicom_test  = calculate_lung_area(path, \"test\",  resize)\n        df_dicom       = pd.concat([df_dicom_train, df_dicom_test]).reset_index(drop=True)\n    return df_dicom","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ADD_DICOM = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nmethod   = \"type2\"\nbasepath = \"../input/osic-pulmonary-fibrosis-progression/\"\nresize   = (256,256)\nif ADD_DICOM:\n    dicom_data = get_dicom_Data(method, basepath, resize)\n    dicom_sub  = dicom_data.query(\"WHERE=='test'\").copy()\n    dicom_sub[\"WHERE\"] = \"sub\"\n    dicom_data = pd.concat([dicom_data, dicom_sub]).reset_index(drop=True)\n    display(dicom_data.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Osic-Multiple-Quantile-Regression"},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression/\"\nBATCH_SIZE = 128","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def preparation(path):\n    # train data\n    train = pd.read_csv(f\"{path}/train.csv\")\n    train.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\n    if USE_ALL:\n        train = train.query(\"Patient not in @BAD_IDS\").reset_index(drop=True)\n    # test data\n    chunk = pd.read_csv(f\"{path}/test.csv\")\n    # submission data\n    sub = pd.read_csv(f\"{path}/sample_submission.csv\")\n    sub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\n    sub['Weeks']   = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\n    sub = sub[['Patient','Weeks','Confidence','Patient_Week']]\n    sub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")\n\n    train['WHERE'] = 'train'\n    chunk['WHERE'] = 'val'\n    sub  ['WHERE'] = 'test'\n    data = train.append([chunk, sub])\n    data = add_last(data, add_other=6)\n    return data\n\ndef make_base(data):\n    data['min_week'] = data['Weeks']\n    data.loc[data.WHERE=='test','min_week'] = np.nan\n    data['min_week'] = data.groupby('Patient')['min_week'].transform('min')\n    \n    base = data.loc[data.Weeks == data.min_week][['Patient','FVC',\"Percent\"]].copy()\n    base.columns = ['Patient','min_FVC','min_Percent']\n    base['nb']   = 1\n    base['nb']   = base.groupby('Patient')['nb'].transform('cumsum')\n    base = base[base.nb==1]\n    base.drop('nb', axis=1, inplace=True)\n    return base","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def feature_engineering(data, target=False):\n    \n    def calculate_height(row):\n        if row['Sex'] == 'Male':\n            return row['min_FVC'] / (27.63 - 0.112 * row['Age'])\n        else:\n            return row['min_FVC'] / (21.78 - 0.101 * row['Age'])\n        \n    def normalize(d, c, c_new, return_vals=False, std=False):\n        d = d.copy()\n        if c not in d.columns: return d\n        if std:\n            c_mean = d[c].mean()\n            c_std  = d[c].std()\n            d[c_new] = (d[c] - d[c].mean()) / d[c].std()\n            if return_vals:\n                return d, c_mean, c_std\n        else:\n            c_min = d[c].min()\n            c_max = d[c].max()\n            d[c_new] = (d[c] - d[c].min()) / (d[c].max() - d[c].min())\n            if return_vals:\n                return d, c_min, c_max\n        return d\n    \n    FE   = []\n    dict = {\"Never smoked\"    : 0,\n            \"Ex-smoker\"       : 0.7,\n            \"Currently smokes\": 1,\n            \"Male\"   : 0,\n            \"Female\" : 1}\n    \n    data = data.copy()\n    # Make features\n    data[\"smoke\"]     = data.SmokingStatus.apply(lambda x: dict[x])\n    data[\"sex\"]       = data.Sex.apply(lambda x: dict[x])\n    data['base_week'] = data['Weeks'] - data['min_week']\n    data['height']    = data.apply(calculate_height, axis=1)\n    # Make fvc features\n    data[\"fvc_expected\"]  = data['min_FVC'] * (100/data['min_Percent'])\n    data[\"fvc_efvc_diff\"] = data['min_FVC'] - data['fvc_expected']\n    data[\"fvc_efvc_div\"]  = data['min_FVC'] / data['fvc_expected']\n    data[\"fvc_per_mul\"]   = data['min_FVC'] * data['min_Percent']\n    data[\"fvc_per_div\"]   = data['min_FVC'] / data['min_Percent']\n    # Normalize\n    data = normalize(data, \"Age\",           \"age\")\n    data = normalize(data, \"base_week\",     \"week\")\n    data = normalize(data, \"height\",        \"height\")\n    data = normalize(data, \"min_FVC\",       \"base_fvc\")\n    data = normalize(data, \"min_Percent\",   \"base_percent\")\n    data = normalize(data, \"fvc_expected\",  \"fvc_expected\")\n    data = normalize(data, \"fvc_efvc_diff\", \"fvc_efvc_diff\")\n    data = normalize(data, \"fvc_efvc_div\",  \"fvc_efvc_div\")\n    data = normalize(data, \"fvc_per_mul\",   \"fvc_per_mul\")\n    data = normalize(data, \"fvc_per_div\",   \"fvc_per_div\")\n    # Normalize for DICOM data\n    data = normalize(data, \"area_cm2\",         \"area_cm2\")\n    data = normalize(data, \"slice_volume_cm3\", \"slice_volume_cm3\")\n    data = normalize(data, \"ps_area\",          \"ps_area\")\n    data = normalize(data, \"volume\",           \"lung_volume\")\n    data = normalize(data, \"mean\",             \"lung_mean\")\n    data = normalize(data, \"skewness\",         \"lung_skew\")\n    data = normalize(data, \"kurtosis\",         \"lung_kurt\")\n    data = normalize(data, \"lung_area\",        \"lung_area\")\n    data = normalize(data, \"chest_diameter\",   \"chest_diameter\")\n    data = normalize(data, \"chest_circ\",       \"chest_circ\")\n    # Make week features\n    data[\"week_age\"]     = data['week'] * data['age']\n    data[\"week_fvc\"]     = data['week'] * data['base_fvc']\n    data[\"week_per\"]     = data['week'] * data['base_percent']\n    data[\"week_fvc_per\"] = data['week'] * data['fvc_per_mul']\n\n    FE += ['age','week',\"height\",'base_fvc',\"base_percent\",\"smoke\",\"sex\",\n           \"fvc_expected\",\"fvc_efvc_diff\",\"fvc_per_mul\",\n           \"week_age\",\"week_fvc\",\"week_per\",\"week_fvc_per\"]\n    if   \"area_cm2\" in data.columns:\n        FE += [\"ps_area\",\"area_cm2\",\"slice_volume_cm3\"]\n    elif \"lung_volume\" in data.columns:\n        FE += [\"lung_volume\"]\n    \n    if target:\n        data, fvc_min, fvc_max = normalize(data, \"FVC\", \"FVC_norm\", return_vals=True, std=True)\n        return data, FE, fvc_min, fvc_max\n    \n    return data, FE","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = preparation(ROOT)\nbase = make_base(data)\n\ndata = data.merge(base, on='Patient', how='left')\nprint(data.shape)\nif ADD_DICOM:\n    data = data.merge(dicom_data, on=['Patient',\"WHERE\"])\n    print(data.shape)\ndata.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data, FE = feature_engineering(data)\n\nsave_cols = [\"Patient\",\"FVC\",\"WHERE\",'Patient_Week','Confidence',\"last3\",\"last6\"] + FE\ntrain = data.loc[data.WHERE=='train', save_cols].sort_values(by=[\"Patient\",\"week\"]).reset_index(drop=True).copy()\nchunk = data.loc[data.WHERE=='val',   save_cols].reset_index(drop=True).copy()\nsub   = data.loc[data.WHERE=='test',  save_cols].reset_index(drop=True).copy()\ndel data, base\ngc.collect()\n\nprint(train.shape, chunk.shape, sub.shape)\ndisplay(train.head(2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(FE)\nprint(train[FE].shape)\ntrain[FE].head()","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 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, C1)\n    delta  = np.minimum(np.abs(actual_fvc - predicted_fvc), C2)\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\ndef calculate_scores(y, pred, return_values=False, show=True):\n    lll = laplace_log_likelihood(y, pred[:, 1], pred[:, 2]-pred[:,0])\n    sigma_opt  = mean_absolute_error(y, pred[:, 1])\n    unc        = pred[:,2] - pred[:, 0]\n    sigma_mean = np.mean(unc)\n    if show:\n        print(\"Laplace log likelihood\", lll)\n        print(\"mean_absolute_error\", sigma_opt)\n        print(\"uncertainty\", sigma_mean, unc.min(), unc.max())\n    \n    if return_values:\n        return lll\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 = 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.25, 0.50, 0.75]\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\ndef make_model(nh):\n    z  = L.Input((nh,), name=\"Patient\")\n    x  = L.Dense(100, activation=\"relu\", kernel_initializer='he_normal', name=\"d1\")(z)\n    x  = L.Dense(100, activation=\"relu\", kernel_initializer='he_normal', name=\"d2\")(x)\n    x  = L.Dense(100, activation=\"relu\", kernel_initializer='he_normal', 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=mloss(0.65),\n                  optimizer=tf.keras.optimizers.Adam(lr=0.1, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False),\n                  metrics=[score])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"er = tf.keras.callbacks.EarlyStopping(\n    monitor=\"val_score\",\n    min_delta=0,\n    patience=500,\n    verbose=1,\n    mode=\"auto\",\n    baseline=None,\n    restore_best_weights=True)\n\nrlp = tf.keras.callbacks.ReduceLROnPlateau(\n    monitor='val_loss', \n    factor=0.6,\n    patience=150, \n    verbose=1, \n    min_lr=1e-8)\n\nlr_schedule = tf.keras.experimental.CosineDecayRestarts(1e-2, 1e1)\nlrs_cdr     = tf.keras.callbacks.LearningRateScheduler(lr_schedule)\n\n# Scheduler didn't work well, so I used only EarlyStopping.\ncallbacks   = [er]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NFOLD  = 5\nEPOCHS = 1200","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ng    = train.Patient\nze   = sub[FE].values\npe   = np.zeros((ze.shape[0], 3))\npred = np.zeros((train.shape[0], 3))\nkf   = GroupKFold(n_splits=NFOLD)\n\nfor cnt, (tr_idx, val_idx) in enumerate(kf.split(train, groups=g)):\n    print(f\"FOLD {cnt+1}\")\n    train_x = (train.iloc[tr_idx, :]).query(\"last6\")[FE].values\n    train_y = (train.iloc[tr_idx, :]).query(\"last6\")[\"FVC\"].values\n    valid_idx = (train.iloc[val_idx,:]).query(\"last3\").index\n    valid_x   = train.loc[valid_idx, FE].values\n    valid_y   = train.loc[valid_idx, \"FVC\"].values\n    net = make_model(len(FE))\n    net.fit(train_x, train_y,\n            batch_size=BATCH_SIZE,\n            epochs=EPOCHS, \n            callbacks=callbacks,\n            validation_data=(valid_x, valid_y),\n            verbose=0)\n    print(\"train\", net.evaluate(train_x, train_y, verbose=0, batch_size=BATCH_SIZE))\n    print(\"valid\", net.evaluate(valid_x, valid_y, verbose=0, batch_size=BATCH_SIZE))\n    pred[valid_idx] = net.predict (valid_x, verbose=0, batch_size=BATCH_SIZE)\n    pe += net.predict(ze, batch_size=BATCH_SIZE, verbose=0) / NFOLD\n    del net\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"last3_idx  = np.where(pred.sum(1)>0)\nlast3_pred = pred[last3_idx]\nlast3_fvc  = train.iloc[last3_idx][\"FVC\"]\n# Reverse normalization\n#last3_pred = last3_pred * (fvc_max-fvc_min) + fvc_min\n#last3_pred = last3_pred * fvc_std + fvc_mean\n\ncalculate_scores(last3_fvc, last3_pred)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Light GBM"},{"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    \nclass BaseModel(object):\n    \"\"\"\n    Base Model Class:\n\n    train_df         : train pandas dataframe\n    test_df          : test pandas dataframe\n    target           : target column name (str)\n    features         : list of feature names\n    categoricals     : list of categorical feature names\n    n_splits         : K in KFold (default is 3)\n    cv_method        : options are .. KFold, StratifiedKFold, TimeSeriesSplit, GroupKFold, or GroupShuffleSplit\n    group            : group feature name when GroupKFold or StratifiedGroupKFold are used\n    task             : options are .. regression, multiclass, or binary\n    param            : dict of parameter, set that if you already define\n    parameter_tuning : bool, only for LGB\n    seed             : seed (int)\n    verbose          : bool\n    \"\"\"\n\n    def __init__(self, train_df, test_df, target, features, \n                 valid_df=None, categoricals=[],\n                 alpha=0, custom_loss=None,\n                 n_splits=3, cv_method=\"KFold\", group=None,\n                 task=\"regression\", params=None, parameter_tuning=False,\n                 seed=42, verbose=True):\n        self.train_df     = train_df\n        if valid_df is not None and valid_df.shape[0]==0:\n            self.valid_df = None            \n        else:\n            self.valid_df = valid_df\n        self.test_df      = test_df\n        self.target       = target\n        self.features     = features\n        self.n_splits     = n_splits\n        self.categoricals = categoricals\n        self.alpha        = alpha\n        self.custom_loss  = custom_loss\n        self.cv_method    = cv_method\n        self.group        = group\n        self.task         = task\n        self.parameter_tuning = parameter_tuning\n        self.seed    = seed\n        self.cv      = self.get_cv()\n        self.verbose = verbose\n        if params is None:\n            self.params = self.get_params()\n        else:\n            self.params = params\n        self.y_pred, self.y_valid, self.score, self.model, self.oof, self.y_val, self.fi_df = self.fit()\n\n    def train_model(self, train_set, val_set):\n        raise NotImplementedError\n\n    def get_params(self):\n        raise NotImplementedError\n\n    def convert_dataset(self, x_train, y_train, x_val, y_val):\n        raise NotImplementedError\n\n    def calc_metric(self, y_true, y_pred): # this may need to be changed based on the metric of interest\n        if   self.task in (\"multiclass\",\"nn_multiclass\"):\n            preds = np.argmax(y_pred, axis=1) if y_true.shape != y_pred.shape else y_pred\n            return f1_score(y_true, preds, average='macro')                \n        elif self.task == \"binary\":\n            return f1_score(y_true, y_pred, average='macro')\n        elif self.task in (\"regression\",\"quantile\"):\n            return np.sqrt(mean_squared_error(y_true, y_pred))\n        elif self.task == \"custom_loss\":\n            if y_pred.ndim==2:\n                return np.sqrt(mean_squared_error(y_true, y_pred[:,0]))\n            else:\n                return np.sqrt(mean_squared_error(y_true, y_pred))\n    \n    def get_cv(self):\n        if self.cv_method == \"KFold\":\n            cv = KFold(n_splits=self.n_splits, shuffle=True, random_state=self.seed)\n            return cv.split(self.train_df)\n        if self.cv_method == \"StratifiedKFold\":\n            cv = StratifiedKFold(n_splits=self.n_splits, shuffle=True, random_state=self.seed)\n            return cv.split(self.train_df, self.train_df[self.target])\n        if self.cv_method == \"TimeSeriesSplit\":\n            cv = TimeSeriesSplit(max_train_size=None, n_splits=self.n_splits)\n            return cv.split(self.train_df)\n        if self.cv_method == \"GroupKFold\":\n            if self.group in self.features:\n                self.features.remove(self.group)\n            cv = GroupKFold(n_splits=self.n_splits)\n            return cv.split(self.train_df[self.features], self.train_df[self.target], self.train_df[self.group])\n        if self.cv_method == \"GroupShuffleSplit\":\n            if self.group in self.features:\n                self.features.remove(self.group)\n            cv = GroupShuffleSplit(n_splits=self.n_splits, random_state=self.seed)\n            return cv.split(self.train_df[self.features], self.train_df[self.target], self.train_df[self.group])\n\n    def fit(self):\n        # Initialize\n        y_vals = np.zeros((self.train_df.shape[0], ))\n        if self.task in (\"multiclass\",\"nn_multiclass\"):\n            col_len = self.train_df[self.target].nunique()\n        elif self.task == \"custom_loss\":\n            col_len = 2\n        else:\n            col_len = 1\n        oof_pred = np.zeros((self.train_df.shape[0], col_len))\n        y_pred   = np.zeros((self.test_df.shape[0],  col_len))\n        y_valid  = np.zeros((self.valid_df.shape[0], col_len)) if self.valid_df is not None else None\n            \n        if self.group is not None:\n            if self.group in self.features:\n                self.features.remove(self.group)\n            if self.group in self.categoricals:\n                self.categoricals.remove(self.group)\n                \n        fi = np.zeros((self.n_splits, len(self.features)))\n        if y_valid is not None:\n            x_valid = self.valid_df[self.features].copy()\n            del self.valid_df\n            gc.collect()\n        x_test = self.test_df[self.features]\n\n        # Fitting with out of fold\n        for fold, (train_idx, val_idx) in enumerate(self.cv):\n            # Prepare train and test dataset\n            x_train = self.train_df.iloc[train_idx, :].query(\"last6\")[self.features]\n            y_train = self.train_df.iloc[train_idx, :].query(\"last6\")[self.target]\n            last3_val_idx = self.train_df.iloc[val_idx, :].query(\"last3\").index\n            x_val   = self.train_df.iloc[last3_val_idx, :][self.features]\n            y_val   = self.train_df.iloc[last3_val_idx, :][self.target]\n            train_set, val_set = self.convert_dataset(x_train, y_train, x_val, y_val)\n            del x_train, y_train\n            gc.collect()\n            \n            # Fit model\n            model, importance = self.train_model(train_set, val_set)\n            fi[fold, :]       = importance\n            y_vals[last3_val_idx] = y_val\n            \n            # Get some scores\n            oof_pred[last3_val_idx] = model.predict(x_val).reshape(oof_pred[last3_val_idx].shape)\n            if y_valid is not None:\n                y_valid += model.predict(x_valid).reshape(y_valid.shape) / self.n_splits\n            y_pred += model.predict(x_test).reshape(y_pred.shape) / self.n_splits\n                \n            print('Partial score of fold {} is: {}'.format(fold, self.calc_metric(y_val, oof_pred[last3_val_idx])))\n        \n        # Create feature importance data frame\n        fi_df = pd.DataFrame()\n        for n in np.arange(self.n_splits):\n            tmp = pd.DataFrame()\n            tmp[\"features\"]   = self.features\n            tmp[\"importance\"] = fi[n, :]\n            tmp[\"fold\"]       = n\n            fi_df = pd.concat([fi_df, tmp], ignore_index=True)\n        gfi   = fi_df[[\"features\", \"importance\"]].groupby([\"features\"]).mean().reset_index()\n        fi_df = fi_df.merge(gfi, on=\"features\", how=\"left\", suffixes=('', '_mean'))\n        \n        # Calculate oof score\n        \n        last3_idx  = np.where(oof_pred.sum(1)>0)\n        last3_pred = oof_pred[last3_idx]\n        last3_fvc  = y_vals[last3_idx]\n        loss_score = self.calc_metric(last3_fvc, last3_pred)\n        print('Our oof loss score is: ', loss_score)\n        \n        return y_pred, y_valid, loss_score, model, oof_pred, y_vals, fi_df\n\n    def plot_feature_importance(self, rank_range=[1, 50]):\n        fig, ax   = plt.subplots(1, 1, figsize=(10, 5))\n        sorted_df = self.fi_df.sort_values(by=\"importance_mean\", ascending=False).reset_index()\n        sns.barplot(data=sorted_df.iloc[self.n_splits*(rank_range[0]-1) : self.n_splits*rank_range[1]],\n                    x=\"importance\", y=\"features\", orient='h')\n        ax.set_xlabel(\"feature importance\")\n        ax.spines['top'].set_visible(False)\n        ax.spines['right'].set_visible(False)\n        return sorted_df\n    \nclass LgbModel(BaseModel):\n    \"\"\"\n    LGB wrapper\n    \"\"\"\n    def train_model(self, train_set, val_set):\n        verbosity = 100 if self.verbose else 0\n        if self.task == \"custom_loss\":\n            model = lgb.train(self.params, train_set, num_boost_round=10000,\n                              valid_sets=[train_set, val_set], verbose_eval=verbosity,\n                              fobj=self.custom_loss.return_grad_and_hess,\n                              feval=self.custom_loss.return_loss)\n        else:\n            model = lgb.train(self.params, train_set, num_boost_round=4000,\n                              valid_sets=[train_set, val_set], verbose_eval=verbosity)\n        fi        = model.feature_importance(importance_type=\"gain\")\n        return model, fi\n\n    def convert_dataset(self, x_train, y_train, x_val=None, y_val=None):\n        train_set   = lgb.Dataset(x_train, y_train, categorical_feature=self.categoricals)\n        if x_val is not None:\n            val_set = lgb.Dataset(x_val,   y_val,   categorical_feature=self.categoricals)\n            return train_set, val_set\n        return train_set\n\n    def get_params(self):\n        # Fast fit parameters\n        params = {\n            'boosting_type'    : \"gbdt\",\n            'objective'        : self.task,\n            \"subsample\"        : 0.4,\n            \"subsample_freq\"   : 1,\n            'max_depth'        : 1,\n            'min_data_in_leaf' : 50,\n            'learning_rate'    : 0.05,\n            'early_stopping_rounds' : 500,\n            'bagging_seed'     : 11,\n            'random_state'     : 42,\n            'verbosity'        : -1\n        }\n\n        # List is here: https://lightgbm.readthedocs.io/en/latest/Parameters.html\n        if   self.task == \"regression\":\n            params[\"metric\"]    = \"regression_l2\"\n        elif self.task == \"quantile\":\n            params[\"metric\"]    = \"quantile\"\n            params[\"alpha\"]     = self.alpha\n        elif self.task == \"custom_loss\":\n            del params[\"objective\"]\n            params[\"num_class\"] = 2\n            params[\"metric\"]    = 'None'\n        elif self.task == \"binary\":\n            params[\"metric\"]    = \"binary_logloss\"\n        elif self.task == \"multiclass\":\n            params[\"metric\"]    = \"multi_logloss\"\n            params[\"num_class\"] = len(self.train_df[self.target].unique())\n            \n        # Bayesian Optimization by Optuna\n        if self.parameter_tuning:\n            # Define objective function\n            def objective(trial):\n                # Split train and test data\n                train_x, test_x, train_y, test_y = train_test_split(self.train_df[self.features], \n                                                                    self.train_df[self.target],\n                                                                    test_size=0.3, random_state=self.seed)\n                \n                dtrain = lgb.Dataset(train_x, train_y, categorical_feature=self.categoricals)\n                dtest  = lgb.Dataset(test_x,  test_y,  categorical_feature=self.categoricals)\n\n                # Parameters to be explored\n                hyperparams = {'max_depth'         : trial.suggest_int('max_depth', 1, 1),\n                               \"subsample\"         : trial.suggest_uniform('subsample', 0.2, 0.5),\n                               \"subsample_freq\"    : trial.suggest_int('subsample_freq', 1, 1),\n                               'min_child_weight'  : trial.suggest_int('min_child_weight', 1, 20),\n                               'feature_fraction'  : trial.suggest_uniform('feature_fraction', 0.4, 1.0),\n                               'bagging_fraction'  : trial.suggest_uniform('bagging_fraction', 0.4, 1.0),\n                               'bagging_freq'      : trial.suggest_int('bagging_freq', 1, 7),\n                               'min_child_samples' : trial.suggest_int('min_child_samples', 5, 100),\n                               'lambda_l1'         : trial.suggest_loguniform('lambda_l1', 1e-8, 10.0),\n                               'lambda_l2'         : trial.suggest_loguniform('lambda_l2', 1e-8, 10.0),\n                               'early_stopping_rounds' : 500}\n                # LGBM\n                params.update(hyperparams)\n                verbosity = 500 if self.verbose else 0\n                if self.task == \"custom_loss\":\n                    model = lgb.train(params, dtrain, valid_sets=dtest, verbose_eval=verbosity,\n                                      fobj=self.custom_loss.return_grad_and_hess,\n                                      feval=self.custom_loss.return_loss)\n                else:\n                    model = lgb.train(params, dtrain, valid_sets=dtest, verbose_eval=verbosity)\n                pred  = model.predict(test_x)\n                return self.calc_metric(test_y, pred)\n\n            # Run optimization\n            study = optuna.create_study(direction='minimize')\n            study.optimize(objective, n_trials=50)\n            print('Number of finished trials: {}'.format(len(study.trials)))\n            print('Best trial:')\n            trial = study.best_trial\n            print('  Value: {}'.format(trial.value))\n            print('  Params: ')\n            for key, value in trial.params.items():\n                print('    {}: {}'.format(key, value))\n\n            params.update(trial.params)\n            #params[\"learning_rate\"] = 0.001\n            # Plot history\n            plot_optimization_history(study)\n            \n        return params","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"exec_custom_loss = True\nexec_regression  = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if exec_custom_loss:\n    lgbm_p1_c = LgbModel(train, sub, \"FVC\", FE,\n                         valid_df=None, task=\"custom_loss\", custom_loss=OSICLossForLGBM(),\n                         cv_method=\"GroupKFold\", n_splits=5, group=\"Patient\", \n                         params=None, parameter_tuning=False, verbose=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if exec_custom_loss:\n    df_imp = lgbm_p1_c.plot_feature_importance()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if exec_regression:\n    lgbm_p1 = LgbModel(train, sub, \"FVC\", FE,\n                       valid_df=None, task=\"regression\",\n                       cv_method=\"GroupKFold\", n_splits=5, group=\"Patient\", \n                       params=None, parameter_tuning=True, verbose=False)    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if exec_regression:\n    df_imp = lgbm_p1.plot_feature_importance()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if exec_regression:\n    lgbm_p0 = LgbModel(train, sub, \"FVC\", FE,\n                       valid_df=None, task=\"quantile\", alpha=0.25,\n                       cv_method=\"GroupKFold\", n_splits=5, group=\"Patient\", \n                       params=None, parameter_tuning=True, verbose=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if exec_regression:\n    lgbm_p2 = LgbModel(train, sub, \"FVC\", FE,\n                       valid_df=None, task=\"quantile\", alpha=0.75,\n                       cv_method=\"GroupKFold\", n_splits=5, group=\"Patient\", \n                       params=None, parameter_tuning=True, verbose=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def change_values(idx, ary):\n    tmp_min = ary[idx].min()\n    tmp_med = np.median(ary[idx])\n    tmp_max = ary[idx].max()\n    ary[idx,0] = tmp_min\n    ary[idx,1] = tmp_med\n    ary[idx,2] = tmp_max\n    print(\"pred_lgbm\", ary[idx])\n    return ary\n\nif exec_regression:\n    pred_lgbm = np.hstack([lgbm_p0.oof,    lgbm_p1.oof,    lgbm_p2.oof])\n    pe_lgbm   = np.hstack([lgbm_p0.y_pred, lgbm_p1.y_pred, lgbm_p2.y_pred])\n\n    for idx in np.where((pred_lgbm[:,2]-pred_lgbm[:,0])<0)[0]:\n        pred_lgbm = change_values(idx, pred_lgbm)\n    for idx in np.where((pe_lgbm[:,2]-pe_lgbm[:,0])<0)[0]:\n        pe_lgbm   = change_values(idx, pe_lgbm)\n        \nif exec_custom_loss:\n    pred_lgbm_c = np.vstack([lgbm_p1_c.oof[:,0]-(lgbm_p1_c.oof[:,1]/2),\n                             lgbm_p1_c.oof[:,0],\n                             lgbm_p1_c.oof[:,0]+(lgbm_p1_c.oof[:,1]/2)]).T\n    pe_lgbm_c   = np.vstack([lgbm_p1_c.y_pred[:,0]-(lgbm_p1_c.y_pred[:,1]/2),\n                             lgbm_p1_c.y_pred[:,0],\n                             lgbm_p1_c.y_pred[:,0]+(lgbm_p1_c.y_pred[:,1]/2)]).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emsambles = []\n\nif exec_regression and exec_custom_loss:\n    for a in np.arange(0,1.1,0.1):\n        a = round(a,1)\n        a_max = 1-a\n        for b in np.arange(0,a_max+0.1,0.1):\n            b = round(b,1)\n            c = round((a_max-b),1)\n            emsambled  = a*pred_lgbm + b*pred_lgbm_c + c*pred\n            last3_pred = emsambled[np.where(emsambled.sum(1))]\n            last3_fvc  = train.query(\"last3\").FVC\n            lll = calculate_scores(last3_fvc, last3_pred, True, False)\n            emsambles.append([a,b,c,lll])\n    df_result_pred = pd.DataFrame(emsambles, columns=[\"a\",\"b\",\"c\",\"lll\"])\n    a,b,c = df_result_pred.sort_values(by=\"lll\", ascending=False).iloc[0,:3]\n    best_pred = a*pred_lgbm + b*pred_lgbm_c + c*pred\n    print(\"Best emsamble parameters are %s,%s,%s\" % (a,b,c))\n\nelse:\n    for a in np.arange(0,1.1,0.1):\n        a = round(a,1)\n        print(a)\n        emsambled  = a*pred_lgbm + (1-a)*pred if exec_regression else a*pred_lgbm_c + (1-a)*pred\n        last3_pred = emsambled[np.where(emsambled.sum(1))]\n        last3_fvc  = train.query(\"last3\").FVC\n        lll = calculate_scores(last3_fvc, last3_pred, True)\n        emsambles.append([a,lll])\n    \n    a = emsambles[np.argmax(np.array(emsambles)[:,1])][0]\n    best_pred = a*pred_lgbm + (1-a)*pred if exec_regression else a*pred_lgbm_c + (1-a)*pred\n    print(\"Best emsamble parameter is %s\" % a)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"last3_pred = best_pred[np.where(best_pred.sum(1))]\nlast3_fvc  = train.query(\"last3\").FVC\ncalculate_scores(last3_fvc, last3_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,5))\nidxs = np.random.randint(0, last3_pred.shape[0], 200)\nplt.plot(last3_fvc.values[idxs], label=\"ground truth\")\nplt.plot(last3_pred[idxs, 0], label=\"q25\")\nplt.plot(last3_pred[idxs, 1], label=\"q50\")\nplt.plot(last3_pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(last3_pred[:, 2]-last3_pred[:, 0])\nplt.title(\"uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# PREDICTION\nif exec_regression and exec_custom_loss:\n    sub['FVC1']        = a*pe_lgbm[:, 1] + b*pe_lgbm_c[:, 1] + c*pe[:, 1]\n    sub['Confidence1'] = a*(pe_lgbm[:, 2] - pe_lgbm[:, 0]) + b*(pe_lgbm_c[:, 2] - pe_lgbm_c[:, 0]) + c*(pe[:, 2] - pe[:, 0])\nelif exec_regression:\n    sub['FVC1']        = a*pe_lgbm[:, 1] + (1-a)*pe[:, 1]\n    sub['Confidence1'] = a*(pe_lgbm[:, 2] - pe_lgbm[:, 0]) + (1-a)*(pe[:, 2] - pe[:, 0])\nelse:\n    sub['FVC1']        = a*pe_lgbm_c[:, 1] + (1-a)*pe[:, 1]\n    sub['Confidence1'] = a*(pe_lgbm_c[:, 2] - pe_lgbm_c[:, 0]) + (1-a)*(pe[:, 2] - pe[:, 0])\n\nsubm = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()\nsubm.loc[~subm.FVC1.isnull(),'FVC'] = subm.loc[~subm.FVC1.isnull(),'FVC1']\nsubm.loc[~subm.FVC1.isnull(),'Confidence'] = subm.loc[~subm.FVC1.isnull(),'Confidence1']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(subm.head())\nsubm[[\"FVC\",\"Confidence\"]].describe().T","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_regression.csv\", index=False)\n\nreg_sub = subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Ensemble (Simple Blend)"},{"metadata":{"trusted":true},"cell_type":"code","source":"if USE_ALL:\n    emsambles = []\n    for a in np.arange(0,1.1,0.05):\n        a = round(a,2)\n        emsambled = a*last3_pred_img + (1-a)*last3_pred\n        lll = calculate_scores(last3_fvc, emsambled, True, False)\n        emsambles.append([a,lll])\n    best_a = emsambles[np.argmax(np.array(emsambles)[:,1])][0]\n    best_pred = best_a*last3_pred_img + (1-best_a)*last3_pred\n    print(\"Best emsamble parameter is %s\" % best_a)\n    calculate_scores(last3_fvc, best_pred, show=True)\nelse:\n    best_a = 0.3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df1 = img_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)\ndf2 = reg_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df1[['Patient_Week']].copy()\ndf['FVC']        = best_a*df1['FVC']        + (1-best_a)*df2['FVC']\ndf['Confidence'] = best_a*df1['Confidence'] + (1-best_a)*df2['Confidence']\nprint(df.shape)\nsample_submission = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\nsubmission = sample_submission[[\"Patient_Week\"]].merge(df, on=\"Patient_Week\")\nprint(submission.shape)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.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}