{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 分類モデルの作成  \n入ってきた画像が、そもそも糸球体(セグメンテーション対象)を含むのかどうかを分類するモデルを作成する。  \n過去コンペの解法の多くは、一度この分類器を通すことで大幅に精度が向上していた。  \n恐らくこちらのモデルにはそこまで精度が要求されないので、軽めのモデル(enetb0-b1)で良いのだと思う。  \n何を正解ラベルとするかだが、mask画像のsumをとったときに、閾値を設けて0or1で分ければ良いと思う。(閾値は色々試す)  \n\n* v1_2 :  \n    confusion matrixの分析コードを追加  \n    enet-b1に変更  \n    4fold学習させる  "},{"metadata":{"trusted":true},"cell_type":"code","source":"DEBUG = False","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## import"},{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\npackage_path = '../input/efficientnet-pytorch/EfficientNet-PyTorch/EfficientNet-PyTorch-master'\nsys.path.append(package_path)\n\nimport os\nimport pandas as pd\nimport numpy as np\nimport random\nimport time\nfrom tqdm import tqdm\nimport json\nimport cv2\nfrom PIL import Image\nimport tifffile\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import accuracy_score, roc_auc_score\nimport torch\nfrom torch import nn, optim\nimport torchvision\nfrom torchvision import transforms\nfrom torch.utils.data import Dataset,DataLoader\nfrom torch.utils.data.sampler import SequentialSampler, RandomSampler\nfrom torch.optim import Adam\nfrom torch.optim.lr_scheduler import CosineAnnealingLR\nfrom scipy.ndimage.interpolation import zoom\nimport albumentations as A\nfrom torch.nn import functional as F\nfrom albumentations.pytorch import ToTensorV2\nfrom efficientnet_pytorch import model as enet","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class CFG:\n    model_name = \"HuBMAP_tachyon_classifier_v1_2\"\n    backbone = 'efficientnet-b1'\n    Progress_Bar = False\n    max_grad_norm = 1000\n    gradient_accumulation_steps = 1\n    init_lr = 1e-3\n    weight_decay = 1e-5\n    image_size = 256\n    batch_size = 64\n    n_epochs = 5 if DEBUG else 25\n    n_fold = 4\n    train_fold = 4 #n_foldあるうち、この数だけしか学習しない\n    n_seed = 1\n    seed = 42\n    num_workers = 4\n    model_save_path = False\n    mask_th = 100 #256x256pixelのうち、これを超えるマスク部分があれば正解ラベルを1とする","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(device)\n\ndef seed_everything(seed=CFG.seed):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    \nseed_everything()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"BASE_PATH = \"../input/hubmap-kidney-segmentation/\"\nTRAIN_PATH = os.path.join(BASE_PATH, \"train\")\n\ntrain_df = pd.read_csv(os.path.join(BASE_PATH, \"train.csv\"))\nhubmap_df = pd.read_csv(os.path.join(BASE_PATH,\"HuBMAP-20-dataset_information.csv\"))\n\n#iafoss dataset\nTRAIN = '../input/hubmap-256x256/train/'\nMASKS = '../input/hubmap-256x256/masks/'\nLABELS = '../input/hubmap-kidney-segmentation/train.csv'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/iafoss/256x256-images\nmean = np.array([0.65459856,0.48386562,0.69428385])\nstd = np.array([0.15167958,0.23584107,0.13146145])\n\ndef img2tensor(img,dtype:np.dtype=np.float32):\n    if img.ndim==2 : img = np.expand_dims(img,2)\n    img = np.transpose(img,(2,0,1))\n    return torch.from_numpy(img.astype(dtype, copy=False))\n\nclass HuBMAPDataset(Dataset):\n    def __init__(self, fold, train=True, tfms=None):\n        ids = pd.read_csv(LABELS).id.values\n        kf = KFold(n_splits=CFG.n_fold,random_state=CFG.seed,shuffle=True)\n        ids = set(ids[list(kf.split(ids))[fold][0 if train else 1]])\n        if DEBUG:\n            self.fnames = [fname for fname in os.listdir(TRAIN) if fname.split('_')[0] in ids][:200]\n        else:\n            self.fnames = [fname for fname in os.listdir(TRAIN) if fname.split('_')[0] in ids]\n        self.train = train\n        self.tfms = tfms\n        \n    def __len__(self):\n        return len(self.fnames)\n    \n    def __getitem__(self, idx):\n        fname = self.fnames[idx]\n        img = cv2.cvtColor(cv2.imread(os.path.join(TRAIN,fname)), cv2.COLOR_BGR2RGB)\n        mask = cv2.imread(os.path.join(MASKS,fname),cv2.IMREAD_GRAYSCALE)\n        if self.tfms is not None:\n            augmented = self.tfms(image=img,mask=mask)\n            img,mask = augmented['image'],augmented['mask']\n        mask = 1 if mask.sum() > CFG.mask_th else 0\n        return img2tensor((img/255.0 - mean)/std),mask #img2tensor(mask)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Augmentation"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_aug(p=1.0):\n    return A.Compose([\n        A.HorizontalFlip(),\n        A.VerticalFlip(),\n        A.RandomRotate90(),\n        A.ShiftScaleRotate(shift_limit=0.0625, scale_limit=0.2, rotate_limit=15, p=0.9, \n                         border_mode=cv2.BORDER_REFLECT),\n        A.OneOf([\n            A.OpticalDistortion(p=0.3),\n            A.GridDistortion(p=.1),\n            A.IAAPiecewiseAffine(p=0.3),\n        ], p=0.3),\n        A.OneOf([\n            A.HueSaturationValue(10,15,10),\n            A.CLAHE(clip_limit=2),\n            A.RandomBrightnessContrast(),            \n        ], p=0.3),\n    ], p=p)\n\ndef get_aug_lastnepo(p=1.0):\n    return A.Compose([\n        A.HorizontalFlip(),\n        A.VerticalFlip(),\n        A.RandomRotate90(),\n    ], p=p)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## model"},{"metadata":{"trusted":true},"cell_type":"code","source":"pretrained_model = {\n        'efficientnet-b0': '../input/efficientnet-pytorch/efficientnet-b0-08094119.pth',\n        'efficientnet-b1': '../input/efficientnet-pytorch/efficientnet-b1-dbc7070a.pth',\n        'efficientnet-b2': '../input/efficientnet-pytorch/efficientnet-b2-27687264.pth',\n        'efficientnet-b3': '../input/efficientnet-pytorch/efficientnet-b3-c8376fa2.pth',\n        'efficientnet-b4': '../input/efficientnet-pytorch/efficientnet-b4-e116e8b3.pth',\n        'efficientnet-b5': '../input/efficientnet-pytorch/efficientnet-b5-586e6cc6.pth',\n        \n    }\n\n\nclass enetv2(nn.Module):\n    def __init__(self, backbone, out_dim=1):\n        super(enetv2, self).__init__()\n        self.enet = enet.EfficientNet.from_name(backbone)\n        self.enet.load_state_dict(torch.load(pretrained_model[backbone]))\n\n        self.myfc = nn.Linear(self.enet._fc.in_features, out_dim)\n        self.enet._fc = nn.Identity()\n        self.sigmoid = nn.Sigmoid()\n\n    def extract(self, x):\n        return self.enet(x)\n\n    def forward(self, x):\n        x = self.extract(x)\n        x = self.myfc(x)\n        #x = self.sigmoid(x)\n        return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## train valid func"},{"metadata":{"trusted":true},"cell_type":"code","source":"def train(model, iterator, optimizer, criterion, device, freeze):\n    \n    epoch_loss = 0\n    model.train()\n    \n    #プログレスバーを表示するか否か\n    bar = tqdm(iterator) if CFG.Progress_Bar else iterator\n    \n    for (x, y) in bar:\n        x = torch.tensor(x, device=device, dtype=torch.float32)\n        y = torch.tensor(y, device=device, dtype=torch.float32)\n        optimizer.zero_grad()\n        y_pred = model(x)\n        loss = criterion(y_pred, y.unsqueeze(1))\n        loss.backward()\n        optimizer.step()\n        loss_np = loss.detach().cpu().numpy()\n        epoch_loss += loss_np\n        \n        if CFG.Progress_Bar:\n            bar.set_description('Training loss: %.5f' % (loss_np))\n        \n    return epoch_loss/len(iterator)\n\ndef evaluate(model, iterator, criterion, device):\n    \n    epoch_loss = 0\n    preds = torch.tensor([])\n    targets = torch.tensor([])\n    model.eval()\n    \n    bar = tqdm(iterator) if CFG.Progress_Bar else iterator\n    \n    with torch.no_grad():\n        for (x, y) in bar:\n            x = torch.tensor(x, device=device, dtype=torch.float32)\n            y = torch.tensor(y, device=device, dtype=torch.float32)\n            \n            y_pred = model(x)\n            loss = criterion(y_pred, y.type_as(y_pred).unsqueeze(1))\n            loss_np = loss.detach().cpu().numpy()\n            epoch_loss += loss_np\n            ###logitです！\n            y_pred = torch.sigmoid(y_pred)\n            ###\n            preds = torch.cat([preds,y_pred.detach().cpu()],dim=0)\n            targets = torch.cat([targets,y.detach().cpu()],dim=0)\n            \n            if CFG.Progress_Bar:\n                bar.set_description('Validation loss: %.5f' % (loss_np))\n    \n    try:\n       val_roc = roc_auc_score(targets, preds)\n    except ValueError:\n       val_roc = -1\n    \n    return epoch_loss/len(iterator), val_roc\n\ndef fit_model(model, name, train_iterator, valid_iterator, optimizer, loss_criterion, device, freeze, epochs):\n    \"\"\" Fits a dataset to model\"\"\"\n    best_valid_score = float('inf')\n    \n    train_losses = []\n    valid_losses = []\n    valid_rocs = []\n    \n    for epoch in range(epochs):\n        #最後の3世代、augmentationなしで学習\n        if (freeze==False) and (epoch == CFG.n_epochs-11): #最後5世代に対して、Augを適用しない\n            print(\"-_-_-_-_-_-_-_-_-\")\n            print(\"No augment mode\")\n            print(\"-_-_-_-_-_-_-_-_-\")\n            #データセット\n            train_data = HuBMAPDataset(fold=fold,train=True,tfms=get_aug_lastnepo())\n            valid_data = HuBMAPDataset(fold=fold,train=False)\n\n            #データローダー\n            #ローカル変数に上書き(あまり清潔な書き方ではない…)\n            train_iterator = DataLoader(train_data,shuffle=True,batch_size=CFG.batch_size,num_workers=CFG.num_workers)\n            valid_iterator = DataLoader(valid_data,shuffle=False,batch_size=8,num_workers=CFG.num_workers)\n        \n        if scheduler:\n            scheduler.step(epoch)\n        start_time = time.time()\n    \n        train_loss = train(model,train_iterator,optimizer,loss_criterion,device,freeze)\n        valid_loss, valid_roc = evaluate(model,valid_iterator,loss_criterion,device)\n        \n        train_losses.append(train_loss)\n        valid_losses.append(valid_loss)\n        valid_rocs.append(valid_roc)\n\n        \"\"\"if valid_dice < best_valid_score:\n            best_valid_score = valid_dice\n            if CFG.model_save_path:\n                torch.save(model.state_dict(), os.path.join(model_save_path,f'{name}.pt'))\n            else:\n                torch.save(model.state_dict(), f'{name}_best.pt')\"\"\"\n        \n        end_time = time.time()\n\n        epoch_mins, epoch_secs = (end_time-start_time)//60,round((end_time-start_time)%60)\n    \n        print(f'Epoch: {epoch+1:02} | Epoch Time: {epoch_mins:.0f}m {epoch_secs}s')\n        print(f'lr:{optimizer.param_groups[0][\"lr\"]:.7f}')\n        print(f'Train Loss: {train_loss:.3f}')\n        print(f'Val. Loss: {valid_loss:.3f} | Val. ROC Score: {valid_roc:.3f}')\n        \n        #最後のAugの効果を見たいため、最終世代のみを出力\n        if not freeze:\n            torch.save(model.state_dict(), f'{name}_final.pt')\n        \n    return train_losses, valid_losses, valid_rocs","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## run training"},{"metadata":{"trusted":true},"cell_type":"code","source":"tr_loss=[]\nval_loss=[]\nval_roc=[]\nmodels = []\n\n\nfor fold in range((1 if DEBUG else CFG.train_fold)):\n    print(f\"Fitting on Fold {fold+1}\")\n    \n    #データセット\n    train_data = HuBMAPDataset(fold=fold,train=True,tfms=get_aug())\n    valid_data = HuBMAPDataset(fold=fold,train=False)\n    \n    #データローダー\n    train_iterator = DataLoader(train_data,shuffle=True,batch_size=CFG.batch_size,num_workers=CFG.num_workers)\n    valid_iterator = DataLoader(valid_data,shuffle=False,batch_size=8,num_workers=CFG.num_workers)\n    \n    #モデルの呼び出し(設計図からインスタンスへ)\n    model = enetv2(CFG.backbone).to(device)\n    name = CFG.model_name + \"_f\" + str(fold)\n    \n    #最初の3世代を出力層以外freezeして学習\n    print(\"+-+-+-+-+-+-+-+-+\")\n    print(\"pretrain mode\")\n    print(\"+-+-+-+-+-+-+-+-+\")\n    loss_criterion = nn.BCEWithLogitsLoss()\n    opt= Adam(model.parameters(), lr=1e-3)\n    scheduler=None\n    \n    head_name = [\"myfc.weight\",\"myfc.bias\"]\n    for hname,param in model.named_parameters():\n        if hname in head_name:\n            param.requires_grad = True\n        else:\n            param.requires_grad = False\n    \n    nouse0,nouse1,nouse2 = fit_model(model, name, train_iterator, valid_iterator, opt, loss_criterion, device,freeze=True,epochs=5)\n    \n    print(\"@*@*@*@*@*@*@*@*@\")\n    print(\"fulltrain mode\")\n    print(\"@*@*@*@*@*@*@*@*@\")\n    #損失関数の定義\n    loss_criterion = nn.BCEWithLogitsLoss()\n    \n    #最適化手法の定義\n    opt= Adam(model.parameters(), lr=CFG.init_lr)\n    \n    #スケジューラーの定義\n    scheduler = CosineAnnealingLR(opt,CFG.n_epochs-5)\n    \n    for hname,param in model.named_parameters():\n        param.requires_grad = True\n    #全ての情報をfit_modelに入れて、学習を開始します\n    temp_tr_loss, temp_val_loss, temp_val_rocs = fit_model(model, name, train_iterator, valid_iterator, opt, loss_criterion, device,freeze=False, epochs=CFG.n_epochs-5)\n    \n    \n    #lossと評価指標に対するスコアを記録します\n    tr_loss.append(temp_tr_loss)\n    val_loss.append(temp_val_loss)\n    val_roc.append(temp_val_rocs)\n    \n    #foldごとにモデルを定義する為、学習し終わったモデルはリストに保持しておきます\n    models.append(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(len(tr_loss)):\n    fig,ax = plt.subplots(nrows=1, ncols=2, figsize=(20,5))\n    ax[0].plot(tr_loss[i])\n    ax[0].set_title('Training and Validation Loss')\n    ax[0].plot(val_loss[i])\n    ax[0].set_xlabel('Epoch')\n\n    ax[1].plot(val_roc[i])\n    ax[1].set_title('Val roc Score')\n    ax[1].set_xlabel('Epoch')\n\n\n    ax[0].legend();\n    ax[1].legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# confusion matrixの作成"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import classification_report","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#(model数*TTA数)回すので注意\ndef get_predictions(model, iterator):\n    model.eval()\n    bar = tqdm(iterator) if CFG.Progress_Bar else iterator\n    \n    with torch.no_grad():\n        res = np.array([])\n        ans = np.array([])\n        for x,y in bar:\n            x = torch.tensor(x, device=device, dtype=torch.float32)\n            y = torch.tensor(y, device=device, dtype=torch.float32)\n            y_pred = model(x)\n            y_pred = torch.sigmoid(y_pred)\n            res = np.append(res, y_pred.detach().cpu().numpy())\n            ans = np.append(ans, y.detach().cpu().numpy())\n    return res,ans","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred = np.array([])\ny = np.array([])\nfor fold in range((1 if DEBUG else CFG.train_fold)):\n    print(f\"Validating on Fold {fold+1}\")\n    \n    #データセット\n    valid_data = HuBMAPDataset(fold=fold,train=False)\n    #データローダー\n    valid_iterator = DataLoader(valid_data,shuffle=False,batch_size=8,num_workers=CFG.num_workers)\n    \n    model = models[fold]\n    pred_res,y_res = get_predictions(model,valid_iterator)\n    \n    pred = np.append(pred, pred_res)\n    y = np.append(y, y_res)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import confusion_matrix, roc_curve, auc\n\ndef youden_index(y_pred,y):\n    fpr, tpr, thres = roc_curve(y, y_pred)\n    auc_ = auc(fpr, tpr)\n    # 特異度\n    sng = 1 - fpr\n    # Youden indexを用いたカットオフ基準\n    cutoff_criterion = tpr + sng - 1\n    #print(f'{model_name}, auc:{auc}')\n    return cutoff_criterion,tpr,fpr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#youden indexを用いた最適閾値の取得\n#閾値の算出\ncutoff,tpr,fpr = youden_index(pred,y) #cutoff:y軸,fpr:x軸 \ncutoff_opt = [max(cutoff),fpr[cutoff.argmax()],tpr[cutoff.argmax()]] #arxmax:最大値のindexを取得\n#これが最適な閾値\nprint(\"[Youden_index,x軸,y軸]\")\nprint(cutoff_opt)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_binary = np.array(pred>cutoff_opt[1],dtype=\"int8\")\nprint(classification_report(pred_binary,y))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import metrics\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfpr, tpr, thres = roc_curve(y,pred)\nauc_ = auc(fpr, tpr)\n\n# ROC曲線をプロット\nplt.plot(fpr, tpr, label='ROC curve (area = %.2f)'%auc_)\nplt.scatter(cutoff_opt[1],cutoff_opt[2])\nplt.legend()\nplt.title('ROC curve')\nplt.xlabel('False Positive Rate')\nplt.ylabel('True Positive Rate')\nplt.text(cutoff_opt[1]+0.1, cutoff_opt[2]-0.1, f'cutoff = {cutoff_opt[0]:.3f}', fontsize=12);\nplt.grid(True)","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}