{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook is based on @ulrich07 https://www.kaggle.com/ulrich07/osic-multiple-quantile-regression-starter and code with pytorch.If you are not familiar with tensorflow, refer to this notebook.\nbtw,Pytorch is a bit slow than tensorflow."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nimport glob\nimport sys\nimport cv2\n\nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold\n#import gdcm\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sys.path.append('../input/efficientnet-pytorch/EfficientNet-PyTorch-master')\nsys.path.append('../input/pretrainedmodels/pretrainedmodels-0.7.4/')\nsys.path.append('../input/segmentation-models-pytorch/')\nimport segmentation_models_pytorch as smp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch.nn as nn\nimport torch\nfrom torch.utils.data.dataset import Dataset\nfrom torch.utils.data import DataLoader\n\nfrom scipy.ndimage import zoom\nfrom scipy.stats import kurtosis\nfrom scipy.stats import skew","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=2020):\n    random.seed(seed)\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    if torch.cuda.is_available():\n        torch.cuda.manual_seed(seed)\n        torch.cuda.manual_seed_all(seed)#set all gpus seed\n        torch.backends.cudnn.deterministic = True\n        torch.backends.cudnn.benchmark = False#if input data type and channels' changes arent' large use it improve train efficient\n        torch.backends.cudnn.enabled = True\n    \nseed_everything(42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\ndevice = torch.device('cuda')\n\ntr = pd.read_csv(f\"{ROOT}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nfeature_ct = pd.read_csv('../input/2020osic/CT_21feature.csv')\n\ntr = tr.merge(feature_ct, on='Patient')\n\nchunk = pd.read_csv(f\"{ROOT}/test.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path,resize_type='no'):\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 = abs(slices[-1].ImagePositionPatient[2] - slices[0].ImagePositionPatient[2])/(len(slices))\n    except:\n        try:\n            slice_thickness = abs(slices[-1].SliceLocation - slices[0].SliceLocation)/(len(slices))\n        except:\n            slice_thickness = slices[0].SliceThickness\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        if resize_type == 'resize':\n            s.PixelSpacing = s.PixelSpacing*(s.Rows/512)  \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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Test_Generate(Dataset):\n    def __init__(self,imgs_dicom):\n        self.imgs_dicom = imgs_dicom\n    def __getitem__(self,index):\n        metainf = self.imgs_dicom[index]\n        slice_img = metainf.pixel_array\n        slice_img = (slice_img-slice_img.min())/(slice_img.max()-slice_img.min())\n        slice_img = (slice_img*255).astype(np.uint8)\n            \n        if metainf.Rows!=512 or metainf.Columns!=512:\n            slice_img = cv2.resize(slice_img,(512,512))\n        \n        slice_img = slice_img[None,:,:]\n        slice_img = (slice_img/255).astype(np.float32)\n        return slice_img\n        \n    def __len__(self):\n        return len(self.imgs_dicom)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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(img,(512,512))\n    return img\n\ndef preprocess_img(img,resize_type):\n    if resize_type == 'resize':\n        img = [resize_image(im) for im in img]\n    if resize_type == 'crop':\n        img = [crop_image(im) for im in img]\n        \n    return np.array(img, dtype=np.int64)\n\n\ndef caculate_lung_volume(patient_scans,patient_masks):\n    \"\"\"\n    caculate volume of lung from mask\n    Parameters: list dicom scans,list patient CT Mask\n    Returns: volume cm³　(float)\n    \"\"\"\n    lung_volume = 0\n    for i in range(len(patient_masks)):\n        \n        pixel_spacing = patient_scans[i].PixelSpacing\n        slice_thickness = patient_scans[i].SliceThickness\n        lung_volume += np.count_nonzero(patient_masks[i])*pixel_spacing[0]*pixel_spacing[1]*slice_thickness\n        \n    return lung_volume*0.001\n\ndef caculate_histgram_statistical(patient_images,patient_masks,thresh = [-600,-250]):\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,Median=0,Skew=0,Kurthosis=0,HAA=0,midMean=0,\n                                      midSkew=0,midKurthosis=0,midMedian=0,midHAA=0)\n    num_slices = len(patient_images)\n    \n    #patient_images = patient_images[int(num_slices*0.1):int(num_slices*0.9)]\n    #patient_masks = patient_masks[int(num_slices*0.1):int(num_slices*0.9)]\n    patient_images = patient_masks*patient_images\n    patient_images_nonzero = patient_images[np.nonzero(patient_images)]\n    s_pixel = patient_images_nonzero.flatten()\n    haa_pixel = s_pixel[np.where((s_pixel>-1000)&(s_pixel<0))]\n    \n    mid_index = np.argsort(np.sum(patient_masks,axis=(1,2)))[-1]\n    \n    mid_image = patient_images[mid_index]\n    mid_images_nonzero = mid_image[np.nonzero(mid_image)]\n    mid_pixel = mid_images_nonzero.flatten()\n    midhaa_pixel = mid_pixel[np.where((mid_pixel>thresh[0])&(mid_pixel<thresh[1]))]\n    \n    \n    statistical_characteristic['Mean'] = np.mean(s_pixel)\n    statistical_characteristic['Median'] = np.median(s_pixel)\n    statistical_characteristic['Skew'] = skew(s_pixel)\n    statistical_characteristic['Kurthosis'] = kurtosis(s_pixel)\n    statistical_characteristic['HAA'] = len(haa_pixel)/len(s_pixel)\n    \n    statistical_characteristic['midMean'] = np.mean(mid_pixel)\n    statistical_characteristic['midMedian'] = np.median(mid_pixel)\n    statistical_characteristic['midSkew'] = skew(mid_pixel)\n    statistical_characteristic['midKurthosis'] = kurtosis(mid_pixel)\n    statistical_characteristic['midHAA'] = len(midhaa_pixel)/len(mid_pixel)\n    \n    for r in range(0,1000,100):\n        area_pixel = s_pixel[np.where((s_pixel>-r-100)&(s_pixel<-r))]\n        statistical_characteristic[f'pro_{r}'] = len(area_pixel)/len(s_pixel)\n    return statistical_characteristic","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"device =  torch.device('cuda:0')\ncheckpoint = '../input/2020osic/best_lung_Unet_densenet121_my.pth'\nModel_M = smp.Unet('densenet121', classes=1, in_channels=1,activation='sigmoid',encoder_weights=None).to(device)\nModel_M.load_state_dict(torch.load(checkpoint))\n#Model_M.eval()\n\ndef Unet_mask(model_m: nn.Module,input_data: DataLoader):\n    model_m.eval()\n    outs = []\n    for idx, sample in enumerate(test_loader):\n        image = sample\n        image = image.to(device)\n        with torch.no_grad():\n            out = model_m(image)\n        out = out.cpu().data.numpy()\n        out = np.where(out>0.5,1,0)\n        out = np.squeeze(out,axis=1)\n        outs.append(out)\n\n    outs = np.concatenate(outs)\n    return outs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ct_root_path = '../input/osic-pulmonary-fibrosis-progression/test/'\n\nlung_stat_pd = pd.DataFrame(columns=['Patient','Volume','Mean','Median','Skew','Kurthosis','HAA',\n                                     'midMean','midMedian','midSkew','midKurthosis','midHAA',\n                                    'pro_0','pro_100','pro_200','pro_300','pro_400','pro_500',\n                                     'pro_600','pro_700','pro_800','pro_900'])\nhaa_thresh = [-600,-250]\nfor i,p in enumerate(tqdm(pd.unique(chunk['Patient']))):\n    \n    lung_stat_pd.loc[i,'Patient'] = p    \n    patient_scans = load_scan(ct_root_path + p)\n    test_db = Test_Generate(patient_scans)\n    test_loader = DataLoader(test_db, batch_size=8, shuffle=False, num_workers=4)\n    masks = Unet_mask(Model_M,test_loader)\n    \n    patient_images = transform_to_hu(patient_scans)\n    if patient_images[0].shape!=(512,512):\n        patient_images = preprocess_img(patient_images,'resize')\n    \n    lung_stat_pd.loc[i,'Volume'] = caculate_lung_volume(patient_scans,masks)                           \n   \n    statistical_characteristic = caculate_histgram_statistical(patient_images,masks,haa_thresh)\n    \n    lung_stat_pd.loc[i,'Mean'] = statistical_characteristic['Mean']\n    lung_stat_pd.loc[i,'Median'] = statistical_characteristic['Median']\n    lung_stat_pd.loc[i,'Skew'] = statistical_characteristic['Skew']\n    lung_stat_pd.loc[i,'Kurthosis'] = statistical_characteristic['Kurthosis']\n    lung_stat_pd.loc[i,'HAA'] = statistical_characteristic['HAA']\n    \n    lung_stat_pd.loc[i,'midMean'] = statistical_characteristic['midMean']\n    lung_stat_pd.loc[i,'midMedian'] = statistical_characteristic['midMedian']\n    lung_stat_pd.loc[i,'midSkew'] = statistical_characteristic['midSkew']\n    lung_stat_pd.loc[i,'midKurthosis'] = statistical_characteristic['midKurthosis']\n    lung_stat_pd.loc[i,'midHAA'] = statistical_characteristic['midHAA']\n    \n    lung_stat_pd.loc[i,'pro_0'] = statistical_characteristic['pro_0']\n    lung_stat_pd.loc[i,'pro_100'] = statistical_characteristic['pro_100']\n    lung_stat_pd.loc[i,'pro_200'] = statistical_characteristic['pro_200']\n    lung_stat_pd.loc[i,'pro_300'] = statistical_characteristic['pro_300']\n    lung_stat_pd.loc[i,'pro_400'] = statistical_characteristic['pro_400']\n    lung_stat_pd.loc[i,'pro_500'] = statistical_characteristic['pro_500']\n    lung_stat_pd.loc[i,'pro_600'] = statistical_characteristic['pro_600']\n    lung_stat_pd.loc[i,'pro_700'] = statistical_characteristic['pro_700']\n    lung_stat_pd.loc[i,'pro_800'] = statistical_characteristic['pro_800']\n    lung_stat_pd.loc[i,'pro_900'] = statistical_characteristic['pro_900']\n\nlung_stat_pd.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"chunk = chunk.merge(lung_stat_pd)\n\nprint(\"add infos\")\nsub = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\ndata = tr.append([chunk, sub])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = (\n    data\n    .loc[data.Weeks == data.min_week][['Patient','FVC']]\n    .rename({'FVC': 'min_FVC'}, axis=1)\n    .groupby('Patient')\n    .first()\n    .reset_index()\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\ndel base\n\nbase = (\n    data\n    .loc[data.Weeks == data.min_week][['Patient','Percent']]\n    .rename({'Percent': 'min_Percent'}, axis=1)\n    .groupby('Patient')\n    .first()\n    .reset_index()\n)\ndata = data.merge(base, on='Patient', how='left')\ndel base\n\ndata.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\ndata['age'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\n#data['BASE'] = (data['min_FVC'] - data['min_FVC'].min() ) / ( data['min_FVC'].max() - data['min_FVC'].min())\ndata['week'] = data['base_week']# - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\ndata['percent'] = (data['min_Percent'] - data['min_Percent'].min() ) / ( data['min_Percent'].max() - data['min_Percent'].min())\n\ndata['volume'] = (data['Volume'] - tr['Volume'].min() ) / ( tr['Volume'].max() - tr['Volume'].min())\ndata['mean'] = (data['Mean'] - tr['Mean'].min()) / (tr['Mean'].max() - tr['Mean'].min())\ndata['skew'] = (data['Skew'] - tr['Skew'].min())/(tr['Skew'].max() - tr['Skew'].min())\ndata['median'] = (data['Median'] - tr['Median'].min()) / (tr['Median'].max() - tr['Median'].min())\ndata['kurthosis'] = (data['Kurthosis'] - tr['Kurthosis'].min())/(tr['Kurthosis'].max() - tr['Kurthosis'].min())\ndata['haa'] = (data['HAA'] - tr['HAA'].min())/(tr['HAA'].max() - tr['HAA'].min())\n\n\ndata['midmedian'] = (data['midMedian'] - tr['midMedian'].min() ) / ( tr['midMedian'].max() - tr['midMedian'].min())\ndata['midmean'] = (data['midMean'] - tr['midMean'].min()) / (tr['midMean'].max() - tr['midMean'].min())\ndata['midskew'] = (data['midSkew'] - tr['midSkew'].min())/(tr['midSkew'].max() - tr['midSkew'].min())\ndata['midkurthosis'] = (data['midKurthosis'] - tr['midKurthosis'].min())/(tr['midKurthosis'].max() - tr['midKurthosis'].min())\ndata['midhaa'] = (data['midHAA'] - tr['midHAA'].min())/(tr['midHAA'].max() - tr['midHAA'].min())\n\ndata['pro_0'] = (data['pro_0']-tr['pro_0'].min())/(tr['pro_0'].max() - tr['pro_0'].min())\ndata['pro_100'] = (data['pro_100']-tr['pro_100'].min())/(tr['pro_100'].max() - tr['pro_100'].min())\ndata['pro_200'] = (data['pro_200']-tr['pro_200'].min())/(tr['pro_200'].max() - tr['pro_200'].min())\ndata['pro_300'] = (data['pro_300']-tr['pro_300'].min())/(tr['pro_300'].max() - tr['pro_300'].min())\ndata['pro_400'] = (data['pro_400']-tr['pro_400'].min())/(tr['pro_400'].max() - tr['pro_400'].min())\ndata['pro_500'] = (data['pro_500']-tr['pro_500'].min())/(tr['pro_500'].max() - tr['pro_500'].min())\ndata['pro_600'] = (data['pro_600']-tr['pro_600'].min())/(tr['pro_600'].max() - tr['pro_600'].min())\ndata['pro_700'] = (data['pro_700']-tr['pro_700'].min())/(tr['pro_700'].max() - tr['pro_700'].min())\ndata['pro_800'] = (data['pro_800']-tr['pro_800'].min())/(tr['pro_800'].max() - tr['pro_800'].min())\ndata['pro_900'] = (data['pro_900']-tr['pro_900'].min())/(tr['pro_900'].max() - tr['pro_900'].min())\n\ndata['res_fvc'] = data['min_FVC']-data['FVC']\n\n\ndata.loc[:,\"Sex\"] = pd.factorize(data.Sex)[0]\ndata.loc[:,\"SmokingStatus\"] = pd.factorize(data.SmokingStatus)[0]\ndata['Sex'] = (data['Sex'] - data['Sex'].min() ) / ( data['Sex'].max() - data['Sex'].min() )\ndata['SmokingStatus'] = (data['SmokingStatus'] - data['SmokingStatus'].min() ) / ( data['SmokingStatus'].max() - data['SmokingStatus'].min())\n\"\"\"                                                                            \n\"\"\"\ndata['age'] = data['Age']\ndata['BASE'] = data['min_FVC']#-data['FVC']\ndata['res_FVC'] = data['min_FVC']-data['FVC']\ndata['week'] = data['base_week']# - data['base_week'].min()\ndata['percent'] = data['Percent']\n\ndata.head()\n\n#FE += ['age','percent','week','BASE','volume']\n#FE += ['age','percent','week','BASE','volume','mean','skew','kurthosis']\n#FE += ['Age','Percent','base_week','min_FVC','Volume','Mean','Skew','Kurthosis']\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndata['age'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\n#data['BASE'] = (data['min_FVC'] - data['min_FVC'].min() ) / ( data['min_FVC'].max() - data['min_FVC'].min())\ndata['week'] = data['base_week']# - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\ndata['percent'] = (data['min_Percent'] - data['min_Percent'].min() ) / ( data['min_Percent'].max() - data['min_Percent'].min())\n\ndata['volume'] = (data['Volume'] - data['Volume'].min() ) / ( data['Volume'].max() - data['Volume'].min())\ndata['mean'] = (data['Mean'] - data['Mean'].min()) / (data['Mean'].max() - data['Mean'].min())\ndata['skew'] = (data['Skew'] - data['Skew'].min())/(data['Skew'].max() - data['Skew'].min())\ndata['median'] = (data['Median'] - data['Median'].min()) / (data['Median'].max() - data['Median'].min())\ndata['kurthosis'] = (data['Kurthosis'] - data['Kurthosis'].min())/(data['Kurthosis'].max() - data['Kurthosis'].min())\ndata['haa'] = (data['HAA'] - data['HAA'].min())/(data['HAA'].max() - data['HAA'].min())\n\n\ndata['midmedian'] = (data['midMedian'] - data['midMedian'].min() ) / ( data['midMedian'].max() - data['midMedian'].min())\ndata['midmean'] = (data['midMean'] - data['midMean'].min()) / (data['midMean'].max() - data['midMean'].min())\ndata['midskew'] = (data['midSkew'] - data['midSkew'].min())/(data['midSkew'].max() - data['midSkew'].min())\ndata['midkurthosis'] = (data['midKurthosis'] - data['midKurthosis'].min())/(data['midKurthosis'].max() - data['midKurthosis'].min())\ndata['midhaa'] = (data['midHAA'] - data['midHAA'].min())/(data['midHAA'].max() - data['midHAA'].min())\n\ndata['pro_0'] = (data['pro_0']-data['pro_0'].min())/(data['pro_0'].max() - data['pro_0'].min())\ndata['pro_100'] = (data['pro_100']-data['pro_100'].min())/(data['pro_100'].max() - data['pro_100'].min())\ndata['pro_200'] = (data['pro_200']-data['pro_200'].min())/(data['pro_200'].max() - data['pro_200'].min())\ndata['pro_300'] = (data['pro_300']-data['pro_300'].min())/(data['pro_300'].max() - data['pro_300'].min())\ndata['pro_400'] = (data['pro_400']-data['pro_400'].min())/(data['pro_400'].max() - data['pro_400'].min())\ndata['pro_500'] = (data['pro_500']-data['pro_500'].min())/(data['pro_500'].max() - data['pro_500'].min())\ndata['pro_600'] = (data['pro_600']-data['pro_600'].min())/(data['pro_600'].max() - data['pro_600'].min())\ndata['pro_700'] = (data['pro_700']-data['pro_700'].min())/(data['pro_700'].max() - data['pro_700'].min())\ndata['pro_800'] = (data['pro_800']-data['pro_800'].min())/(data['pro_800'].max() - data['pro_800'].min())\ndata['pro_900'] = (data['pro_900']-data['pro_900'].min())/(data['pro_900'].max() - data['pro_900'].min())\n\ndata['res_fvc'] = data['min_FVC']-data['FVC']\n\n\ndata.loc[:,\"Sex\"] = pd.factorize(data.Sex)[0]\ndata.loc[:,\"SmokingStatus\"] = pd.factorize(data.SmokingStatus)[0]\ndata['Sex'] = (data['Sex'] - data['Sex'].min() ) / ( data['Sex'].max() - data['Sex'].min() )\ndata['SmokingStatus'] = (data['SmokingStatus'] - data['SmokingStatus'].min() ) / ( data['SmokingStatus'].max() - data['SmokingStatus'].min())\n                                                                                  \n\n\"\"\"\ndata['age'] = data['Age']\ndata['BASE'] = data['min_FVC']#-data['FVC']\ndata['res_FVC'] = data['min_FVC']-data['FVC']\ndata['week'] = data['base_week']# - data['base_week'].min()\ndata['percent'] = data['Percent']\n\ndata.head()\n\n#FE += ['age','percent','week','BASE','volume']\n#FE += ['age','percent','week','BASE','volume','mean','skew','kurthosis']\n#FE += ['Age','Percent','base_week','min_FVC','Volume','Mean','Skew','Kurthosis']\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Start get lung volume feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nFE = ['percent','SmokingStatus','Sex','age','week','volume','mean','skew','kurthosis','haa','median',\n     'midmean','midmedian','midskew','midkurthosis','midhaa',\n     'pro_0','pro_100','pro_200','pro_300','pro_400','pro_500',\n     'pro_600','pro_700','pro_800','pro_900']\n\n#FE = ['percent','age','week','volume','haa','pro_0','pro_700','pro_800','pro_900']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = data[data.WHERE=='train']\ntest = data[data.WHERE=='test']\nsub = data.loc[data.WHERE=='test']\ndel data\n#tr = tr.merge(feature_ct, on='Patient')\n#data = data.merge(lung_stat_pd, on='Patient', how='left')\ntr.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### BASELINE NN "},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch.nn.functional as F\n\nclass MishFunction(torch.autograd.Function):\n    @staticmethod\n    def forward(ctx, x):\n        ctx.save_for_backward(x)\n        return x * torch.tanh(F.softplus(x))   # x * tanh(ln(1 + exp(x)))\n\n    @staticmethod\n    def backward(ctx, grad_output):\n        x = ctx.saved_variables[0]\n        sigmoid = torch.sigmoid(x)\n        tanh_sp = torch.tanh(F.softplus(x)) \n        return grad_output * (tanh_sp + x * sigmoid * (1 - tanh_sp * tanh_sp))\n\nclass Mish(nn.Module):\n    def forward(self, x):\n        return MishFunction.apply(x)\n\ndef to_Mish(model):\n    for child_name, child in model.named_children():\n        if isinstance(child, nn.ReLU):\n            setattr(model, child_name, Mish())\n        else:\n            to_Mish(child)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = torch.tensor(70,dtype=torch.float),torch.tensor(1000,dtype=torch.float)\nC1, C2 = C1.to(device),C2.to(device)\n\ndef score_np(fvc_true, y_pred):\n    sigma = y_pred[:,2] - y_pred[:,0]\n    fvc_pred = y_pred[:,1]\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\n#=============================#\ndef score(y_true, y_pred):\n    y_true = y_true.to(torch.float)\n    y_pred = y_pred.to(torch.float)\n    \n    sigma = y_pred[:,2] - y_pred[:,0]\n    fvc_pred = y_pred[:,1]\n    #sigma_clip = sigma + C1\n    sigma_clip = torch.max(sigma, C1)\n    delta = torch.abs(y_true[:,0] - fvc_pred)\n    delta = torch.min(delta, C2)\n    sq2 = torch.sqrt(torch.tensor(2, dtype=torch.float))\n    metric = (delta / sigma_clip)*sq2 + torch.log(sigma_clip* sq2)\n    return torch.mean(metric)\n\n#============================#\ndef qloss(y_true, y_pred):\n    # Pinball loss for multiple quantiles\n    device = y_true.device\n    qs = [0.25, 0.50, 0.75]\n    q = torch.tensor(np.array([qs]), dtype=torch.float32)\n    q = q.to(device)\n    e = y_true - y_pred\n    v = torch.max(q*e, (q-1)*e)\n    return torch.mean(v)\n\n\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#=================\n\"\"\"\nclass make_model(nn.Module):\n    def __init__(self, in_ch, out_ch=3):\n        super(make_model, self).__init__()\n        self.fc1 = nn.Sequential(\n            nn.Linear(in_ch, 160),\n            Mish(),\n            #nn.BatchNorm1d(160),\n            nn.Dropout(0.3),\n        )\n        self.fc2 = nn.Sequential(\n            nn.Linear(160, 128),\n            Mish(),\n            #nn.BatchNorm1d(128),\n            nn.Dropout(0.25),\n        )\n        #self.fc3_p1 = nn.Linear(128, out_ch)\n        self.fc3_p1 = nn.Sequential(\n            nn.Linear(128, out_ch),\n            Mish()\n        )\n        \n        self.fc3_p2 = nn.Sequential(\n            nn.Linear(128, out_ch),\n            Mish()\n        )\n        #self.dropout = nn.Dropout(0.5)\n        \n    def forward(self, x):\n        x = self.fc1(x)\n        x = self.fc2(x)\n        x1 = self.fc3_p1(x)\n        x2 = self.fc3_p2(x)\n        x = x1 + torch.cumsum(x2,dim=1)\n        return x\n\"\"\"\n\nclass GaussianNoise(nn.Module):\n    \"\"\"Gaussian noise regularizer.\n\n    Args:\n        sigma (float, optional): relative standard deviation used to generate the\n            noise. Relative means that it will be multiplied by the magnitude of\n            the value your are adding the noise to. This means that sigma can be\n            the same regardless of the scale of the vector.\n        is_relative_detach (bool, optional): whether to detach the variable before\n            computing the scale of the noise. If `False` then the scale of the noise\n            won't be seen as a constant but something to optimize: this will bias the\n            network to generate vectors with smaller values.\n    \"\"\"\n    def __init__(self, sigma=0.1, is_relative_detach=True):\n        super().__init__()\n        self.sigma = sigma\n        self.is_relative_detach = is_relative_detach\n        self.register_buffer('noise', torch.tensor(0))\n\n    def forward(self, x):\n        if self.training and self.sigma != 0:\n            scale = self.sigma * x.detach() if self.is_relative_detach else self.sigma * x\n            sampled_noise = self.noise.expand(*x.size()).float().normal_() * scale\n            x = x + sampled_noise\n        return x\n    \nclass Model(nn.Module):\n    def __init__(self, num_features, num_targets, hidden_size=512):\n        super(Model, self).__init__()\n        self.batch_norm1 = nn.BatchNorm1d(num_features)\n        self.dropout1 = nn.Dropout(0.2)\n        self.dense1 = nn.utils.weight_norm(nn.Linear(num_features, hidden_size))\n        \n        self.batch_norm2 = nn.BatchNorm1d(hidden_size)\n        self.dropout2 = nn.Dropout(0.5)\n        self.dense2 = nn.utils.weight_norm(nn.Linear(hidden_size, hidden_size))\n        \n        self.batch_norm3 = nn.BatchNorm1d(hidden_size)\n        self.dropout3 = nn.Dropout(0.5)\n        self.dense3 = nn.utils.weight_norm(nn.Linear(hidden_size, num_targets))\n    \n    def forward(self, x):\n        \n        x = self.batch_norm1(x)\n        x = self.dropout1(x)\n        x = F.relu(self.dense1(x))\n        \n        x = self.batch_norm2(x)\n        x = self.dropout2(x)\n        x = F.relu(self.dense2(x))\n        \n        x = self.batch_norm3(x)\n        x = self.dropout3(x)\n        x = self.dense3(x)\n        \n        return x\n    \nclass make_big_model(nn.Module):\n    def __init__(self, in_ch, out_ch=3):\n        super(make_big_model, self).__init__()\n        self.gasus_nosie = GaussianNoise()\n        self.fc1 = nn.Sequential(\n            nn.Linear(in_ch, 256),\n            Mish()\n        )\n        self.fc2 = nn.Sequential(\n            nn.Linear(256, 256),\n            Mish(),\n            nn.Dropout(0.5)\n        )\n        self.fc3_p1 = nn.Sequential(\n            nn.Linear(256, 128),\n            Mish(),\n            nn.Dropout(0.5)\n        )\n    \n        self.fc3_p2 = nn.Linear(128, out_ch)\n        \n        \n        \n    def forward(self, x):\n        x = self.gasus_nosie(x)\n        x = self.fc1(x)\n        x = self.fc2(x)\n        x = self.fc3_p1(x)\n        #x = self.dropout(x)\n        x = self.fc3_p2(x)\n        return x\n    \nclass make_model(nn.Module):\n    def __init__(self, in_ch, out_ch=3):\n        super(make_model, self).__init__()\n        self.fc1 = nn.Sequential(\n            nn.Linear(in_ch, 160),\n            Mish()\n        )\n        self.fc2 = nn.Sequential(\n            nn.Linear(160, 128),\n            Mish()\n        )\n        self.fc3_p1 = nn.Sequential(\n            nn.Linear(128, 64),\n            Mish()\n        )\n        #self.dropout = nn.Dropout(0.2)\n        self.fc3_p2 = nn.Linear(128, out_ch)\n        \n        \n        \n    def forward(self, x):\n        x = self.fc1(x)\n        x = self.fc2(x)\n        #x = self.fc3_p1(x)\n        #x = self.dropout(x)\n        x = self.fc3_p2(x)\n        return x\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NFOLD = 5\nkf = KFold(n_splits=NFOLD)\npd_patient = pd.DataFrame({\"Patient\":tr[\"Patient\"].unique()})\n\nfor idx, (tr_idx, val_idx) in enumerate(kf.split(pd_patient)):\n    pd_patient.loc[val_idx,\"fold\"] = idx\npd_patient.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = tr.merge(pd_patient, on='Patient', how='left')\ntr.head()\n\n#tr['fold'] = -1\n#for i in range(len(pd_patient)):\n#    tr.fold[tr.Patient==pd_patient.loc[i,\"Patient\"]] = pd_patient.loc[i,\"fold\"]\n#tr.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nbatch = 128","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Data_Generate(Dataset):\n    def __init__(self,data,label=None):\n        self.data = data\n        self.label = label\n        \n    def __getitem__(self,index):\n        z_ = self.data[index]\n        if self.label is not None:\n            y_ = self.label[index]\n            y_ = y_[None,]\n            return z_,y_\n        else:\n            return z_\n         \n    def __len__(self):\n        return len(self.data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class EarlyStopping:\n    \"\"\"Early stops the training if validation loss doesn't improve after a given patience.\"\"\"\n    def __init__(self, patience=7, verbose=False, delta=0, path='checkpoint.pt'):\n        \"\"\"\n        Args:\n            patience (int): How long to wait after last time validation loss improved.\n                            Default: 7\n            verbose (bool): If True, prints a message for each validation loss improvement. \n                            Default: False\n            delta (float): Minimum change in the monitored quantity to qualify as an improvement.\n                            Default: 0\n            path (str): Path for the checkpoint to be saved to.\n                            Default: 'checkpoint.pt'\n        \"\"\"\n        self.patience = patience\n        self.verbose = verbose\n        self.counter = 0\n        self.best_score = None\n        self.early_stop = False\n        self.val_loss_min = np.Inf\n        self.delta = delta\n        self.path = path\n\n    def __call__(self, val_loss, model):\n\n        score = -val_loss\n\n        if self.best_score is None:\n            self.best_score = score\n            self.save_checkpoint(val_loss, model)\n        elif score < self.best_score + self.delta:\n            self.counter += 1\n            #print(f'EarlyStopping counter: {self.counter} out of {self.patience}')\n            if self.counter >= self.patience:\n                self.early_stop = True\n        else:\n            self.best_score = score\n            self.save_checkpoint(val_loss, model)\n            self.counter = 0\n\n    def save_checkpoint(self, val_loss, model):\n        '''Saves model when validation loss decrease.'''\n        if self.verbose:\n            print(f'Validation loss decreased ({self.val_loss_min:.6f} --> {val_loss:.6f}).  Saving model ...')\n        torch.save(model.state_dict(), self.path)\n        self.val_loss_min = val_loss","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FE","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n%%time\ncnt = 0\nEPOCHS = 1000\ncriterion = mloss(0.3)\n#for tr_idx, val_idx in kf.split(z):\nfor fold in range(NFOLD):\n    list_train_loss,list_val_loss,list_train_score,list_val_score = [],[],[],[]\n    val_out = []\n    print(f\"FOLD {fold+1}\")\n    #==================load data kfold==========================#\n    tr_z = tr[FE][tr.fold!=fold].values.astype(np.float32)\n    tr_y = tr.res_fvc[tr.fold!=fold].values.astype(np.float32)\n    val_z = tr[FE][tr.fold==fold].values.astype(np.float32)\n    val_y = tr.res_fvc[tr.fold==fold].values.astype(np.float32)\n    base_tr = tr['min_FVC'][tr.fold!=fold].values.astype(np.float32)\n    base_val = tr['min_FVC'][tr.fold!=fold].values.astype(np.float32)\n    train_db = Data_Generate(tr_z,tr_y)\n    train_loader = DataLoader(train_db, batch_size=batch, shuffle=True, num_workers=4)\n    val_db = Data_Generate(val_z,val_y)\n    val_loader = DataLoader(val_db, batch_size=batch, shuffle=False, num_workers=4)\n    \n    #==================prepare model==========================#\n    tr_num_batch = len(train_loader)\n    val_num_batch = len(val_loader)\n    net = make_big_model(len(FE),3).to(device)\n    optimizer = torch.optim.Adam(filter(lambda p: p.requires_grad, net.parameters()), lr=0.001 ,weight_decay=5e-4)\n    scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='min', factor=0.5, patience=20, min_lr=1e-8, verbose=False)\n    early_stopping = EarlyStopping(patience=80,path=f'Osic-NN-fold_{fold}.pth',verbose=False)\n\n    for epoch in tqdm(range(EPOCHS)):\n        train_loss,train_score,val_loss,val_score = 0,0,0,0\n        #==================train ==========================#\n        net.train()\n        for idx, sample in enumerate(train_loader):\n            feature, label = sample\n          \n            feature, label = feature.to(device), label.to(device)\n            out = net(feature)\n            loss = criterion(label, out)\n            score_ = score(label ,out)\n            optimizer.zero_grad()\n            loss.backward()\n            optimizer.step()\n            train_loss += loss.item()/tr_num_batch\n            train_score += score_.item()/tr_num_batch\n        list_train_loss.append(train_loss)\n        list_train_score.append(train_score)\n        #==================val ==========================#\n        net.eval()   \n        for idx, sample in enumerate(val_loader):\n            feature, label = sample\n            feature, label = feature.to(device), label.to(device)\n            with torch.no_grad():\n                out = net(feature)\n            #val_out.append(out.cpu().numpy())\n            loss = criterion(label, out)\n            score_ = score(label, out)\n            val_loss += loss.item()/val_num_batch\n            val_score += score_.item()/val_num_batch\n        list_val_loss.append(val_loss)\n        list_val_score.append(val_score)\n        early_stopping(val_score, net)\n        if early_stopping.early_stop:\n            print(\"Early stopping\")\n            break\n        scheduler.step(val_loss)\n             \n    print(f\"train loss: {min(list_train_loss)}  train score: {min(list_train_score)}\\n \\\n          val loss: {min(list_val_loss)} val score: {min(list_val_score)}\\n \\\n          final lr: {optimizer.param_groups[0]['lr']}\"\n         )\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score_np(fvc_true, fvc_pred, sigma):\n    sigma_clip = np.maximum(sigma, 70) \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\n\nscores,sigma,mean_FVC,true_y = [],[],[],[]\n\n#tr = new_tr\ntr['pred_FVC']=-1\ntr['Confidence']=-1\npreds = []\nfor fold in range(NFOLD):\n    pred = []\n    net = make_big_model(in_ch=len(FE)).to(device)\n    net.load_state_dict(torch.load(f\"Osic-NN-fold_{fold}.pth\"))\n    net.eval()\n    \n    tr_z = tr[FE][tr.fold==fold].values.astype(np.float32)\n    tr_y = tr.FVC[tr.fold==fold].values.astype(np.float32)\n    base_val = tr['min_FVC'][tr.fold==fold].values.astype(np.float32)\n    valid_db = Data_Generate(tr_z,tr_y)\n    valid_loader = DataLoader(valid_db, batch_size=batch, shuffle=False, num_workers=4)\n    \n    for idx, sample in enumerate(valid_loader):\n        feature, label = sample\n        feature, label = feature.to(device), label.to(device)\n        with torch.no_grad():\n            out = net(feature)\n        out = out.cpu().numpy() \n        pred.append(out)\n    pred = np.concatenate(pred)# + base_val[:,None]\n    preds.append(pred)\n    tr.loc[tr.fold==fold,'pred_FVC'] = base_val - pred[:,1]\n    tr.loc[tr.fold==fold,'pred_Confidence'] = abs(pred[:,2]-pred[:,0])\n    \n    \norder_tr = tr.groupby(tr['Patient']).apply(lambda x:x.sort_values(\"Weeks\",ascending=False))\np_id = pd.unique(tr.Patient)\n\nfor i in range(len(tr)):\n    sigma.append(tr.iloc[i].pred_Confidence)\n    mean_FVC.append(tr.iloc[i].pred_FVC)\n    true_y.append(tr.iloc[i].FVC)\n\nscore_np(np.array(true_y),np.array(mean_FVC),np.array(sigma))\n#np.array(scores).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_mean,last_mean = [],[]\nfor i in tr.Patient.unique():\n    sigma = tr[tr.Patient==i].pred_Confidence\n    mean_FVC = tr[tr.Patient==i].pred_FVC\n    true_y = tr[tr.Patient==i].FVC\n    all_mean.append(score_np(np.array(true_y),np.array(mean_FVC),np.array(sigma)))\n    last_mean.append(score_np(np.array(true_y[-3:]),np.array(mean_FVC[-3:]),np.array(sigma[-3:])))\n    \nprint(np.mean(all_mean),np.mean(last_mean))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr_1 = tr.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\nscores,sigma,mean_FVC,true_y = [],[],[],[]\nfor i in p_id:\n    last3_index = order_tr.loc[i,'Weeks'].index[:3]\n    sigma.append(order_tr.iloc[last3_index].pred_Confidence)\n    mean_FVC.append(order_tr.iloc[last3_index].pred_FVC)\n    true_y.append(order_tr.iloc[last3_index].FVC)\n    \n    scores.append(score_np(np.array(true_y),np.array(mean_FVC),np.array(sigma)))\n\nnp.array(scores).mean()\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv(f'{ROOT}/test.csv')\nfor i in test.Patient:\n    sub.loc[sub.Patient==i,'week'] = sub[sub.Patient==i].Weeks.values-test[test.Patient==i].Weeks.values\n    \nsub","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ze = (sub[FE].values).astype(np.float32)\nbase_eval = (sub['min_FVC'].values).astype(np.float32)\npe = np.zeros((ze.shape[0], 3))\n\ntest_db = Data_Generate(ze)\ntest_loader = DataLoader(test_db, batch_size=batch, shuffle=False, num_workers=4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"predict test...\")\npe = 0\nfor k in range(NFOLD):\n    if k == 2:\n        continue\n    pred = []\n    base_fvc = sub.FVC.values.astype(np.float32)\n    net = make_big_model(len(FE)).to(device)\n    net.load_state_dict(torch.load(f\"Osic-NN-fold_{k}.pth\"))\n    net.eval()\n    for idx, sample in enumerate(test_loader):\n        data = sample\n        data = data.to(device)\n        with torch.no_grad():\n            out = net(data)\n        out = out.cpu().numpy()\n        pred.append(out)\n    pred = base_fvc[:,None] - np.concatenate(pred)\n    pe += pred / 4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC1'] = pe[:, 1]\nsub['Confidence1'] = abs(pe[:, 2] - pe[:, 0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm = sub[['Patient_Week','FVC1','Confidence1']].copy()\nsubm.rename(columns={'FVC1':'FVC','Confidence1':'Confidence'},inplace=True) \nsubm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.ensemble import RandomForestRegressor,GradientBoostingRegressor\nfrom sklearn.multioutput import MultiOutputRegressor\nimport joblib","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\n\ntr = pd.read_csv(f\"{ROOT}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nfeature_ct = pd.read_csv('../input/2020osic/CT_21feature.csv')\n\ntr = tr.merge(feature_ct, on='Patient')\nchunk = pd.read_csv(f\"{ROOT}/test.csv\")\n\nprint(\"add infos\")\nsub = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")\n\ntr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\n\ndata = tr.append([chunk, sub])\ndata['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')\n\nbase = (\n    data\n    .loc[data.Weeks == data.min_week][['Patient','FVC']]\n    .rename({'FVC': 'min_FVC'}, axis=1)\n    .groupby('Patient')\n    .first()\n    .reset_index()\n)\n\ndata = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\ndel base\n\nbase = (\n    data\n    .loc[data.Weeks == data.min_week][['Patient','Percent']]\n    .rename({'Percent': 'min_Percent'}, axis=1)\n    .groupby('Patient')\n    .first()\n    .reset_index()\n)\ndata = data.merge(base, on='Patient', how='left')\ndel base\n\ndata.loc[:,\"Sex\"] = pd.factorize(data.Sex)[0]\ndata.loc[:,\"SmokingStatus\"] = pd.factorize(data.SmokingStatus)[0]\n\ntrain_df = data[data.WHERE=='train']\n\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_reg = pd.DataFrame(columns=['Patient', 'q2_Slope', 'q5_Slope','q8_Slope','slp'])\n\ndef model_slope(x):\n    b = x[0][0]\n    slope_1 = [(t[0]-b)/t[1] for t in x[1:]]\n    sort_slope = sorted(slope_1,reverse=True)\n    length = len(slope_1)\n    return np.mean(slope_1),np.mean(sort_slope[:length//2]),np.mean(sort_slope[length//2:])\n    \nfor idx,i in enumerate(pd.unique(train_df['Patient'])):\n    FVC_Weeks = train_df[train_df.Patient==i][['FVC','base_week']].values\n    slope_m,slope_f,slope_l = model_slope(FVC_Weeks)\n    \n \n    sub = train_df.loc[train_df.Patient == i, :] \n    fvc = sub.FVC.values\n    weeks = sub.Weeks.values\n    c = np.vstack([weeks, np.ones(len(weeks))]).T\n    a, b = np.linalg.lstsq(c, fvc)[0]\n\n    \n    train_reg.loc[idx] = [i,slope_m,slope_f,slope_l,a]\ntrain_reg.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"final_df = pd.DataFrame(columns=train_df.columns)\n\nfor idx,i in enumerate(pd.unique(train_df['Patient'])):\n    temp_df = train_df[train_df.Patient==i]\n    final_df = final_df.append(temp_df[temp_df.Weeks==temp_df.min_week].iloc[0,:])\n\ntrain_reg = train_reg.merge(final_df)\ntrain_reg = train_reg.drop(['min_week','base_week'],axis=1)\ntrain_reg.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_reg.loc[:,\"Sex\"] = pd.factorize(train_reg.Sex)[0]\ntrain_reg.loc[:,\"SmokingStatus\"] = pd.factorize(train_reg.SmokingStatus)[0]\n\ntrain_reg['SmokingStatus'] = (train_reg['SmokingStatus'] - train_reg['SmokingStatus'].min() ) / ( train_reg['SmokingStatus'].max() - train_reg['SmokingStatus'].min())\n\ntrain_reg.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NFOLD = 5\nkf = KFold(n_splits=NFOLD)\ntrain_reg['fold']=-1\nfor idx, (tr_idx, val_idx) in enumerate(kf.split(train_reg)):\n    train_reg.loc[val_idx,\"fold\"] = idx\ntrain_reg.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_reg = train_reg.merge(feature_ct)\ntrain_reg_bp = train_reg.copy()\ntrain_reg.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_reg['age'] = (train_reg['Age'] - train_reg['Age'].min() ) / ( train_reg['Age'].max() - train_reg['Age'].min() )\ntrain_reg['BASE'] = (train_reg['min_FVC'] - train_reg['min_FVC'].min() ) / ( train_reg['min_FVC'].max() - train_reg['min_FVC'].min())\ntrain_reg['percent'] = (train_reg['min_Percent'] - train_reg['min_Percent'].min() ) / ( train_reg['min_Percent'].max() - train_reg['min_Percent'].min())\ntrain_reg['volume'] = (train_reg['Volume'] - train_reg['Volume'].min() ) / ( train_reg['Volume'].max() - train_reg['Volume'].min())\ntrain_reg['mean'] = (train_reg['Mean'] - train_reg['Mean'].min()) / (train_reg['Mean'].max() - train_reg['Mean'].min())\ntrain_reg['skew'] = (train_reg['Skew'] - train_reg['Skew'].min())/(train_reg['Skew'].max() - train_reg['Skew'].min())\ntrain_reg['median'] = (train_reg['Median'] - train_reg['Median'].min()) / (train_reg['Median'].max() - train_reg['Median'].min())\ntrain_reg['kurthosis'] = (train_reg['Kurthosis'] - train_reg['Kurthosis'].min())/(train_reg['Kurthosis'].max() - train_reg['Kurthosis'].min())\ntrain_reg['haa'] = (train_reg['HAA'] - train_reg['HAA'].min())/(train_reg['HAA'].max() - train_reg['HAA'].min())\n\n\ntrain_reg['midmedian'] = (train_reg['midMedian'] - train_reg['midMedian'].min() ) / ( train_reg['midMedian'].max() - train_reg['midMedian'].min())\ntrain_reg['midmean'] = (train_reg['midMean'] - train_reg['midMean'].min()) / (train_reg['midMean'].max() - train_reg['midMean'].min())\ntrain_reg['midskew'] = (train_reg['midSkew'] - train_reg['midSkew'].min())/(train_reg['midSkew'].max() - train_reg['midSkew'].min())\ntrain_reg['midkurthosis'] = (train_reg['midKurthosis'] - train_reg['midKurthosis'].min())/(train_reg['midKurthosis'].max() - train_reg['midKurthosis'].min())\ntrain_reg['midhaa'] = (train_reg['midHAA'] - train_reg['midHAA'].min())/(train_reg['midHAA'].max() - train_reg['midHAA'].min())\n\ntrain_reg['pro_0'] = (train_reg['pro_0']-train_reg['pro_0'].min())/(train_reg['pro_0'].max() - train_reg['pro_0'].min())\ntrain_reg['pro_100'] = (train_reg['pro_100']-train_reg['pro_100'].min())/(train_reg['pro_100'].max() - train_reg['pro_100'].min())\ntrain_reg['pro_200'] = (train_reg['pro_200']-train_reg['pro_200'].min())/(train_reg['pro_200'].max() - train_reg['pro_200'].min())\ntrain_reg['pro_300'] = (train_reg['pro_300']-train_reg['pro_300'].min())/(train_reg['pro_300'].max() - train_reg['pro_300'].min())\ntrain_reg['pro_400'] = (train_reg['pro_400']-train_reg['pro_400'].min())/(train_reg['pro_400'].max() - train_reg['pro_400'].min())\ntrain_reg['pro_500'] = (train_reg['pro_500']-train_reg['pro_500'].min())/(train_reg['pro_500'].max() - train_reg['pro_500'].min())\ntrain_reg['pro_600'] = (train_reg['pro_600']-train_reg['pro_600'].min())/(train_reg['pro_600'].max() - train_reg['pro_600'].min())\ntrain_reg['pro_700'] = (train_reg['pro_700']-train_reg['pro_700'].min())/(train_reg['pro_700'].max() - train_reg['pro_700'].min())\ntrain_reg['pro_800'] = (train_reg['pro_800']-train_reg['pro_800'].min())/(train_reg['pro_800'].max() - train_reg['pro_800'].min())\ntrain_reg['pro_900'] = (train_reg['pro_900']-train_reg['pro_900'].min())/(train_reg['pro_900'].max() - train_reg['pro_900'].min())\n\ntrain_reg['Sex'] = (train_reg['Sex'] - train_reg['Sex'].min() ) / ( train_reg['Sex'].max() - train_reg['Sex'].min() )\ntrain_reg['SmokingStatus'] = (train_reg['SmokingStatus'] - train_reg['SmokingStatus'].min() ) / ( train_reg['SmokingStatus'].max() - train_reg['SmokingStatus'].min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n######################xgboost##############################\nfrom xgboost import XGBRegressor\nimport xgboost as xgb\nfrom sklearn.metrics import mean_absolute_error\n\nFE = ['SmokingStatus','Sex','age','BASE','volume','mean','skew','kurthosis','haa','median',\n     'midmean','midmedian','midskew','midkurthosis','midhaa',\n     'pro_0','pro_100','pro_200','pro_300','pro_400','pro_500',\n     'pro_600','pro_700','pro_800','pro_900']\n\npred_xgboost = []\ntrain_reg['pre_s_xgb'] = -1\nm,last3 = [],[]\nparas={\n    'booster':'gbtree',\n \n    'objective':'reg:squarederror',\n    'gamma':0.05,#树的叶子节点下一个区分的最小损失，越大算法模型越保守\n    'lambda':8,#L2正则项权重\n    'subsample':0.8,#采样训练数据，设置为0.5\n    'colsample_bytree':0.7,#构建树时的采样比率\n    'min_child_weight':1,#节点的最少特征数\n    'eta':0.05,#类似学习率\n    'seed':0,\n    'nthread':4,#cpu线程数\n    'eval_metric':'mae'\n}\n\nfor n in range(NFOLD):\n        \n    print(f\"FOLD {n+1}\")\n    #==================load data kfold==========================#\n    tr_x,val_x = train_reg[train_reg.fold!=n][FE].values,train_reg[train_reg.fold==n][FE].values\n    tr_y,val_y = train_reg[train_reg.fold!=n]['slp'].values,train_reg[train_reg.fold==n]['slp'].values\n\n \n    my_model = XGBRegressor(**paras,n_estimators=1000)\n    dtrain = xgb.DMatrix(tr_x, label=tr_y)\n    dval = xgb.DMatrix(val_x, label=val_y)\n    evallist = [(dval, 'eval'), (dtrain, 'train')]\n    num_round = 10\n    my_model.fit(tr_x, tr_y, early_stopping_rounds=10, \n             eval_set=[(val_x, val_y)], verbose=False)\n    y_pred = my_model.predict(val_x)\n    pred_xgboost.append(y_pred)\n    print('The rmse of prediction is:', mean_absolute_error(val_y, y_pred))\n    index = train_reg[train_reg.fold==n].index\n    \n    train_reg.loc[index,['pre_s_xgb']] = my_model.predict(val_x)\n    for p in train_reg[train_reg.fold==n].Patient:\n        \n        percent_true = train_df.Percent.values[train_df.Patient == p]\n        fvc_true = train_df.FVC.values[train_df.Patient == p]\n        weeks_true = train_df.Weeks.values[train_df.Patient == p]\n        a = train_reg[train_reg.Patient==p].pre_s_xgb.values\n        \n        fvc = a * (weeks_true - weeks_true[0]) + fvc_true[0]\n        percent = percent_true[0] - a * abs(weeks_true - weeks_true[0])\n        train_df.loc[train_df.Patient==p,'img_fvc_xgb'] = fvc\n        train_df.loc[train_df.Patient==p,'img_con_xgb'] = percent\n        m.append(score_np(fvc_true, fvc,percent))\n        last3.append(score_np(fvc_true[-3:], fvc[-3:],percent[-3:]))\n    print(f'all_mean:{(np.mean(m))},last3_mean:{np.mean(last3)}')\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FE = ['percent','SmokingStatus','Sex','age','BASE','volume','mean','skew','kurthosis','haa','median',\n     'midmean','midmedian','midskew','midkurthosis','midhaa',\n     'pro_0','pro_100','pro_200','pro_300','pro_400','pro_500',\n     'pro_600','pro_700','pro_800','pro_900']\n\ntrain_reg['pre_s'] = -1\nm,last3 = [],[]\nfor n in range(NFOLD):\n    tr_x,val_x = train_reg[train_reg.fold!=n][FE].values,train_reg[train_reg.fold==n][FE].values\n    tr_y,val_y = train_reg[train_reg.fold!=n]['slp'].values,train_reg[train_reg.fold==n]['slp'].values\n\n \n    randomforest = RandomForestRegressor(n_estimators = 500, oob_score = True, n_jobs = -1,\n                                    max_features = \"sqrt\", min_samples_leaf = 5, max_samples = 0.8)\n       \n    randomforest.fit(tr_x, tr_y)\n\n    joblib.dump(randomforest, f'./all_rf_{n}.pkl')\n    \n    index = train_reg[train_reg.fold==n].index\n    \n    train_reg.loc[index,['pre_s']] = randomforest.predict(val_x)\n    \n    s = np.mean(abs(train_reg.loc[index,['pre_s']].values-train_reg.loc[index,['slp']].values))\n    print(f\"{n}th model's out-of-bag score : \", s)\n    \n    for p in train_reg[train_reg.fold==n].Patient:\n        \n        percent_true = train_df.Percent.values[train_df.Patient == p]\n        fvc_true = train_df.FVC.values[train_df.Patient == p]\n        weeks_true = train_df.Weeks.values[train_df.Patient == p]\n        a = train_reg[train_reg.Patient==p].pre_s.values\n        \n        fvc = a * (weeks_true - weeks_true[0]) + fvc_true[0]\n        percent = percent_true[0] - a * abs(weeks_true - weeks_true[0])\n        train_df.loc[train_df.Patient==p,'img_fvc'] = fvc\n        train_df.loc[train_df.Patient==p,'img_con'] = percent\n        m.append(score_np(fvc_true, fvc,percent))\n        last3.append(score_np(fvc_true[-3:], fvc[-3:],percent[-3:]))\n    print(f'all_mean:{(np.mean(m))},last3_mean:{np.mean(last3)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv(f\"{ROOT}/test.csv\")\ntest_df = test.merge(lung_stat_pd,)\ntest_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_reg = train_reg_bp\ntest_df['age'] = (test_df['Age'] - train_reg['Age'].min() ) / ( train_reg['Age'].max() - train_reg['Age'].min() )\ntest_df['BASE'] = (test_df['FVC'] - train_reg['min_FVC'].min() ) / ( train_reg['min_FVC'].max() - train_reg['min_FVC'].min())\ntest_df['percent'] = (test_df['Percent'] - train_reg['min_Percent'].min() ) / ( train_reg['min_Percent'].max() - train_reg['min_Percent'].min())\ntest_df['volume'] = (test_df['Volume'] - train_reg['Volume'].min() ) / ( train_reg['Volume'].max() - train_reg['Volume'].min())\ntest_df['mean'] = (test_df['Mean'] - train_reg['Mean'].min()) / (train_reg['Mean'].max() - train_reg['Mean'].min())\ntest_df['skew'] = (test_df['Skew'] - train_reg['Skew'].min())/(train_reg['Skew'].max() - train_reg['Skew'].min())\ntest_df['median'] = (test_df['Median'] - train_reg['Median'].min()) / (train_reg['Median'].max() - train_reg['Median'].min())\ntest_df['kurthosis'] = (test_df['Kurthosis'] - train_reg['Kurthosis'].min())/(train_reg['Kurthosis'].max() - train_reg['Kurthosis'].min())\ntest_df['haa'] = (test_df['HAA'] - train_reg['HAA'].min())/(train_reg['HAA'].max() - train_reg['HAA'].min())\n\n\ntest_df['midmedian'] = (test_df['midMedian'] - train_reg['midMedian'].min() ) / ( train_reg['midMedian'].max() - train_reg['midMedian'].min())\ntest_df['midmean'] = (test_df['midMean'] - train_reg['midMean'].min()) / (train_reg['midMean'].max() - train_reg['midMean'].min())\ntest_df['midskew'] = (test_df['midSkew'] - train_reg['midSkew'].min())/(train_reg['midSkew'].max() - train_reg['midSkew'].min())\ntest_df['midkurthosis'] = (test_df['midKurthosis'] - train_reg['midKurthosis'].min())/(train_reg['midKurthosis'].max() - data['midKurthosis'].min())\ntest_df['midhaa'] = (test_df['midHAA'] - train_reg['midHAA'].min())/(train_reg['midHAA'].max() - train_reg['midHAA'].min())\n\ntest_df['pro_0'] = (test_df['pro_0']-train_reg['pro_0'].min())/(train_reg['pro_0'].max() - train_reg['pro_0'].min())\ntest_df['pro_100'] = (test_df['pro_100']-train_reg['pro_100'].min())/(train_reg['pro_100'].max() - train_reg['pro_100'].min())\ntest_df['pro_200'] = (test_df['pro_200']-train_reg['pro_200'].min())/(train_reg['pro_200'].max() - train_reg['pro_200'].min())\ntest_df['pro_300'] = (test_df['pro_300']-train_reg['pro_300'].min())/(train_reg['pro_300'].max() - train_reg['pro_300'].min())\ntest_df['pro_400'] = (test_df['pro_400']-train_reg['pro_400'].min())/(train_reg['pro_400'].max() - train_reg['pro_400'].min())\ntest_df['pro_500'] = (test_df['pro_500']-train_reg['pro_500'].min())/(train_reg['pro_500'].max() - train_reg['pro_500'].min())\ntest_df['pro_600'] = (test_df['pro_600']-train_reg['pro_600'].min())/(train_reg['pro_600'].max() - train_reg['pro_600'].min())\ntest_df['pro_700'] = (test_df['pro_700']-train_reg['pro_700'].min())/(train_reg['pro_700'].max() - train_reg['pro_700'].min())\ntest_df['pro_800'] = (test_df['pro_800']-train_reg['pro_800'].min())/(train_reg['pro_800'].max() - train_reg['pro_800'].min())\ntest_df['pro_900'] = (test_df['pro_900']-train_reg['pro_900'].min())/(train_reg['pro_900'].max() - train_reg['pro_900'].min())\n\n\ntest_df.loc[:,\"Sex\"] = pd.factorize(test_df.Sex)[0]\ntest_df.loc[:,\"SmokingStatus\"] = pd.factorize(test_df.SmokingStatus)[0]\ntest_df['Sex'] = (test_df['Sex'] - train_reg['Sex'].min() ) / ( train_reg['Sex'].max() - train_reg['Sex'].min() )\ntest_df['SmokingStatus'] = (test_df['SmokingStatus'] - train_reg['SmokingStatus'].min() ) / ( train_reg['SmokingStatus'].max() - train_reg['SmokingStatus'].min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")\nsub.head()\n\nFE = ['percent','SmokingStatus','Sex','age','BASE','volume','mean','skew','kurthosis','haa','median',\n     'midmean','midmedian','midskew','midkurthosis','midhaa',\n     'pro_0','pro_100','pro_200','pro_300','pro_400','pro_500',\n     'pro_600','pro_700','pro_800','pro_900']\n\nY = ['slp']\npred_y = []\n\nfor n in range(NFOLD):\n    val_x = test_df[FE]\n\n    randomforest = joblib.load( f'./all_rf_{n}.pkl') \n    \n    pred_y.append(randomforest.predict(val_x))\n\npred_y = np.array(pred_y).mean(0)\ntest_df.loc[:,'preslp'] = pred_y\n\nfor p in test_df.Patient:\n    a = test_df[test_df.Patient==p].preslp.values\n    weeks_base = test_df[test_df.Patient==p].Weeks.values\n    weeks_true = sub[sub.Patient==p].Weeks.values\n    percent_true = sub[sub.Patient==p].Percent.values\n    fvc_true = sub[sub.Patient==p].FVC.values\n    \n    fvc = a * (weeks_true - weeks_base[0]) + fvc_true[0]\n    percent = weeks_base[0] - a * abs(weeks_true - weeks_base[0])\n       \n    sub.loc[sub.Patient==p,'FVC'] = fvc\n    sub.loc[sub.Patient==p,'Confidence'] = percent\n\nsub.head()\n\nsubs = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"l = 0.\nm,last3 = [],[]\nfor i in train_df.Patient:\n    true_fvc = tr_1[tr_1.Patient==i].FVC\n    pred_fvc = (train_df[train_df.Patient==i].img_fvc)*l + (1-l)*tr_1[tr_1.Patient==i].pred_FVC\n    pred_sig = (train_df[train_df.Patient==i].img_con)*l + (1-l)*tr_1[tr_1.Patient==i].pred_Confidence\n    m.append(score_np(true_fvc, pred_fvc,pred_sig))\n    last3.append(score_np(true_fvc[-3:], pred_fvc[-3:],pred_sig[-3:]))\nprint(f'all_mean:{(np.mean(m))},last3_mean:{np.mean(last3)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}