{"cells":[{"metadata":{},"cell_type":"markdown","source":"# XY's solution\n\n## Features\n1. Basic features like age,sex,..\n2. Features generated by using lungmask (https://github.com/JoHof/lungmask/) as disclosed in https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/164883\n3. Segmented images using method as described in https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed\n\n## Models to predict FVC\n1. Custom Neural Network to predict slant of FVC decay\n2. 3 LGBMs to predict slant of FVC decay ( differs in calculation of slant )\n3. LGBM to predict FVC directly\n\n## Procedure to predict Confidence\n1. OOF prediction FVC is calculated.\n2. OOF prediction error is calculated.\n3. Calculate optimal confidence value to get best score when get that prediction error.\n4. Train fitter to predict that optimal confidnece value form features.\n5. Use that confience fitter to predict test confidence from test features.\n\n## Custom Neural Network\n- Inputs are features, original images, and segmented images.\n- Sequence of images are encoded using part of u-net(R231), which is the bottom of u-net conv layer having the lowest resolution output.\n- Encoded values are feed to LSTM\n- Final output if calculated by using table features and LSTM outputs\n\n## How this notebook works\n- Precalculated features are loaded (Datasets lungmasks5, classicsegmentedjpg )\n- Pretrained custom neural network's weights are loaded (Datasets osicnn0927segonlyseed0, 1, 2, 3, 4) ; pls don't care about the name \"segonly\" which actualy uses both orignal images and segmented images.\n- The other predictors are trained to calculate parameters\n- To generate submission file, first calculate features, then prediction by using above models.\n\n## Two submission file\n- One submission file is generated by above models, another is generated by ensemble of above models by changing 5 different fold splits.\n\n\n---\n\n"},{"metadata":{},"cell_type":"markdown","source":"## Config"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\n\nCNF={\n    'lstm_mem':64,\n    'mixup':True,\n    'lung_area_th':0.2, # maxの何割以下を落とすか\n    'calcAmethod':'simpleLinearRegress',#'huber'# last3 # last3からtrueAを計算する\n    'jpgPath':\"../input/classicsegmentedjpg/classic_segmented/classic_segmented\",\n    'isTrain':False,\n    'verbose':True,\n    'DBG':True,\n    'num_sequence':12,\n    'LUNG_MODEL':'R231',\n    'batch_size':3,  \n    'mskDfPath':'../input/lungmasks5/LTRC_feat_df.csv',\n    'trnPath':'../input/osic-pulmonary-fibrosis-progression/train',\n    'mskjpgSizePklPath':'../input/lungmasks/cropped_orig_size_dict.pkl',\n    'useLungVolume':True,\n    'maxEpoch':20,\n    'EvalOnlyLast3':True,\n    'aug':False, # Falseでも一応アスペクト対応だけはする\n    'kl_weight':0.01,\n    'EMA':True,\n    \n    \n    # -------------------------------------------\n    \n    'prefixes':[[f\"osicnn0927segonlyseed{_seed}\",f\"A0LGBMv1seed{_seed}\",f\"ARLGBMv1seed{_seed}\",f\"LGBMv1seed{_seed}\",f\"weekLGBMv1seed{_seed}\"] for _seed in range(5)],\n    'ensembleWeight':[[3,1,2,1,2] for _ in range(5)],\n    \n    # -------------------------------------------\n    \n    \n    'n_tta':1,\n    'useLGBM':True,\n    'quality':75, # jpg quality\n \n}\n\nCNF['ensembleWeight']=[ np.array(_ws)/np.sum(_ws) for _ws in CNF['ensembleWeight'] ]\n\nassert len(CNF['prefixes'])==len(CNF['ensembleWeight'])\n\nLUNG_CNF={\n    'R231':{'pth':'../input/pthunetlungmask/unet_r231-d5d2fc3d.pth','nc':3},\n    'COVID19':{'pth':'../input/pthunetlungmask/unet_r231covid-0de78a7e.pth','nc':3},\n    'LTRC':{'pth':'../input/pthunetlungmask/unet_ltrclobes-3a07043d.pth','nc':6}\n}\n\n!cp {CNF['jpgPath']+'/../../createClassicSegmented.py'} .\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# %%capture\n# !conda install -c conda-forge -y gdcm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%capture\n%cd ../input/python3gdcm\n!dpkg -i build_1-1_amd64.deb\n!apt-get install -f\n!cp /usr/local/lib/gdcm.py /opt/conda/lib/python3.7/site-packages/.\n!cp /usr/local/lib/gdcmswig.py /opt/conda/lib/python3.7/site-packages/.\n!cp /usr/local/lib/_gdcmswig.so /opt/conda/lib/python3.7/site-packages/.\n!cp /usr/local/lib/libgdcm* /opt/conda/lib/python3.7/site-packages/.\n!ldconfig","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%cd -","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# !pip install pydicom==2.0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport torch\nfrom torch.utils.data import Dataset, DataLoader\nimport torch.optim as optim\nimport os,sys,random,gc\nimport seaborn as sns\n\nimport cv2\nfrom matplotlib import pyplot as plt\nimport matplotlib.cm as cm\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_modality_lut\nimport pandas as pd\nfrom osicutil import UNet\nfrom scipy.stats import median_absolute_deviation as mad\nimport matplotlib.pyplot as plt\n\nfrom sklearn.model_selection import KFold,GroupKFold\n\nfrom osicutil import ModelEMA\nfrom osicutil import LinearARD,ELBOLoss,get_ard_reg,_get_dropped_params_cnt,_get_params_cnt,get_dropped_params_ratio # ARD\nfrom osicutil import calcScore\nfrom osicutil import _save,_load\n\ndef seed_everything(seed: int):\n    random.seed(seed);os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    np.random.seed(seed);torch.manual_seed(seed);torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n\nclass HorizontalDisplay:\n    def __init__(self, *args):\n        self.args = args\n    def _repr_html_(self):\n        template = '<div style=\"float: left; padding: 10px;\">{0}</div>'\n        return \"\\n\".join(template.format(arg._repr_html_()) for arg in self.args)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## create lung mask features"},{"metadata":{"trusted":true},"cell_type":"code","source":"def ceateVerticalVolumeFeatDf(path,prefixXorY):\n    \n    print(f\"{path}\")\n    ret_df = pd.DataFrame()\n    tmp_dict = dict()\n    \n    _eps=1e-6\n    \n    for fname in os.listdir(path):\n        \n        if not fname.endswith('.pkl'):\n            continue\n        \n        pid = fname[:-4]\n        \n        volume = _load(f\"{path}/{fname}\")\n        \n        tmp_dict = {'Patient':pid}\n        \n        \n        # mask6 & masked image6 & lung area1 total 13 class\n        \n        # 最後は肺領域 range(7) \n        for cls_id in [1,2]: # 3,4,5は意味ないっぽい\n            \n            volume = volume[:,50:-50]\n            # ltrc mask\n            tmp_dict[f\"vert{prefixXorY}_msk_cls{cls_id}_mad\"]  = mad(volume[cls_id],axis=-1)\n            \n            # ltrc mask*image\n            tmp_dict[f\"vert{prefixXorY}_mskdimage_cls{cls_id}_mad\"]  = mad(volume[cls_id+6],axis=-1)\n            \n            # ratio\n            tmp_dict[f\"vert{prefixXorY}_ratio_cls{cls_id}_mad\"]  = mad(volume[cls_id+6]/(volume[cls_id]+_eps),axis=-1)\n            \n        \n        ret_df = pd.concat([ret_df,pd.DataFrame(tmp_dict,index=[-1])]) # 間違いないようにindexに-1を設定\n    \n    return ret_df,ret_df.columns.drop('Patient').tolist()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Confidence fitter"},{"metadata":{"trusted":true},"cell_type":"code","source":"class ConfidenceFitterOptimizedLinear:\n    \n    def __init__(self,feat,allFVC,allPredFVC):\n        \n        from scipy.optimize import minimize\n        from sklearn.linear_model import LinearRegression\n        \n        self.reg = LinearRegression().fit(feat,abs(allFVC-allPredFVC))\n        \n        init_param = np.concatenate([self.reg.coef_, [self.reg.intercept_] ])\n        \n        def confidence_loss_func(param):\n            estimated_simga = feat@param[:-1]+param[-1]\n            return calcScore(allFVC, allPredFVC, estimated_simga)\n        \n        self.res = minimize(confidence_loss_func, init_param, method='Nelder-Mead', tol=1e-9)\n            \n    def predict(self,feat):\n        return feat@self.res.x[:-1]+self.res.x[-1]\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Add Feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.linear_model import HuberRegressor\nimport glob\n\n\ndef addFeatureDf(df,mskDfPath,dicomPath,is_test=False): \n    \n    feat_cols=['Age','Female', 'Male', 'Currently smokes', 'Ex-smoker', 'Never smoked','cat_coef','FVC0','Week0is0']\n    \n    df['Week0'] = df['Weeks'] # ----------------- append Week0,FVC0,rWeeks(week-week0)\n    df['Week0'] = df.groupby('Patient')['Week0'].transform('min')\n    df['Week0is0'] = (df['Week0']==0).astype('int')\n    \n    base = df.loc[df.Weeks == df.Week0]\n    base = base[['Patient', 'FVC','Percent']].copy()\n    base.columns = ['Patient', 'FVC0','Percent']\n    base['nb'] = 1\n    base['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\n    base['cat_coef'] = base['FVC0'].values/base['Percent'].values\n    base = base[base.nb == 1]\n    base.drop(['nb','Percent'], axis=1, inplace=True)\n    df = df.merge(base, on='Patient', how='left')\n\n    df['rWeeks'] = df['Weeks'] - df['Week0'] \n    \n    df = pd.concat([df, pd.get_dummies(df.Sex), pd.get_dummies(df.SmokingStatus)], axis=1).drop(columns=['Sex','SmokingStatus']) #  one-hot enc\n    for col in ['Female', 'Male', 'Currently smokes', 'Ex-smoker', 'Never smoked']: # if there is no such values in df, create and set 0\n        df[col] = 0 if col not in df.columns else df[col]\n        df[col] = df[col].astype('int')\n    \n    df.Age = ( df.Age - 30 ) / 30 \n    \n    if CNF['useLungVolume']: # -------------------- Create Lung Volume Features\n        \n        _df = pd.read_csv(mskDfPath)\n        if not is_test:\n            _df = _df.drop(index=_df.loc[(_df.Patient=='ID00078637202199415319443')&(_df.file_id>509)].index).reset_index(drop=True) # remove duplicated sequence \n        \n        # drop non lung image\n        _df['lung_area'] = _df['lung_area']*_df['dx']*_df['dy']\n        _lung_area_df = _df[['Patient','file_id','lung_area']]\n        #_df = _df.drop(index=_df.loc[(_df.lung_area<400)].index).reset_index(drop=True) \n        \n        for p in _lung_area_df.Patient.unique(): # ----------------------------------------------------  Drop small area\n            _th = _lung_area_df[_lung_area_df.Patient==p].lung_area.max()*CNF['lung_area_th']\n            _df = _df.drop(index=_df.loc[(_df.Patient==p)&(_df.lung_area<_th)].index).reset_index(drop=True) # new th for new mask\n        \n        \n        cls_ids = [ int(col[-1]) for col in _df.columns.tolist() if col.startswith('area_cls')]\n        \n        for cls_id in cls_ids:\n            _df[f\"area_cls{cls_id}\"]*=_df['dx']*_df['dy']\n            _df[f\"mskd_img_sum_cls{cls_id}\"]*=_df['dx']*_df['dy']\n            _df[f\"mskd_ratio_cls{cls_id}\"]=_df[f\"mskd_img_sum_cls{cls_id}\"]/_df[f\"area_cls{cls_id}\"]\n\n        mask_feat_cols = [ f\"area_cls{cls_id}\" for cls_id in cls_ids ]+[ f\"mskd_img_sum_cls{cls_id}\" for cls_id in cls_ids ]+[ f\"mskd_ratio_cls{cls_id}\" for cls_id in cls_ids ] \n        \n        mask_feat = _df.groupby('Patient')[mask_feat_cols].agg(['median','mean',mad]).reset_index()\n        \n        mask_feat.columns = [\"_\".join(pair) if '' not in pair else pair[0] for pair in mask_feat.columns ] #マルチインデックス解消\n        \n        df=df.merge(mask_feat,how='left',on='Patient')\n        feat_cols+=mask_feat.columns[1:-1].tolist() # 最後にはr231_lung_area\n        \n        # dicom info \n        aux_df = _df[['Patient','file_id']].groupby('Patient')['file_id'].apply(max).reset_index().rename(columns={'file_id':'max_file_id'})\n        df=df.merge(aux_df ,how='left',on='Patient')\n        feat_cols+=['max_file_id']\n        \n        \n    if not is_test: # ------------- Create Reg. Coef. 'A' ; FVC = A*(week-week0)+FVC0\n    \n        df['trueA']=np.nan # slant model fvc = trueA*(week-week0)+fvc0 where trueA is solved as A of y=Ax+B\n        df['trueB']=np.nan\n        \n        df['trueAR']=np.nan # Ratio model fvc = fvc0*trueAR*(week-week0)+fvc0\n        df['trueA0']=np.nan # week0 model fvc = trueA0*(week-week0)+fvc0 かならず ( weeek0, fvc0 )を通る\n        \n        print(CNF['calcAmethod'])\n        \n        for i, p in enumerate(df.Patient.unique()):\n        \n            if CNF['calcAmethod']=='last3':\n                \n                sub = df.loc[(df.Patient == p)&(df.is_last3 == 1), :] \n                fvc = sub.FVC.values - sub.FVC0\n                c = np.vstack([sub.rWeeks.values]).T\n                df.loc[df.Patient == p, 'trueA'] = np.linalg.lstsq(c, fvc, rcond=None)[0]\n\n            elif CNF['calcAmethod']=='simpleLinearRegress':\n                \n                sub = df.loc[(df.Patient == p), :] \n                fvc = sub.FVC.values \n                c = np.vstack([sub.Weeks.values,np.ones(len(sub.Weeks.values))]).T\n                df.loc[df.Patient == p, 'trueA'],df.loc[df.Patient == p, 'trueB'] = np.linalg.lstsq(c, fvc, rcond=None)[0]\n                \n                df.loc[df.Patient == p, 'w0'] = 1/np.sqrt(np.mean(np.abs(sub.FVC.values -  df.loc[df.Patient == p,'trueA'].values*sub.Weeks.values-df.loc[df.Patient == p, 'trueB'].values  )))\n                \n                \n                # AR Retio model\n                c = np.vstack([sub.rWeeks.values]).T\n                df.loc[df.Patient == p, 'trueAR'] = np.linalg.lstsq(c, (fvc-sub.FVC0.values)/sub.FVC0.values, rcond=None)[0]\n                \n                df.loc[df.Patient == p, 'w0'] = 1/np.sqrt(np.mean(np.abs(sub.FVC.values -  df.loc[df.Patient == p,'FVC0'].values*df.loc[df.Patient == p,'trueAR'].values*sub.Weeks.values )))\n                \n                # A0  model\n                c = np.vstack([sub.rWeeks.values]).T\n                df.loc[df.Patient == p, 'trueA0'] = np.linalg.lstsq(c, (fvc-sub.FVC0.values), rcond=None)[0]\n                \n                df.loc[df.Patient == p, 'w0'] = 1/np.sqrt(np.mean(np.abs(sub.FVC.values -  df.loc[df.Patient == p,'trueA'].values*sub.Weeks.values  )))\n                \n                \n            elif CNF['calcAmethod']=='huber':\n                \n                sub = df.loc[(df.Patient == p), :] \n                linear = HuberRegressor().fit(sub.Weeks.values[:,np.newaxis],sub.FVC.values )\n                \n                df.loc[df.Patient == p, 'trueA'] = linear.coef_\n                df.loc[df.Patient == p, 'trueB'] = linear.intercept_\n                \n            else:\n                sub = df.loc[(df.Patient == p), :] \n                fvc = sub.FVC.values - sub.FVC0\n                c = np.vstack([sub.rWeeks.values]).T\n                df.loc[df.Patient == p, 'trueA'] = np.linalg.lstsq(c, fvc, rcond=None)[0]\n    \n\n    for XorY in ['X','Y']:# ------------ vertical features\n        \n        if is_test:\n            vert_df,_vertical_feat_cols = ceateVerticalVolumeFeatDf(f\"./test_vertical{XorY}\",prefixXorY=XorY)            \n        else:\n            vert_df,_vertical_feat_cols = ceateVerticalVolumeFeatDf(f\"{'/'.join(mskDfPath.split('/')[:-1])}/vertical{XorY}/vertical{XorY}\",prefixXorY=XorY)\n        \n        df=df.merge(vert_df ,how='left',on='Patient')\n        feat_cols+=_vertical_feat_cols\n\n    return df,feat_cols\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def markLast3(df):\n    \n    # trainはweek順にきちんと整列してある（確認済）ことを利用して、各Patientの最初の最後の３つだけ残す\n    df['nb']=1\n    df['nb'] = df.groupby('Patient')['nb'].transform('cumsum')\n    df['max_nb'] = df.groupby('Patient')['nb'].transform('max')\n    df['is_last3'] = (df['nb']-df['max_nb']>-3).astype(int)\n    df.drop(['nb','max_nb'], axis=1, inplace=True)\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"device = torch.device(\"cuda:0\") if torch.cuda.is_available() else torch.device(\"cpu\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## OsicDataset for NN"},{"metadata":{"trusted":true},"cell_type":"code","source":"# 患者ごとに　画像スタック、特徴量、ラベル（A）を返す\n\nfrom createClassicSegmented import separate_lungs\n\nclass OsicDataset4NN(Dataset):\n    \n    def __init__(self,table_df,mskDfPath,target_col,is_train=True,num_slice=12,dbg=False,tta_id=0,image_path='train',transform=None,mixup=False):\n        \n        # tta_idはスライスの最初の番号、全スライスが100で10分割するなら0～9\n        \n        assert (image_path=='train')|(image_path=='test')\n        assert transform!=None\n        \n        self.table_df = table_df.copy() ; self.dbg = dbg ; self.num_slice = num_slice \n        self.is_train = is_train ; self.tta_id = tta_id ; self.image_path = image_path \n        self.transform = transform\n        self.mixup = mixup\n        self.target_col = target_col\n        \n        if 'cat_coef' in self.table_df.columns:\n            self.table_df['cat_coef']=self.table_df['cat_coef']/50\n        \n        self.p2fids = self.mkPatient2Files(mskDfPath,num_slice=self.num_slice) # 有効肺画像のdict\n        self.dcmlist=[]\n        self.curryMask={}\n        \n        if self.dbg:\n            allDir = [\"ID00419637202311204720264\",\"ID00132637202222178761324\",\"ID00009637202177434476278\"]\n        else:\n            \n            allDir = self.table_df.Patient.unique()\n        \n        for Patient in allDir:\n            \n            paths=[]\n            jpgpaths=[]\n            for cnt,fids in enumerate(self.p2fids[Patient]):\n\n                paths.append( [ f\"../input/osic-pulmonary-fibrosis-progression/{self.image_path}/{Patient}/{fid}.dcm\" for fid in fids ] )\n                \n                #\n                # CNFがグローバル変数であることに注意\n                #\n                \n                jpgpaths.append( [ f\"{CNF['jpgPath']}/{Patient}/{fid}.jpg\" for fid in fids ] )\n                \n                dcm = pydicom.read_file(f\"../input/osic-pulmonary-fibrosis-progression/{self.image_path}/{Patient}/{fids[0]}.dcm\")\n                \n                # カレーちゃんのマスク\n                edge_pixel_value = dcm.pixel_array[0, 0]\n                if cnt==0:\n                   self.curryMask[Patient]=(dcm.pixel_array != edge_pixel_value)\n                else:\n                   self.curryMask[Patient]|=(dcm.pixel_array != edge_pixel_value)\n    \n            self.dcmlist.append({'Patient':Patient,'paths':paths,'jpgpaths':jpgpaths})\n        \n        #\n        # jpgがない場合は作る\n        #\n\n        for dcm_info in self.dcmlist:\n            \n            for _jpgpaths in dcm_info['jpgpaths']:\n                \n                jpgpath = _jpgpaths[0]\n                file_id = int(jpgpath.split('/')[-1][:-4])\n                p = dcm_info['Patient']\n\n                if os.path.isfile(jpgpath)==False: \n                    \n                    _dir = '/'.join( jpgpath.split('/')[:-1] )\n                    os.makedirs(f\"{_dir}\", exist_ok=True)\n\n                    dcm = pydicom.read_file(f\"../input/osic-pulmonary-fibrosis-progression/{self.image_path}/{p}/{file_id}.dcm\")\n                    \n                    _img,_,_ = self.getNormalized4LungMaskImageFromDCM(dcm,self.curryMask[p],out_size=512)\n\n                    _img = separate_lungs(_img) \n                    _img = (_img-(-1500))/(600-(-1500))\n                    _img = ((255*_img)).clip(0,255).astype('uint8')\n\n                    cv2.imwrite(jpgpath,_img,[cv2.IMWRITE_JPEG_QUALITY, CNF['quality']])\n    \n    def __len__(self):\n        return len(self.dcmlist)\n\n    def seq_aug(self,images):\n        \n        # とりあえずDCMもJPGもtransformするか\n        \n        #DCM\n        _arg ={};_arg['image']=images[0,0]\n        for _i in range(1,CNF['num_sequence']):\n            _arg[f\"image{_i-1}\"]=images[0,_i]\n        \n        #JPG\n        for cnt,_i in enumerate(range(CNF['num_sequence'],2*CNF['num_sequence'])):\n            _arg[f\"image{_i-1}\"]=images[1,cnt]\n        \n        _imgs = self.transform(**_arg)\n        \n        # 最初がDCMで次がJPGの順番\n        transformed_images =  np.stack( [ _imgs['image'] ] + [_imgs[f\"image{_i-1}\"] for _i in range(1,2*CNF['num_sequence']) ] )\n        transformed_images = np.stack([transformed_images[:CNF['num_sequence']],transformed_images[CNF['num_sequence']:]  ])\n        \n        return transformed_images\n        \n    def __getitem__(self, idx):\n        \n        idx = [idx] if type(idx)==int else idx\n        idx = idx.tolist() if torch.is_tensor(idx) else idx\n        # augmentaion \n        \n        images = np.stack([ self.seq_aug(self.readImages(self.dcmlist[_idx]))  for _idx in idx ])\n        \n        Patients = [ self.dcmlist[_idx]['Patient'] for _idx in idx]\n        tables = pd.concat([ self.table_df[self.table_df.Patient==p] for p in Patients ])\n        tables.drop(columns=['Patient'],inplace=True)\n        labels = tables.pop(self.target_col)\n        \n        return torch.FloatTensor(images),torch.FloatTensor(tables.values.astype('float')),torch.FloatTensor(labels.values.astype('float'))\n    \n    def _window_image4LungMask(self,img, img_min, img_max,out_size=256):\n   \n        img[img<img_min] = img_min\n        img[img>img_max] = img_max\n        \n        # LungMask 特有\n        if out_size!=512: # out_size==512 つまりseg用の時は正規化をしない\n            img = np.divide( ( img + 1024 ), 1624 ) # -1024,600でclipした場合は0-1に正規化することになる\n        \n        if self.transform:            # ---- 一枚ごとに異なるaug\n             img = self.transform(image=img)['image']\n        else:\n            img = cv2.resize(img, (out_size, out_size))\n            \n        if out_size==512: # seg用の時はint化\n            img = np.int16(img)\n\n        return img\n\n    def getNormalized4LungMaskImageFromDCM(self,dcm,curry_mask,out_size=256):\n        \n        img = apply_modality_lut(dcm.pixel_array, dcm)\n        \n        # 異常値対応\n        if dcm.PatientID in [\"ID00026637202179561894768\",\"ID00128637202219474716089\"]:\n            img=img+1024  \n        elif dcm.PatientID in [\"ID00132637202222178761324\"]:\n            img=img+4096-2048\n        else:\n            pass\n        \n        # カレーちゃんのクロップ\n        img = img[np.ix_(curry_mask.any(1),curry_mask.any(0))]\n        \n        croppedImageSize = img.shape\n        \n        img = self._window_image4LungMask(img,-1024,600,out_size=out_size)\n        \n        dxyz = np.array([float(dcm.PixelSpacing[0]),float(dcm.PixelSpacing[1]),float(dcm.SliceThickness)])\n        dxyz[0]*=(croppedImageSize[0]/out_size)\n        dxyz[1]*=(croppedImageSize[1]/out_size)\n        \n\n        return img.astype('float'), croppedImageSize, dxyz\n\n    def readImages(self,dcmInfo):\n        \n        Patient = dcmInfo['Patient']\n\n        imgs=[]\n        jpgimgs=[]\n        \n        if self.is_train:\n            \n            # random choice for train とはいってもスライス順番は守る\n            seq_path = [random.choice(paths) for paths in dcmInfo['paths']]\n        \n        else:\n            seq_path = [paths[self.tta_id%len(paths)] for paths in dcmInfo['paths']]\n        \n        for path in seq_path:\n            \n            jpgfname = ( path.split('/')[-1][:-3]+'jpg' )\n            \n            jpgpath = f\"{CNF['jpgPath']}/{Patient}/{jpgfname}\"\n            \n            dcm = pydicom.read_file(path)\n            \n            img, _, _ = self.getNormalized4LungMaskImageFromDCM(dcm,self.curryMask[Patient])\n            \n            img = np.stack(img)\n            jpgimg = np.stack(cv2.resize(cv2.imread(jpgpath, cv2.IMREAD_GRAYSCALE), (256, 256)))\n            \n            imgs.append(img)\n            jpgimgs.append(jpgimg)\n            \n        imgs = np.stack([imgs,jpgimgs])\n            \n        return imgs\n    \n    # 計算済のLungMaskDfから面積の小さい肺スライス画像を除く\n    def mkPatient2Files(self,mskDfPath,num_slice=12):\n\n        _dict={}\n        \n        _df = pd.read_csv(mskDfPath)\n        _df['lung_area'] = _df['lung_area']*_df['dx']*_df['dy']\n        _df = _df.drop(index=_df.loc[(_df.lung_area<400)].index).reset_index(drop=True) \n        \n        # 2回繰り返しも片方だけにする\n        if self.is_train:\n            _df = _df.drop(index=_df.loc[(_df.Patient=='ID00078637202199415319443')&(_df.file_id>509)].index).reset_index(drop=True) # remove duplicated sequence \n        \n        \n        for p in _df.Patient.unique():\n\n            file_ids = np.sort(_df.loc[_df.Patient==p,'file_id'].values)\n            \n            if len(file_ids)<num_slice:\n                _dict[p] = [ [fid] for fid in file_ids[np.linspace(0,len(file_ids)-1,num=num_slice).round(0).astype('int')] ]\n            \n            else:\n                _dict[p] = [ file_ids[vid] for tid,vid in KFold(n_splits=num_slice).split(file_ids) ] # KFoldでfile_idsをnum_sliceに分割\n        \n        # sub時はたくさんsegできないので12スライスのみsegしてjpg生成しとく\n        if self.image_path=='test': \n            \n            print(\"submit mode : dcm is gonna be reducedto 12\")\n            \n            for p in _df.Patient.unique():\n                _dict[p] =  [ np.array([_paths[int(len(_paths)/2)]]) for _paths in _dict[p] ]\n\n        return _dict\n\n    def collate_fn(self,batch):\n        \n        if self.mixup:\n            \n            # 最後にmixupを付け加える\n            images, tables, labels = list(zip(*batch))\n            \n            mix_image = (images[-1]+images[-2])/2.0\n            mix_table = (tables[-1]+tables[-2])/2.0\n            mix_label = (labels[-1]+labels[-2])/2.0\n            \n            images = torch.stack(images+(mix_image,)).squeeze(1)\n            tables = torch.stack(tables+(mix_table,)).squeeze(1)\n            labels = torch.stack(labels+(mix_label,)).squeeze(1)\n            \n        else:\n            \n            images, tables, labels = list(zip(*batch))\n            images = torch.stack(images).squeeze(1)\n            tables = torch.stack(tables).squeeze(1)\n            labels = torch.stack(labels).squeeze(1)\n\n        return images, tables, labels    \n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cusutom Neural Network"},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nfrom torch import nn\nimport torch.nn.functional as F\n\nclass OSICNetR231(nn.Module):\n    \n    def __init__(self,num_sequence,num_table_ft,table_only=False,initialize_with_r231=True):\n    \n        super(OSICNetR231, self).__init__()\n        \n        self.dual = True\n        \n        self.table_only = table_only\n        if table_only:\n            print(f\"TABLE ONLY MODE\")\n\n        self.base = UNet(n_classes=3, padding=True, depth=5, up_mode='upsample', batch_norm=True, residual=False)\n        if initialize_with_r231:\n            self.base.load_state_dict(torch.load('../input/pthunetlungmask/unet_r231-d5d2fc3d.pth'))\n        \n        del self.base.softmax, self.base.last, self.base.up_path\n        \n        \n        self.lstm = nn.LSTM(2*1024, CNF['lstm_mem'], 1,bidirectional=True) \n        \n        self.dropout = nn.Dropout(p=0.5)\n        \n        if self.table_only:\n            self.linear = nn.Linear(0+num_table_ft, 1)\n            self.bn0    = nn.BatchNorm1d(0+num_table_ft)\n        else:\n            self.linear = nn.Linear(2*num_sequence+num_table_ft, 1)\n            self.bn0    = nn.BatchNorm1d(2*num_sequence+num_table_ft)\n\n            \n    def forward(self, dual_images, table_ft):\n        \n        # dual_images[:,0]:dcm ,dual_images[:,1]:classic_segmentated\n        \n        images = torch.stack( (  dual_images[:,0,0,:,:], dual_images[:,1,0,:,:],\n                                 dual_images[:,0,1,:,:], dual_images[:,1,1,:,:],\n                                 dual_images[:,0,2,:,:], dual_images[:,1,2,:,:],\n                                 dual_images[:,0,3,:,:], dual_images[:,1,3,:,:],\n                                 dual_images[:,0,4,:,:], dual_images[:,1,4,:,:],\n                                 dual_images[:,0,5,:,:], dual_images[:,1,5,:,:],\n                                 dual_images[:,0,6,:,:], dual_images[:,1,6,:,:],\n                                 dual_images[:,0,7,:,:], dual_images[:,1,7,:,:],\n                                 dual_images[:,0,8,:,:], dual_images[:,1,8,:,:],\n                                 dual_images[:,0,9,:,:], dual_images[:,1,9,:,:],\n                                 dual_images[:,0,10,:,:], dual_images[:,1,10,:,:],\n                                 dual_images[:,0,11,:,:], dual_images[:,1,11,:,:]) )\n        \n        images = images.permute(1, 0, 2,3)\n        images = images.unsqueeze(2) \n\n        seq=[]\n        for x in images:\n\n            blocks=[]\n            for i, down in enumerate(self.base.down_path):\n\n                x = down(x)\n\n                if i != len(self.base.down_path) - 1:\n                    blocks.append(x)\n                    x = F.avg_pool2d(x, 2)\n\n                # この時点で[1024, 16, 16]\n\n            x = torch.cat( ( F.avg_pool2d(x, (16,16)).squeeze(),F.max_pool2d(x, (16,16)).squeeze() ),dim=-1 )\n\n            seq.append(x)\n\n        lstm_input = torch.stack(seq)\n        \n        lstm_input = lstm_input.squeeze(2)\n        lstm_out,_ = self.lstm(lstm_input)\n        image_ft = torch.mean(lstm_out,-1) # seqごとの特徴量の平均値\n        \n        if self.table_only:\n        \n            ft = self.bn0(table_ft)\n        \n        else:\n            \n            ft = torch.cat([image_ft,table_ft],1)\n            ft = self.bn0(ft)\n\n        output=self.linear(ft).squeeze(1)\n        \n        \n        return output\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## FVCFitter using LGBM ( includes Weeks as feature )"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\nclass FVCFitterWeekLGBM():\n    \n    # Aを推定するLGBM\n    \n    def __init__(self,CNF,fold_id,prefix=\"osicLGBMweek\"):\n        \n        self.target_col = 'FVC'\n        print(f\"Estimation mid target {self.target_col}\")\n        \n        self.CNF = CNF\n        self.fold_id = fold_id\n        self.reg = lgb.LGBMRegressor(objective='l1', num_leaves=5, importance_type='gain') \n        self.prefix = prefix\n        \n        \n    # train \n    def fit(self,trn_df,val_df,target):\n        \n        seed_everything(2020)\n        \n        # Estimate A, so we don't need Weeks\n        \n        val_df = val_df[['Patient']+self.CNF['feat_cols']+[self.target_col]+['Weeks']].reset_index(drop=True)\n        \n            \n        trn_df = trn_df[['Patient']+self.CNF['feat_cols']+[self.target_col]+['Weeks']].reset_index(drop=True)\n\n        self.reg.fit(trn_df[self.CNF['feat_cols']+['Weeks']],trn_df[self.target_col],\n                     eval_set=(val_df[self.CNF['feat_cols']+['Weeks']],val_df[self.target_col]),early_stopping_rounds=100,verbose=False)\n\n        _save(self.reg,f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pkl\")\n        \n        del self.reg\n        \n    def predict(self,tst_df,image_path):\n        \n        self.reg = _load(f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pkl\")\n        \n        seed_everything(2020)\n        \n        _df = tst_df[['Patient']+self.CNF['feat_cols']+['Weeks']].reset_index(drop=True).copy()\n        _df[self.target_col] = np.nan # dummmy\n        \n        pred_df = _df[['Patient']].copy()\n        \n        pred_df['_predFVC'] = self.reg.predict(_df[self.CNF['feat_cols']+['Weeks']])\n        \n        gc.collect()\n        \n        return pred_df['_predFVC'].values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## FVCFitter using LGBM"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\nclass FVCFitterLGBM():\n    \n    # Aを推定するLGBM\n    \n    def __init__(self,CNF,fold_id,prefix=\"osicLGBM\"):\n        \n        if ('A0' in prefix)|('a0' in prefix):\n            self.target_col = 'trueA0'\n        elif ('AR' in prefix)|('ar' in prefix):\n            self.target_col = 'trueAR'\n        else:\n            self.target_col = 'trueA'\n            \n        if 'w0_' in prefix:\n            self.weight_col = 'w0'\n        else:\n            self.weight_col = None\n            \n        assert self.target_col in ['trueA','trueAR','trueA0']\n        \n        print(f\"Estimation mid target {self.target_col}\")\n        \n        self.CNF = CNF\n        self.fold_id = fold_id\n        self.reg = lgb.LGBMRegressor(objective='l1', num_leaves=5, importance_type='gain') \n        self.prefix = prefix\n        \n        \n    # train \n    def fit(self,trn_df,val_df,target):\n        \n        seed_everything(2020)\n        \n        # Estimate A, so we don't need Weeks\n        \n        val_df = val_df[['Patient']+self.CNF['feat_cols']+[self.target_col]].drop_duplicates().reset_index(drop=True)\n        \n        if self.weight_col!=None:\n            \n            trn_df = trn_df[['Patient']+self.CNF['feat_cols']+[self.target_col]+[self.weight_col]].drop_duplicates().reset_index(drop=True)\n            \n            self.reg.fit(trn_df[self.CNF['feat_cols']],trn_df[self.target_col],sample_weight=trn_df[self.weight_col],\n                         eval_set=(val_df[self.CNF['feat_cols']],val_df[self.target_col]),early_stopping_rounds=100,verbose=False)\n        else:\n            \n            trn_df = trn_df[['Patient']+self.CNF['feat_cols']+[self.target_col]].drop_duplicates().reset_index(drop=True)\n            \n            self.reg.fit(trn_df[self.CNF['feat_cols']],trn_df[self.target_col],\n                         eval_set=(val_df[self.CNF['feat_cols']],val_df[self.target_col]),early_stopping_rounds=100,verbose=False)\n\n        _save(self.reg,f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pkl\")\n        \n        del self.reg\n        \n    def predict(self,tst_df,image_path):\n        \n        self.reg = _load(f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pkl\")\n        \n        seed_everything(2020)\n        \n        _df = tst_df[['Patient']+self.CNF['feat_cols']].drop_duplicates().reset_index(drop=True).copy()\n        _df[self.target_col] = np.nan # dummmy\n        \n        pred_df = _df[['Patient']].copy()\n        \n        pred_df['predA'] = self.reg.predict(_df[self.CNF['feat_cols']])\n        \n        gc.collect()\n        \n        _tst_df = tst_df.merge(pred_df,on='Patient').copy()\n        \n        if self.target_col=='trueAR':\n            return ( _tst_df['FVC0'] + _tst_df['FVC0']*_tst_df['predA']*(_tst_df['Weeks']-_tst_df['Week0']) ).values\n        else: \n            return ( _tst_df['FVC0'] + _tst_df['predA']*(_tst_df['Weeks']-_tst_df['Week0']) ).values  # 'trueA','trueA0' use same eq\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## FVCFitter using custom Neural Network"},{"metadata":{"trusted":true},"cell_type":"code","source":"class FVCFitterNN01():\n    \n    # Aを推定するNN\n    \n    def __init__(self,CNF,fold_id,prefix=\"osicNN\"):\n        \n        if ('A0' in prefix)|('a0' in prefix):\n            self.target_col = 'trueA0'\n        elif ('AR' in prefix)|('ar' in prefix):\n            self.target_col = 'trueAR'\n        else:\n            self.target_col = 'trueA'\n            \n        assert self.target_col in ['trueA','trueAR','trueA0']\n        \n        print(f\"Estimation mid target {self.target_col}\")\n        \n        self.CNF = CNF\n        self.fold_id = fold_id\n        self.prefix = prefix\n        self.num_table_ft = len(CNF['feat_cols'])\n        \n        from albumentations import Compose, Resize, HorizontalFlip,RandomBrightnessContrast,ShiftScaleRotate,RandomResizedCrop,VerticalFlip\n\n        self.trnTransform = Compose([ Compose([HorizontalFlip(always_apply=False, p=1),VerticalFlip(always_apply=False, p=1)],p=0.1),#時々180度回転してるやつがある\n                                 ShiftScaleRotate(shift_limit=0.0625, scale_limit=(0.1,0.3), rotate_limit=10, interpolation=1, border_mode=cv2.BORDER_CONSTANT, value=0, mask_value=None,p=1),\n                                 RandomResizedCrop(256, 256, scale=(1, 1), ratio=(1.0, 1.0), interpolation=1, always_apply=True, p=1.0), # アスペクト１じゃない場合はセンターでアスペクト１でクロップ \n                               ],additional_targets = { f\"image{_i}\":'image' for _i in range(0,CNF['num_sequence']-1) } )\n\n        self.valTransform = Compose([ ShiftScaleRotate(shift_limit=0, scale_limit=(0.2,0.2), rotate_limit=0, interpolation=1, border_mode=cv2.BORDER_CONSTANT, value=0, mask_value=None,p=1),\n                                 RandomResizedCrop(256, 256, scale=(1, 1), ratio=(1.0, 1.0), interpolation=1, always_apply=True, p=1.0),\n                               ],additional_targets = { f\"image{_i}\":'image' for _i in range(0,CNF['num_sequence']-1) })\n\n        self.noTransform = Compose([ RandomResizedCrop(256, 256, scale=(1, 1), ratio=(1.0, 1.0), interpolation=1, always_apply=True, p=1.0),\n                               ],additional_targets = { f\"image{_i}\":'image' for _i in range(0,CNF['num_sequence']-1) })\n        \n    # train NN then save pth\n    def fit(self,trn_df,val_df,target):\n        \n        if not self.CNF['isTrain']:\n            print('isTrain is set False,skip training')\n            return True\n        \n        print(f\"preparing data\")\n        \n        seed_everything(2020)\n        \n        # Estimate A, so we don't need Weeks\n        trn_df = trn_df[['Patient']+self.CNF['feat_cols']+[self.target_col]].drop_duplicates().reset_index(drop=True)\n        val_df = val_df[['Patient']+self.CNF['feat_cols']+[self.target_col]].drop_duplicates().reset_index(drop=True)\n        \n        # mixup有効時はcollate_fun内でmixupする5+1\n        trn_dset    = OsicDataset4NN(trn_df,self.CNF['mskDfPath'],target_col=self.target_col,is_train=True,num_slice=self.CNF['num_sequence'],dbg=False,image_path='train',transform=self.trnTransform if self.CNF['aug'] else self.noTransform,mixup=self.CNF['mixup'])\n        trn_dloader = DataLoader(trn_dset, batch_size=self.CNF['batch_size']-int(self.CNF['mixup']),  shuffle=True, num_workers=2,collate_fn=trn_dset.collate_fn,drop_last=True) \n\n        val_dset    = OsicDataset4NN(val_df,self.CNF['mskDfPath'],target_col=self.target_col,is_train=False,num_slice=self.CNF['num_sequence'],dbg=False,image_path='train',tta_id=0,transform=self.valTransform if self.CNF['aug'] else self.noTransform,mixup=False)\n        val_dloader = DataLoader(val_dset, batch_size=self.CNF['batch_size'],  shuffle=False, num_workers=2,collate_fn=val_dset.collate_fn,drop_last=False,)\n\n        model = OSICNetR231(num_sequence=self.CNF['num_sequence'],num_table_ft=self.num_table_ft,table_only=False).to(device)\n        plist = [{'params': model.parameters(), 'lr': 0.001}]\n        optimizer = optim.Adam(plist, lr=0.001) \n        \n        # ---------------------------- comment out for MSE LOSS -------------------------------------\n        \n        criterion = torch.nn.L1Loss()\n            \n        val_criterion = torch.nn.L1Loss()\n        \n        if self.CNF['EMA']:\n            ema = ModelEMA(model)\n\n        trnLoss=[]\n        valLoss=[]\n\n        best_vloss=np.inf\n        \n        print(f\"start training\")\n        \n        for ep in range(2):\n        #for ep in range(self.CNF['maxEpoch']):\n\n            model.train()\n\n            tloss=[]\n            for images,tables,labels in trn_dloader:\n                \n                images = images.to(device); tables = tables.to(device); labels = labels.to(device)\n                \n                pred = model(images,tables)\n                loss = criterion(pred,labels)\n                loss.backward()\n                optimizer.step()\n                optimizer.zero_grad()\n                if self.CNF['EMA']:\n                    ema.update(model)\n\n                tloss.append(loss.detach().cpu().numpy())\n            \n            trnLoss.append(np.mean(tloss))\n\n            model.eval()\n\n            with torch.no_grad():\n\n                vloss=[]\n                for images,tables,labels in val_dloader:\n\n                    images = images.to(device); tables = tables.to(device); labels = labels.to(device)\n\n                    pred = model(images,tables)\n                    loss = val_criterion(pred,labels)\n\n                    vloss.append(loss.detach().cpu().numpy())\n\n            if self.CNF['verbose']:\n                print(f\"EP{ep:02d} trnLoss:{np.mean(tloss):.4f} valLoss:{np.mean(vloss):.4f}\")\n            \n            valLoss.append(np.mean(vloss))\n\n            if ep>0:\n\n                if valLoss[-1] < best_vloss:\n\n                    best_vloss = valLoss[-1]\n                    os.makedirs(f\"{self.prefix}\", exist_ok=True)\n\n                    torch.save(model.state_dict(), f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pth\")\n\n\n        del model, plist, optimizer, criterion,trn_dset,val_dset,trn_dloader,val_dloader\n        gc.collect()\n        torch.cuda.empty_cache()\n        \n        print(f\"FOLD{self.fold_id} {np.min(valLoss):.4f}@ep{np.argmin(valLoss)}\")            \n        \n        plt.plot(trnLoss,label='trn')\n        plt.plot(valLoss,label='val')\n        plt.xlabel('epoch')\n        plt.ylabel('L1 loss')\n        plt.title(f\"Fold{self.fold_id} best_vloss {best_vloss:.4f}@ep{np.argmin(valLoss)}\")\n        plt.legend()\n        plt.show()\n        \n            \n    def predict(self,tst_df,image_path):\n        \n        print(\"start prediction\")\n        seed_everything(2020)\n        \n        model = OSICNetR231(num_sequence=self.CNF['num_sequence'],num_table_ft=self.num_table_ft,table_only=False).to(device)\n        \n        _df = tst_df[['Patient']+self.CNF['feat_cols']].drop_duplicates().reset_index(drop=True).copy()\n        \n        _df[self.target_col] = np.nan # dummmy\n        \n        pred_df = _df[['Patient']].copy()\n        pred_df['predA'] = np.nan\n        \n        \n\n        if self.CNF['isTrain']:\n            model.load_state_dict(torch.load(f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pth\"))\n        else:\n            \n            #\n            # awsの場合はコメントアウト\n            #\n            capital_prefix = self.prefix.replace('a0','A0').replace('ar','AR')\n            model.load_state_dict(torch.load(f\"../input/{self.prefix}/{capital_prefix}_model_fold{self.fold_id}.pth\"))\n            \n            #model.load_state_dict(torch.load(f\"../input/{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pth\"))\n\n        model.eval()\n\n        pred_over_tta=[]\n        \n        for tta_id in range(CNF['n_tta']):\n            \n            tst_dset    = OsicDataset4NN(_df,self.CNF['mskDfPath'],target_col=self.target_col,is_train=False,num_slice=self.CNF['num_sequence'],dbg=False,image_path=image_path,tta_id=tta_id,transform=self.valTransform if CNF['aug'] else self.noTransform,mixup=False)\n            tst_dloader = DataLoader(tst_dset, batch_size=self.CNF['batch_size'],  shuffle=False, num_workers=2,collate_fn=tst_dset.collate_fn,drop_last=False)\n\n            fold_pred=[]\n            with torch.no_grad():\n\n                for images,tables,labels in tst_dloader:\n                    \n                    images = images.to(device); tables = tables.to(device); labels = labels.to(device)\n\n                    pred = model(images,tables)\n                    pred=pred.cpu().numpy()\n                    fold_pred.append(pred)\n                    \n\n            pred_over_tta.append( np.concatenate(fold_pred) )\n            \n            del tst_dset,tst_dloader\n            gc.collect()\n\n        pred_df['predA'] = np.mean(pred_over_tta,axis=0)\n        \n        del model\n        gc.collect()\n        torch.cuda.empty_cache()\n        \n        _tst_df = tst_df.merge(pred_df,on='Patient').copy()\n        \n        if self.target_col == 'trueAR':\n        \n            return ( _tst_df['FVC0'] + _tst_df['predA']*(_tst_df['Weeks']-_tst_df['Week0']) ).values\n            \n            \n        else: # 'trueA','trueA0' use same eq\n        \n            return ( _tst_df['FVC0'] + _tst_df['predA']*(_tst_df['Weeks']-_tst_df['Week0']) ).values\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load train and add features"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv') \ntrain = markLast3(train)\n# dupおとして\ntrain = train.drop_duplicates(subset=['Patient','Weeks'],keep='last')\ntrain = train.reset_index(drop=True)\ntrain, feat_cols = addFeatureDf(train,CNF['mskDfPath'],CNF['trnPath'],is_test=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Single Models Fit and Predict"},{"metadata":{"trusted":true},"cell_type":"code","source":"import warnings\nwarnings.simplefilter('ignore')\n\nseed_everything(2020)\n\nfeature_importance_df = pd.DataFrame()\n\nnMegaFold = 5  \n\nCNF['feat_cols']=feat_cols\n\ntrain_orig = train.copy()\n\ntrainModelGroup = [ train.copy() for _ in range(len(CNF['prefixes'])) ]\n\n\nfor prefixsModelGroup,weightsModelGroup,train in zip(CNF['prefixes'],CNF['ensembleWeight'],trainModelGroup):\n    \n    _isImportancePlot=False\n    \n    for prefix in prefixsModelGroup+['']:\n        train['predFVC'+prefix] = 0\n\n    for prefix,ensemble_weight in zip(prefixsModelGroup,weightsModelGroup):\n\n        if CNF['isTrain']|('LGBM' in prefix):\n\n            random_seed = int(prefix.split('seed')[1])\n\n            print(f\"SEED:{random_seed}\")\n\n            unique_patients = train.Patient.unique()\n            valPatSplit5 = [ unique_patients[val_patient] for trn_patient,val_patient in KFold(n_splits=5,shuffle=True,random_state=random_seed).split(unique_patients) ]\n\n            os.makedirs(f\"{prefix}\", exist_ok=True)\n            _save(valPatSplit5,f\"{prefix}/valPatSplit5.pkl\")\n\n        else:\n\n            valPatSplit5 = _load(f\"../input/{prefix}/valPatSplit5.pkl\")\n\n        #\n        # train or validate FVCfitter\n        #\n\n        for megaFold in range(nMegaFold):\n\n            print(f\"megaFold:{megaFold}\")\n\n            patTrn = np.concatenate([ valPatSplit5[(_i+megaFold)%5] for _i in range(4) ])\n            patVal = np.concatenate([ valPatSplit5[(_i+megaFold)%5] for _i in [4] ])\n\n            dataTrn,dataVal = train[train.Patient.isin(patTrn)],train[train.Patient.isin(patVal)]\n\n\n            # train fvc_fitter on dataTrn\n            if 'LGBM' in prefix:\n                if 'week' in prefix:\n                    fvc_fitter = FVCFitterWeekLGBM(CNF,fold_id=megaFold,prefix=prefix)\n                else:\n                    fvc_fitter = FVCFitterLGBM(CNF,fold_id=megaFold,prefix=prefix)\n                \n            else:\n                fvc_fitter = FVCFitterNN01(CNF,fold_id=megaFold,prefix=prefix)\n\n            fvc_fitter.fit(trn_df=dataTrn, val_df=dataVal, target=dataTrn['FVC']) \n\n            _pred = fvc_fitter.predict(dataVal.drop(columns='FVC'),image_path='train')\n\n            train.loc[train.Patient.isin(patVal),'predFVC'] += ensemble_weight*_pred\n            train.loc[train.Patient.isin(patVal),'predFVC'+prefix] += _pred\n            \n            if ('LGBM' in prefix)&('week' not in prefix):\n                _isImportancePlot=True\n                fold_importance_df = pd.DataFrame()\n                fold_importance_df[\"Feature\"] = feat_cols\n                fold_importance_df[\"importance\"] = fvc_fitter.reg.feature_importances_\n                fold_importance_df[\"fold\"] = f\"{prefix}_fold{megaFold}\"\n                feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n\n\n            del fvc_fitter\n            gc.collect()\n        \n    if _isImportancePlot:\n        plt.figure(figsize=(8, 20))\n        sns.barplot(x=\"importance\", y=\"Feature\", data=feature_importance_df.sort_values(by=\"importance\", ascending=False));","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train and predict preliminary confidence value for just calculate oof "},{"metadata":{"trusted":true},"cell_type":"code","source":"#\n# train and validate Pre-Confidenvefitter\n#\n\nscore_dfs = []\n\nfor prefixsModelGroup,weightsModelGroup,train in zip(CNF['prefixes'],CNF['ensembleWeight'],trainModelGroup):\n    \n    scores = {}\n    ensembleW = {}\n\n    conf_feat = ['mskd_img_sum_cls4_mean','Male'] # 仮のconf feat\n\n\n    for prefix in prefixsModelGroup+['']: # 必ず最後に''を入れること じゃないとtestの時にバグるはず\n\n        pred_df = train[['rWeeks','FVC','predFVC'+prefix,'is_last3']]\n        pred_df['dFVC'+prefix]=pred_df['predFVC'+prefix]-pred_df['FVC']\n        err_df = pred_df.groupby('rWeeks')['dFVC'+prefix].apply(lambda x: np.sqrt(np.mean(x**2))).reset_index()\n        err_df.columns = ['rWeeks', 'rmseFVC']\n\n        train['dFVC'+prefix] = abs(train['predFVC'+prefix]-train['FVC'])\n\n        # features correlated to FVC error\n        if prefix=='':\n            #pass\n            print(train.loc[train.is_last3==1].corr()['dFVC'+prefix].sort_values(ascending=False)[:20])\n\n        # features used for confidence fit\n\n\n        confidence_fitter = ConfidenceFitterOptimizedLinear(train.loc[train.is_last3==1,conf_feat].values,train.loc[train.is_last3==1,'FVC'].values,train.loc[train.is_last3==1,'predFVC'+prefix].values)\n\n        #\n        # plot result\n        #\n\n        # FVC vs predFVC\n        pred_df['sigma'+prefix] = confidence_fitter.predict(train[conf_feat])\n        if prefix=='':\n            pass\n            #plt.figure(figsize=(16,5))\n            #plt.subplot(122);plt.plot([500,7000],[500,7000],c='orange',alpha=0.5)\n            #sc=plt.scatter( train.loc[train.is_last3==1,'FVC'],  train.loc[train.is_last3==1,'predFVC'+prefix],c=pred_df.loc[pred_df.is_last3==1,'sigma'+prefix],cmap=cm.seismic,alpha=0.2 );plt.colorbar(sc)\n            #plt.xlabel('FVC');plt.ylabel('predFVC'+prefix);plt.show()\n\n        #\n        # CV\n        #\n\n        if CNF['EvalOnlyLast3']:\n            pred_df = pred_df[ pred_df.is_last3==1 ]\n\n        scores[prefix] = calcScore(pred_df['FVC'].values, pred_df['predFVC'+prefix].values, pred_df['sigma'+prefix].values)\n        \n    _score_df = pd.DataFrame.from_dict(scores, orient='index').rename(index={'': 'Ensemble'}).reset_index()\n    _score_df.columns=['model','score']\n    _score_df['weight']=np.concatenate([weightsModelGroup,np.array([0])])\n\n    _score_df.reindex(columns=['model', 'weight', 'score'])\n    \n    score_dfs.append(_score_df)\n\ndisplay(HorizontalDisplay(*score_dfs))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Last Confidence Fitter used to generate final confidence value"},{"metadata":{"trusted":true},"cell_type":"code","source":"# fold sanity check and get each model group's split\n\nvalPatSplit5ModelGroup = []\n\nfor prefixsModelGroup,weightsModelGroup,train in zip(CNF['prefixes'],CNF['ensembleWeight'],trainModelGroup):\n    \n    all_splits=[]\n    for prefix in prefixsModelGroup:\n        path = f\"{prefix}/valPatSplit5.pkl\" if CNF['isTrain']|('LGBM' in prefix) else f\"../input/{prefix}/valPatSplit5.pkl\"\n        all_splits.append(_load(path))\n\n    sanity_ok=True\n    from itertools import permutations\n    for fold in range(5):\n        for _i,_j in permutations(np.arange(len(prefixsModelGroup)),2):\n            if set(all_splits[_i][fold])!=set(all_splits[_j][fold]):\n                assert False,'invalid fold splits'\n\n    valPatSplit5ModelGroup.append( all_splits[0] )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# find ideal confidence\n\nimport scipy as sp\nfrom tqdm.notebook import tqdm\nfrom functools import partial\nimport math\n\ndef loss_func(weight, row):\n    confidence = weight\n    sigma_clipped = max(confidence, 70)\n    diff = abs(row['FVC'] - row['predFVC'])\n    delta = min(diff, 1000)\n    score = -math.sqrt(2)*delta/sigma_clipped - np.log(math.sqrt(2)*sigma_clipped)\n    return -score\n\nfor train in trainModelGroup:\n    \n    results = [] # ----------------------------------- 正解のConfidenceを作成する\n    #tk0 = tqdm(train.iterrows(), total=len(train))\n    for _, row in train.iterrows():\n        loss_partial = partial(loss_func, row=row)\n        weight = [100]\n        result = sp.optimize.minimize(loss_partial, weight, method='SLSQP')\n        x = result['x']\n        results.append(x[0])\n\n    train['trueConfidence']=results\n    train['trueConfidence']=np.maximum(train['trueConfidence'].values,70)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### check if confidence values are perfectly predicted, how good CV is"},{"metadata":{"trusted":true},"cell_type":"code","source":"for cnt,train in enumerate(trainModelGroup):\n    print(f\"ModelGroup{cnt} Possible CV {calcScore(train['FVC'].values, train['predFVC'].values, train['trueConfidence'].values):.4f}\") # 完璧にConfidenceが推定できた場合","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.linear_model import Ridge,Lasso\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import NuSVR\nfrom sklearn.ensemble import RandomForestRegressor\n\nclass ConfidenceFitter():\n    \n    # trueConfidenceを推定するRidge or SVR\n    \n    def __init__(self,CNF,fold_id,model_name,prefix=\"confRidge\",target_col='trueConfidence'):\n        \n        self.target_col = target_col\n        self.scaler = StandardScaler()\n        \n        self.CNF = CNF\n        self.fold_id = fold_id\n        \n        if model_name == 'ridge':\n            self.reg = Ridge(alpha=1.0) \n        elif model_name == 'svr':\n            self.reg =NuSVR()\n        else:\n            assert False\n            \n        self.prefix = prefix\n        self.feat_cols = None\n        \n    # train \n    def fit(self,trn_df,val_df,feat_cols):\n        \n        seed_everything(2020)\n        self.feat_cols = feat_cols\n        \n        trn_df = trn_df.reset_index(drop=True)\n        val_df = val_df.reset_index(drop=True)\n        \n        self.scaler.fit(trn_df[self.feat_cols])\n        \n        self.reg.fit(self.scaler.transform( trn_df[self.feat_cols] ),  trn_df[self.target_col])\n        \n        os.makedirs(f\"{self.prefix}\", exist_ok=True)\n        _save(self.reg,f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pkl\")\n        _save(self.scaler,f\"{self.prefix}/{self.prefix}_scaler_fold{self.fold_id}.pkl\")\n        _save(self.feat_cols,f\"{self.prefix}/{self.prefix}_feat_fold{self.fold_id}.pkl\")\n        \n        del self.reg\n        del self.scaler\n        \n    def predict(self,tst_df):\n        \n        tst_df= tst_df.reset_index(drop=True).copy()\n        \n        self.reg = _load(f\"{self.prefix}/{self.prefix}_model_fold{self.fold_id}.pkl\")\n        self.scaler = _load(f\"{self.prefix}/{self.prefix}_scaler_fold{self.fold_id}.pkl\")\n        self.feat_cols = _load(f\"{self.prefix}/{self.prefix}_feat_fold{self.fold_id}.pkl\")\n        \n        seed_everything(2020)\n        \n        _df = tst_df[['Patient']+self.feat_cols].copy()\n        \n        _df['predConfidence'] = self.reg.predict(self.scaler.transform(_df[self.feat_cols]))\n        \n        gc.collect()\n        \n        return _df['predConfidence'].values\n        \n        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train last confidence fitter and predict confidence"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nconf_cols = ['mskd_img_sum_cls4_mean','Male','mskd_img_sum_cls2_median_absolute_deviation','Currently smokes','rWeeks']\n\n\nfor cnt,(valPatSplit5,train) in enumerate(zip(valPatSplit5ModelGroup,trainModelGroup)):\n\n    confidence_fitter_prefix=f\"lastconfModelGroup{cnt}\"\n    \n    train['predConfidence']=0\n\n    for fold in range(5):\n\n        patTrn = np.concatenate([ valPatSplit5[(_i+fold)%5] for _i in range(4) ])\n        patVal = np.concatenate([ valPatSplit5[(_i+fold)%5] for _i in [4] ])\n\n        dataTrn,dataVal = train[(train.is_last3==1)&(train.Patient.isin(patTrn))],train[(train.is_last3==1)&(train.Patient.isin(patVal))]\n\n        conf_fitterRidge = ConfidenceFitter(CNF,fold_id=fold,model_name='ridge',prefix=confidence_fitter_prefix+'RIDGE')\n        conf_fitterRidge.fit(trn_df=dataTrn.drop(columns='FVC'), val_df=dataVal.drop(columns='FVC'),feat_cols=conf_cols) \n\n        conf_fitterSVR = ConfidenceFitter(CNF,fold_id=fold,model_name='svr',prefix=confidence_fitter_prefix+'SVR')\n        conf_fitterSVR.fit(trn_df=dataTrn.drop(columns='FVC'), val_df=dataVal.drop(columns='FVC'),feat_cols=conf_cols) \n\n        _predRidge = conf_fitterRidge.predict(dataVal)\n        _predSVR   = conf_fitterSVR.predict(dataVal)\n\n        _pred = np.sqrt( 0.2*(_predSVR**2) + 0.8*(_predRidge**2))\n\n        train.loc[(train.is_last3==1)&(train.Patient.isin(patVal)),'predConfidence'] = _pred\n\n        del conf_fitterRidge,conf_fitterSVR\n        gc.collect()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_group_scores=[]\nfor cnt,(valPatSplit5,train) in enumerate(zip(valPatSplit5ModelGroup,trainModelGroup)):\n\n    plt.figure(figsize=(16,5))\n    plt.subplot(121);plt.plot([500,7000],[500,7000],c='orange',alpha=0.5)\n    sc=plt.scatter( train.loc[train.is_last3==1,'FVC'],  train.loc[train.is_last3==1,'predFVC'],c=train.loc[train.is_last3==1,'predConfidence'],cmap=cm.seismic,alpha=0.2 );plt.colorbar(sc)\n    plt.xlabel('FVC');plt.ylabel('predFVC'+prefix)\n    plt.subplot(122);plt.scatter(train.loc[train.is_last3==1,'predConfidence'],train.loc[train.is_last3==1,'trueConfidence']);plt.xlabel('predConfidence');plt.ylabel('trueConfidence')\n    plt.plot([0,600],[0,600],c='orange');plt.show()\n    \n    model_group_scores.append( calcScore(train.loc[train.is_last3==1,'FVC'].values, train.loc[train.is_last3==1,'predFVC'].values, train.loc[train.is_last3==1,'predConfidence'].values) )\n    print(f\"CV {model_group_scores[-1]}\")\n    \na=np.array(model_group_scores)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# CV"},{"metadata":{"trusted":true},"cell_type":"code","source":"final_oof_df = pd.DataFrame({'predFVC':np.average([train.predFVC for train in trainModelGroup],0),\n              'FVC':np.mean([train.FVC for train in trainModelGroup],0),\n              'predConfidence':np.sqrt( np.average([train.predConfidence**2 for train in trainModelGroup],0)),\n              'is_last3':np.mean([train.is_last3 for train in trainModelGroup],0),\n                        })\n\nprint(f\"CV {calcScore(final_oof_df.loc[final_oof_df.is_last3==1,'FVC'].values, final_oof_df.loc[final_oof_df.is_last3==1,'predFVC'].values, final_oof_df.loc[final_oof_df.is_last3==1,'predConfidence'].values)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()\n\nif CNF['isTrain']==False:\n    \n    !cp {'/'.join(CNF['mskDfPath'].split('/')[:-1])}/createLungMaskFeat.py .\n    !python createLungMaskFeat.py -path ../input/osic-pulmonary-fibrosis-progression/test -out LTRC_feat_test_df.csv -outvertical ./test_vertical\n    \n    # read test csv\n    test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv') \n    \n\n    # add Featreus\n    test,_ = addFeatureDf(test,'./LTRC_feat_test_df.csv','../input/osic-pulmonary-fibrosis-progression/test',is_test=True)\n    \n    test = test.drop(columns=['Weeks','rWeeks'])\n\n    # add Weeks from sub\n    _sub  = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\n    _sub['Patient'] = _sub['Patient_Week'].apply(lambda x:x.split('_')[0])\n    _sub['Weeks']   = _sub['Patient_Week'].apply(lambda x:x.split('_')[1]).astype('int')\n    _sub.drop(columns=['FVC','Confidence'],inplace=True)\n\n    # update test\n    test = _sub.merge(test,on='Patient')\n    test['rWeeks'] = test['Weeks']-test['Week0']\n    \n    # make prediction\n    test['predFVC']=0\n    \n    \n    #\n    # predict FVC and Confidence\n    #\n    \n    test_orig = test.copy()\n    testModelGroup = [ test.copy() for _ in range(len(CNF['prefixes'])) ]\n\n    CNF['jpgPath']='./segmented'\n    for prefixsModelGroup,weightsModelGroup,test in zip(CNF['prefixes'],CNF['ensembleWeight'],testModelGroup):\n        \n        for prefix,ensemble_weight in zip(prefixsModelGroup,weightsModelGroup):\n\n            for megaFold in range(nMegaFold):\n\n                # fitterをtest用に設定\n                if 'LGBM' in prefix:\n                    if 'week' in prefix:\n                        fvc_fitter = FVCFitterWeekLGBM(CNF,fold_id=megaFold,prefix=prefix)\n                    else:\n                        fvc_fitter = FVCFitterLGBM(CNF,fold_id=megaFold,prefix=prefix)\n            \n                else:\n                    fvc_fitter = FVCFitterNN01(CNF,fold_id=megaFold,prefix=prefix)\n\n                fvc_fitter.CNF['mskDfPath'] = './LTRC_feat_test_df.csv'\n\n                test['predFVC'] += ensemble_weight * ( fvc_fitter.predict( test.drop(columns='FVC'),image_path='test' ) )/5\n                \n                del fvc_fitter\n                gc.collect()\n\n    \n    for cnt,(prefixsModelGroup,weightsModelGroup,test) in enumerate( zip(CNF['prefixes'],CNF['ensembleWeight'],testModelGroup) ):\n        \n        confidence_fitter_prefix=f\"lastconfModelGroup{cnt}\"\n        \n        test['predConfidence']=0\n        for megaFold in range(nMegaFold):\n\n            conf_fitterRidge = ConfidenceFitter(CNF,fold_id=fold,model_name='ridge',prefix=confidence_fitter_prefix+'RIDGE')\n            conf_fitterSVR   = ConfidenceFitter(CNF,fold_id=fold,model_name='svr',prefix=confidence_fitter_prefix+'SVR')\n\n            _predRidge = conf_fitterRidge.predict(test[['Patient']+conf_cols])\n            _predSVR   = conf_fitterSVR.predict(test[['Patient']+conf_cols])\n\n            _pred = np.sqrt( 0.2*(_predSVR**2) + 0.8*(_predRidge**2) )\n\n            test['predConfidence'] += (_pred**2)/5\n\n        test['predConfidence'] = np.sqrt(test['predConfidence']) \n        \n    # モデルグループをアンサンブル\n    \n    _predFVC = np.mean([test.predFVC for test in testModelGroup],0)\n    _predConfidence = np.sqrt( np.mean([test.predConfidence**2 for test in testModelGroup],0) )\n    \n    # ここで繰り返し変数的に使っていたtestを元に戻す\n    \n    test = test_orig\n    \n    test['predFVC'] = _predFVC\n    test['predConfidence'] = _predConfidence \n    \n        \n    #\n    # save submission file \n    #\n    \n    sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\n\n    for k in sub.Patient_Week.values:\n\n        p, w = k.split('_')\n        w = int(w) \n\n        sub.loc[sub.Patient_Week == k, 'FVC'] = test.loc[ (test.Patient==p)&(test.Weeks==w), 'predFVC'].values\n        sub.loc[sub.Patient_Week == k, 'Confidence'] = test.loc[ (test.Patient==p)&(test.Weeks==w), 'predConfidence'].values\n\n    sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)\n    \n    # check submission file\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    subPatientList = sub.Patient.unique().tolist()\n\n    if len(subPatientList)==5:\n\n        for p in subPatientList:\n\n            plt.scatter( sub.loc[sub.Patient==p,'Weeks'], sub.loc[sub.Patient==p,'FVC'])\n            plt.show()\n\n    \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!rm -rf test_verticalX\n!rm -rf test_verticalY","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for cnt,prefixsModelGroup in enumerate( CNF['prefixes'] ):\n    \n    for prefix in prefixsModelGroup:\n\n        if 'LGBM' in prefix:\n\n            !rm -rf {prefix}\n        \n!rm -rf ./segmented\n\nfor cnt,prefixsModelGroup in enumerate( CNF['prefixes'] ):\n        \n    confidence_fitter_prefix=f\"lastconfModelGroup{cnt}\"\n    \n    !rm -rf {confidence_fitter_prefix+'RIDGE'}\n    !rm -rf {confidence_fitter_prefix+'SVR'}","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}