{"cells":[{"metadata":{},"cell_type":"markdown","source":"## **Quantile Regression with CT + tabular features [pytorch]**\n\n<p align=\"center\">\n  <img src=\"https://i0.wp.com/post.medicalnewstoday.com/wp-content/uploads/sites/3/2020/08/What-is-idiopathic-pulmonary-fibrosis_header-1024x575.jpg?w=1155&h=1528\" />\n</p>\n\n\n### **What is IPF?**\n\nIdiopathic pulmonary fibrosis (IPF) is a chronic and progressive lung condition. IPF causes scar tissue to build up in the lungs, leading to shortness of breath and a persistent cough.\n\nPulmonary fibrosis is a scarring of the lungs.\n\nThe NHLBI state that IPF occurs when the lung tissue that surrounds the air sacs thickens and becomes stiff, leading to permanent scarring of the lungs.\n\nThis scarring affects the function of the lungs, making it difficult to breathe.\n\n- https://www.medicalnewstoday.com/articles/idiopathic-pulmonary-fibrosis#definition\n\n### **What are we doing?**\n\n\nIn this competiton, we're given a CT scan (which contains different number of slices), multiple observation of patients' FVC with tabular features. We need to predict past and future FVC (Forced vital capacity) from that data.\n\n### **How are we doing it?**\n\n\n1. We're going to extract useful features from the CT scans (volume, histogram mean, skew, kurtosis).\n2. We'll merge the CT features with tabular features.\n3. We'll use quantile regression model to predict the three percentiles with quantile loss function that'll give us the FVC and confidence.\n","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# necessary imports\nimport os\nimport cv2\nimport random\nimport pickle\nimport glob\nfrom pathlib import Path\nimport gc\nfrom tqdm.notebook import tqdm\nimport pandas as pd\nimport numpy as np\nfrom PIL import Image\nfrom tqdm.auto import tqdm\n\nimport pydicom\nfrom pydicom.tag import Tag\n\nimport scipy.ndimage as ndimage\nfrom scipy.ndimage import zoom\nfrom scipy.stats import kurtosis\nfrom scipy.stats import skew\n\nfrom sklearn.cluster import KMeans\n\nfrom skimage import measure, morphology, segmentation","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Everything we need to extract features from CT slices.","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# aggregation of ct features and tabular data\n\ndef load_scan(path):\n    \"\"\"\n    Loads scans from a folder and into a list.\n    \n    Parameters: path (Folder path)\n    \n    Returns: slices (List of slices)\n    \"\"\"\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    \n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        try:\n            slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        except:\n            try:\n                slice_thickness = slices[0].SliceThickness\n            except:\n                slice_thickness = 0.5 # random value for now\n            \n    # print(slice_thickness)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\ndef transform_to_hu(slices):\n    \"\"\"\n    transform dicom.pixel_array to Hounsfield.\n    Parameters: list dicoms\n    Returns:numpy Hounsfield\n    \"\"\"\n    \n    images = np.stack([file.pixel_array for file in slices])\n    images = images.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    for n in range(len(slices)):\n        \n        intercept = slices[n].RescaleIntercept\n        slope = slices[n].RescaleSlope\n        \n        if slope != 1:\n            images[n] = slope * images[n].astype(np.float64)\n            images[n] = images[n].astype(np.int16)\n            \n        images[n] += np.int16(intercept)\n    \n    return np.array(images, dtype=np.int16)\n\ndef generate_internal_mask(image):\n    \"\"\"\n    Generates markers for a given image.\n    Parameters: image\n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    \n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    \n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    \n    marker_internal = marker_internal_labels > 0\n    \n    return marker_internal\n\ndef generate_markers(image):\n    \"\"\"\n    Generates markers for a given image.\n    Parameters: image\n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    \n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    \n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    \n    marker_internal = marker_internal_labels > 0\n    \n    # Creation of the External Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    \n    # Creation of the Watershed Marker\n    marker_watershed = np.zeros((image.shape[0], image.shape[1]), dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed\n\ndef seperate_lungs_Watershed(image, iterations = 1):\n    \"\"\"\n    Segments lungs using various techniques.\n    \n    Parameters: image (Scan image), iterations (more iterations, more accurate mask)\n    \n    Returns: \n        - Segmented Lung\n        - Lung Filter\n        - Outline Lung\n        - Watershed Lung\n        - Sobel Gradient\n    \"\"\"\n    \n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n    \n    \n    '''\n    Creation of Sobel Gradient\n    '''\n    \n    # Sobel-Gradient\n    sobel_filtered_dx = ndimage.sobel(image, 1)\n    sobel_filtered_dy = ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n    \n    \n    '''\n    Using the watershed algorithm\n    \n    \n    We pass the image convoluted by sobel operator and the watershed marker\n    to morphology.watershed and get a matrix matrix labeled using the \n    watershed segmentation algorithm.\n    '''\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    \n    \n    return watershed\n\ndef crop_image(img: np.ndarray):\n    edge_pixel_value = img[0, 0]\n    mask = img != edge_pixel_value\n    return img[np.ix_(mask.any(1),mask.any(0))]\n\ndef resize_image(img: np.ndarray,reshape=(512,512)):\n    img = [cv2.resize(im,(512,512)) for im in img]\n    return img\n\ndef preprocess_img(img,local_pd):\n    #if local_pd.resize_type == 'resize':\n    #    img = [resize_image(im) for im in img]\n    if local_pd.resize_type == 'crop':\n        img = [crop_image(im) for im in img]\n        \n    return np.array(img, dtype=np.int16)\n\ndef func_volume(patient_scan,patient_mask):\n    pixel_spacing = patient_scan.PixelSpacing\n    slice_thickness = patient_scan.SliceThickness\n    slice_volume = np.count_nonzero(patient_mask)*pixel_spacing[0]*pixel_spacing[1]*slice_thickness\n    return slice_volume\n\ndef caculate_lung_volume(patient_scans,patient_images):\n    \"\"\"\n    caculate volume of lung from mask\n    Parameters: list dicom scans,list patient CT image\n    Returns: volume cm³　(float)\n    \"\"\"\n    # patient_masks = pool.map(generate_internal_mask,patient_images)\n    patient_masks = list(map(generate_internal_mask,patient_images)) # non-parallel version\n    lung_volume = np.array(list(map(func_volume,patient_scans,patient_masks))).sum()\n\n    return lung_volume*0.001\n\ndef caculate_histgram_statistical(patient_images,thresh = [-500,-50]):\n    \"\"\"\n    caculate hisgram kurthosis of lung hounsfield\n    Parameters: list patient CT image 512*512,thresh divide lung\n    Returns: histgram statistical characteristic(Mean,Skew,Kurthosis)\n    \"\"\"\n    statistical_characteristic = dict(Mean=0,Skew=0,Kurthosis=0)\n    num_slices = len(patient_images)\n    patient_images = patient_images[int(num_slices*0.1):int(num_slices*0.9)]\n    patient_images_mean = np.mean(patient_images,0)\n    \n    s_pixel = patient_images_mean.flatten()\n    s_pixel = s_pixel[np.where((s_pixel)>thresh[0]&(s_pixel<thresh[1]))]\n    \n    statistical_characteristic['Mean'] = np.mean(s_pixel)\n    statistical_characteristic['Skew'] = skew(s_pixel)\n    statistical_characteristic['Kurthosis'] = kurtosis(s_pixel)\n    \n    return statistical_characteristic\n\ndef aggregate_dicom_features(dicom_path = 'data_download/train/', csv_path = 'data_download/train.csv', show = True): # pass path to datafolder\n    \n    Patients_id = os.listdir(dicom_path)\n\n    n_dicom_dict = {\"Patient\":[],\"n_dicom\":[],\"list_dicom\":[]}\n\n    for Patient_id in Patients_id:\n        dicom_id_path = glob.glob(dicom_path + Patient_id + \"/*\")\n        n_dicom_dict[\"n_dicom\"].append(len(dicom_id_path))\n        n_dicom_dict[\"Patient\"].append(Patient_id)\n        list_dicom_id = sorted([int(i.split(\"/\")[-1][:-4]) for i in dicom_id_path])\n        n_dicom_dict[\"list_dicom\"].append(list_dicom_id)\n\n    dicom_pd = pd.DataFrame(n_dicom_dict)\n    \n    if show:\n        dicom_pd.head()\n    \n    dicom_pd['height'],dicom_pd['width'] = -1,-1\n    for Patient_id in Patients_id:\n        dicom_id_path = glob.glob(dicom_path + Patient_id + \"/*\")\n        for patient_dicom_id_path in dicom_id_path:\n            try:\n                dicom = pydicom.dcmread(patient_dicom_id_path)\n                dicom_pd.loc[dicom_pd.Patient==Patient_id,'height'] = dicom.Rows\n                dicom_pd.loc[dicom_pd.Patient==Patient_id,'width'] = dicom.Columns\n            except:\n                dicom_pd.loc[dicom_pd.Patient==Patient_id,'height'] = 512\n                dicom_pd.loc[dicom_pd.Patient==Patient_id,'width'] = 512\n                \n            \n            \n    reshape_dicom_pd = dicom_pd.loc[(dicom_pd.height!=512) | (dicom_pd.width!=512),:]\n    reshape_dicom_pd = reshape_dicom_pd.reset_index(drop=True)\n    \n    if show:\n        reshape_dicom_pd.head()\n    \n    crop_id = list(reshape_dicom_pd[reshape_dicom_pd['height']!=reshape_dicom_pd['width']]['Patient'])\n    reshape_dicom_pd['resize_type'] = 'resize'\n    reshape_dicom_pd.loc[reshape_dicom_pd.Patient.isin(crop_id),'resize_type'] = 'crop'\n    \n    dicom_pd['resize_type'] = 'no'\n    for idx,i in enumerate(reshape_dicom_pd['Patient']):\n        dicom_pd.loc[dicom_pd.Patient==i,'resize_type'] = reshape_dicom_pd.loc[idx,'resize_type']\n    \n    if show:\n        dicom_pd.head()\n    \n    train_pd = pd.read_csv(csv_path)\n    temp_pd = pd.DataFrame(columns=train_pd.columns)\n    \n    for i in range(len(dicom_pd)):\n        patient_pd = train_pd[train_pd.Patient==dicom_pd.iloc[i].Patient]\n        zeroweek = patient_pd['Weeks'].min()\n        if sum(patient_pd.Weeks==zeroweek)>1:\n            print(pd.unique(patient_pd.Patient))\n        temp_pd = temp_pd.append(patient_pd[patient_pd.Weeks==zeroweek].iloc[0])\n    dicom_pd = pd.merge(dicom_pd, temp_pd, on=['Patient'])\n    dicom_pd.head()\n    \n    lung_stat_pd = pd.DataFrame(columns=['Patient','Volume','Mean','Skew','Kurthosis'])\n\n\n    for i in tqdm(range(len(dicom_pd))):\n        path = os.path.join(dicom_path,dicom_pd.iloc[i].Patient)\n        lung_stat_pd.loc[i,'Patient'] = dicom_pd.iloc[i].Patient\n        try:\n            patient_scans = load_scan(path)\n            patient_images = transform_to_hu(patient_scans)\n            patient_images = preprocess_img(patient_images,dicom_pd.iloc[i])\n            lung_stat_pd.loc[i,'Volume'] = caculate_lung_volume(patient_scans,patient_images)\n\n            patient_images = resize_image(patient_images) if dicom_pd.iloc[i].resize_type=='resize' else patient_images\n\n            statistical_characteristic = caculate_histgram_statistical(patient_images)\n            lung_stat_pd.loc[i,'Mean'] = statistical_characteristic['Mean']\n            lung_stat_pd.loc[i,'Skew'] = statistical_characteristic['Skew']\n            lung_stat_pd.loc[i,'Kurthosis'] = statistical_characteristic['Kurthosis']\n        except:\n            # fill up with mean features\n            lung_stat_pd.loc[i,'Volume'] = 3554.55 # training mean\n            lung_stat_pd.loc[i,'Mean'] = 81.75 \n            lung_stat_pd.loc[i,'Skew'] = 1.90\n            lung_stat_pd.loc[i,'Kurthosis'] = 5.69\n            \n            \n\n    \n    dicom_feature = pd.merge(dicom_pd, lung_stat_pd, on=['Patient'])\n    \n    dicom_feature = dicom_feature.drop(['list_dicom', 'height','width','resize_type','n_dicom'], axis=1)\n    \n    dicom_feature.head()\n    \n    return dicom_feature","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Using pre-computed features as it's slow.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# change in inference\nuse_saved = True\n\nfolder_path = '../input/osic-pulmonary-fibrosis-progression'\ndicom_path = f'{folder_path}/train/'\ncsv_path = f'{folder_path}/train.csv'\n\n# dummy data == real data, when only tabular\ntrain_csv = f'{folder_path}/train.csv'\nmodel_data_path = \"processed-features\" # may be changed\n\ntrain_data_tab = pd.read_csv(train_csv)\n\nif not use_saved:\n    merged_data = aggregate_dicom_features(dicom_path = dicom_path, csv_path = csv_path, show = True)\nelse:\n    merged_data = pd.read_csv(f'../input/{model_data_path}/train_data_ct_tab.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_tab.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Let's merge features from CT scans and tabular ones into a single dataframe.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# merging both\ntrain_data_tab['Volume'] = merged_data['Volume'].mean()\ntrain_data_tab['Mean'] = merged_data['Mean'].mean()\ntrain_data_tab['Skew'] = merged_data['Skew'].mean()\ntrain_data_tab['Kurthosis'] = merged_data['Kurthosis'].mean()\n\nfor i in range(len(train_data_tab)):\n    pid = train_data_tab.iloc[i]['Patient']\n    train_data_tab.at[i, 'Volume'] = merged_data[merged_data['Patient']==pid].iloc[0]['Volume']\n    train_data_tab.at[i, 'Mean'] = merged_data[merged_data['Patient']==pid].iloc[0]['Mean']\n    train_data_tab.at[i, 'Skew'] = merged_data[merged_data['Patient']==pid].iloc[0]['Skew']\n    train_data_tab.at[i, 'Kurthosis'] = merged_data[merged_data['Patient']==pid].iloc[0]['Kurthosis']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_tab.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Adding baseline week into the feature to let the model know about the offset (in weeks).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#getting base week for patient\ndef get_baseline_week(data):\n    df = data.copy()\n    df['Weeks'] = df['Weeks'].astype(int)\n    df['min_week'] = df.groupby('Patient')['Weeks'].transform('min')\n    df['baseline_week'] = df['Weeks'] - df['min_week']\n    return df\n\n#getting FVC for base week and setting it as bas_FVC of patient\ndef get_base_FVC(data):\n    df = data.copy()\n    base = df.loc[df.Weeks == df.min_week][['Patient','FVC']].copy()\n    base.columns = ['Patient','base_FVC']\n    \n    base['nb']=1\n    base['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\n    \n    base = base[base.nb==1]\n    base.drop('nb',axis =1,inplace=True)\n    df = df.merge(base,on=\"Patient\",how='left')\n    df.drop(['min_week'], axis = 1)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler,LabelEncoder,OneHotEncoder,PowerTransformer\nfrom sklearn.model_selection import train_test_split, cross_val_score,cross_validate, KFold\nfrom sklearn.metrics import make_scorer\nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.pipeline import Pipeline","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Now, it's time to load the test data, aggregate features.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# inference\nfolder_path = '../input/osic-pulmonary-fibrosis-progression'\ndicom_path = f'{folder_path}/test/'\ncsv_path = f'{folder_path}/test.csv'\ntest_merged = aggregate_dicom_features(dicom_path = dicom_path, csv_path = csv_path, show = True)\ntest_data_tab = pd.read_csv(csv_path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data_tab.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# merging both\ntest_data_tab['Volume'] = test_merged['Volume'].mean()\ntest_data_tab['Mean'] = test_merged['Mean'].mean()\ntest_data_tab['Skew'] = test_merged['Skew'].mean()\ntest_data_tab['Kurthosis'] = test_merged['Kurthosis'].mean()\n\nfor i in range(len(test_data_tab)):\n    pid = test_data_tab.iloc[i]['Patient']\n    try:\n        test_data_tab.at[i, 'Volume'] = test_merged[test_merged['Patient']==pid].iloc[0]['Volume']\n        test_data_tab.at[i, 'Mean'] = test_merged[test_merged['Patient']==pid].iloc[0]['Mean']\n        test_data_tab.at[i, 'Skew'] = test_merged[test_merged['Patient']==pid].iloc[0]['Skew']\n        test_data_tab.at[i, 'Kurthosis'] = test_merged[test_merged['Patient']==pid].iloc[0]['Kurthosis']\n    except:\n        print('failed at here: 1')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data_tab.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# alisaing\ntrain_data = train_data_tab\ntest_data = test_data_tab\n\nsample_csv = f'{folder_path}/sample_submission.csv'\n\ntrain_data.drop_duplicates(keep=False,inplace=True,subset=['Patient','Weeks'])\ntrain_data = get_baseline_week(train_data)\ntrain_data = get_base_FVC(train_data)\n\nsample = pd.read_csv(sample_csv)\nsample.drop(\"FVC\",axis=1,inplace=True)\nsample[[\"Patient\",\"Weeks\"]] = sample[\"Patient_Week\"].str.split(\"_\",expand=True) \nsample = sample.merge(test_data.drop(\"Weeks\",axis=1),on=\"Patient\",how=\"left\")\n\n#we have to predict for all weeks \nsample[\"min_Weeks\"] = np.nan\nsample = get_baseline_week(sample)\nsample = get_base_FVC(sample)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_columns = ['baseline_week','base_FVC','Percent','Age','Sex','SmokingStatus', 'Volume', 'Mean', 'Skew', 'Kurthosis'] # more features\ntrain_label = ['FVC']\nsub_columns = ['Patient_Week','FVC','Confidence']\n\ntrain = train_data[train_columns]\ntest = sample[train_columns]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# pre-processing, need to improve this\ntransformer = ColumnTransformer([('s',StandardScaler(),[0,1,2,3,6,7,8,9]),('o',OneHotEncoder(),[4,5])])\ntarget = train_data[train_label].values\ntrain = transformer.fit_transform(train)\ntest = transformer.transform(test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Now, a simple pytorch model with quantile loss (+ a weighted version of the metric).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.functional as F\nimport torch.optim as optim\nfrom torch.utils.data import TensorDataset, DataLoader\n\nclass Model(nn.Module):\n    def __init__(self,n):\n        super(Model,self).__init__()\n        self.layer1 = nn.Linear(n,128)\n        self.relu1 = nn.LeakyReLU()\n        self.layer2 = nn.Linear(128,128)\n        self.relu2 = nn.LeakyReLU()\n        #self.layer3 = nn.Linear(64,32)\n        #self.relu3 = nn.LeakyReLU()\n        \n        self.out1 = nn.Linear(128,3) # 3 outputs\n            \n    def forward(self,xb):\n        x1 = self.relu1(self.layer1(xb))\n        x1 = self.relu2(self.layer2(x1))\n        # x1 = self.relu3(self.layer3(x1))\n        \n        o1 = self.out1(x1)\n        return o1 ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### We'll be saving the best models with k-fold validation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"best_metrics_so_far = [99999999]*10 # the smaller the better, at most 10 folds\n\ndef run():\n    \n    def score(outputs,target):\n        confidence = outputs[:,2] - outputs[:,0] # output -> 0-> 20% percentile, 50% percentile, 80% percentile\n        clip = torch.clamp(confidence,min=70)\n        target=torch.reshape(target,outputs[:,1].shape)\n        delta = torch.abs(outputs[:,1] - target)\n        delta = torch.clamp(delta,max=1000)\n        sqrt_2 = torch.sqrt(torch.tensor([2.])).to(device)\n        metrics = (delta*sqrt_2/clip) + torch.log(clip*sqrt_2)\n        return torch.mean(metrics)\n    \n    def qloss(outputs,target):\n        qs = [0.2,0.5,0.8]\n        qs = torch.tensor(qs,dtype=torch.float).to(device)\n        e =  outputs - target\n        e.to(device)\n        v = torch.max(qs*e,(qs-1)*e)\n        return torch.mean(v)\n\n    \n    def loss_fn(outputs,target,l):\n        return l * qloss(outputs,target) + (1- l) * score(outputs,target)\n        \n    def train_loop(train_loader,model,loss_fn,device,optimizer,lr_scheduler=None):\n        model.train()\n        losses = list()\n        metrics = list()\n        for i, (inputs, labels) in enumerate(train_loader):\n            inputs = inputs.to(device)\n            labels = labels.to(device)\n            \n            optimizer.zero_grad()\n            with torch.set_grad_enabled(True):           \n                outputs = model(inputs)                 \n                metric = score(outputs,labels)\n\n                loss = loss_fn(outputs,labels,0.77)\n                metrics.append(metric.cpu().detach().numpy())\n                losses.append(loss.cpu().detach().numpy())\n\n                loss.backward()\n\n                optimizer.step()\n                if lr_scheduler != None:\n                    lr_scheduler.step()\n            \n        return losses,metrics\n    \n    def valid_loop(k,valid_loader,model,loss_fn,device):\n        model.eval()\n        losses = list()\n        metrics = list()\n        for i, (inputs, labels) in enumerate(valid_loader):\n            inputs = inputs.to(device)\n            labels = labels.to(device)\n            \n            outputs = model(inputs)                 \n            metric = score(outputs,labels)\n            \n            loss = loss_fn(outputs,labels,0.77)\n            metrics.append(metric.cpu().detach().numpy())\n            losses.append(loss.cpu().detach().numpy())\n        \n        # print(np.mean(metrics))\n        global best_metrics_so_far\n        if np.mean(metrics) <= best_metrics_so_far[k]:\n            torch.save(model.state_dict(),f'model{k}.bin')\n            print('model saved')\n            best_metrics_so_far[k] = np.mean(metrics)\n            \n        return losses,metrics    \n\n    NFOLDS = 2\n    kfold = KFold(NFOLDS,shuffle=True,random_state=1997)\n    \n    #kfold\n    for k , (train_idx,valid_idx) in enumerate(kfold.split(train)):\n        batch_size = 128\n        epochs = 666\n        device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n        print(f\"{device} is used\")\n        x_train,x_valid,y_train,y_valid = train[train_idx,:],train[valid_idx,:],target[train_idx],target[valid_idx]\n        n = x_train.shape[1]\n        model = Model(n)\n        model.to(device)\n        lr = 0.2\n        optimizer = optim.Adam(model.parameters(),lr=lr)\n        lr_scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=25, gamma=0.44)\n\n        train_tensor = torch.tensor(x_train,dtype=torch.float)\n        y_train_tensor = torch.tensor(y_train,dtype=torch.float)\n\n        train_ds = TensorDataset(train_tensor,y_train_tensor)\n        train_dl = DataLoader(train_ds,\n                             batch_size = batch_size,\n                             num_workers=4,\n                             shuffle=True\n                             )\n\n        valid_tensor = torch.tensor(x_valid,dtype=torch.float)\n        y_valid_tensor = torch.tensor(y_valid,dtype=torch.float)\n\n        valid_ds = TensorDataset(valid_tensor,y_valid_tensor)\n        valid_dl = DataLoader(valid_ds,\n                             batch_size = batch_size,\n                             num_workers=4,\n                             shuffle=False\n                             )\n        \n        print(f\"Fold {k}\")\n        for i in range(epochs):\n            losses,metrics = train_loop(train_dl,model,loss_fn,device,optimizer,lr_scheduler)\n            valid_losses,valid_metrics = valid_loop(k,valid_dl,model,loss_fn,device)\n            if (i+1)%50==0:\n                print(f\"epoch:{i} Training | loss:{np.mean(losses)} score: {np.mean(metrics)}| \\n Validation | loss:{np.mean(valid_losses)} score:{np.mean(valid_metrics)}|\")\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"run()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def inference(model_path = ''):\n    device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n    nfold = 2\n    all_prediction = np.zeros((test.shape[0],3))\n    \n    model_weights = [0.5, 0.5]\n    for k in range(nfold):\n        n = train.shape[1]\n        \n        model = Model(n)\n        model.load_state_dict(torch.load(f\"model{k}.bin\"))\n        predictions = list()\n        model.to(device)\n        test_tensor = torch.tensor(test,dtype=torch.float)\n        test_dl = DataLoader(test_tensor,\n                        batch_size=64,\n                        num_workers=2,\n                        shuffle=False)\n    \n        with torch.no_grad():\n            for i, inputs in enumerate(test_dl):\n                inputs = inputs.to(device, dtype=torch.float)\n                outputs= model(inputs) \n                predictions.extend(outputs.cpu().detach().numpy())\n\n        all_prediction += np.array(predictions) * model_weights[k]\n        \n    return all_prediction  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction = inference()\nsample[\"Confidence\"] = np.abs(prediction[:,2] - prediction[:,0])\nsample[\"FVC\"] = prediction[:,1]\nsub = sample[sub_columns]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plotting\nimport matplotlib.pyplot as plt\nimport numpy as np\n\ndef plot_prediction(pid = 'ID00419637202311204720264'):\n\n    fvc = np.array(sub[sub['Patient_Week'].str.contains(pid)]['FVC'])\n    conf = np.array(sub[sub['Patient_Week'].str.contains(pid)]['Confidence'])\n\n    month = np.array(range(-12,134))\n\n    fig, ax = plt.subplots()\n    ax.plot(month,fvc)\n    ax.fill_between(month, (fvc-conf/2), (fvc+conf/2), color='b', alpha=.1)\n    ax.set_title(f'PID: {pid}')\n    ax.set_xlabel('month')\n    ax.set_ylabel('FVC')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# let's see our predictions\npid_list = ['ID00419637202311204720264', 'ID00421637202311550012437', 'ID00422637202311677017371', 'ID00423637202312137826377', 'ID00426637202313170790466']\n\nfor pid in pid_list:\n    plot_prediction(pid)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.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}