{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction\n![CNN Schema](https://raw.githubusercontent.com/jxtrbtk/kaggle/master/aptos2019-blindness-detection/CNNSchema.png)\n\nThis kernel is a playground around the following idea. CNN are made of two part, convultions that acts like a features extractor and a fully connected neural network used for classification. What if we replace the neural network classifier by an XGBoost based one ? XGBoost is most often very efficient and makes it easy to achieve very good performances without a lot of tuning effort.\n\nAt the very beginning, I had better results with XGB, but the more optimized the process, the more similar the performance were...\nHope you'll enjoy the journey !"},{"metadata":{},"cell_type":"markdown","source":"# Imports, settings and references"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np \nimport time\n\nimport cv2\nfrom PIL import Image\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset\nfrom torch.utils.data.sampler import SubsetRandomSampler\n\nimport torchvision\nfrom torchvision import transforms\n\nimport xgboost as xgb\nfrom sklearn.metrics import cohen_kappa_score\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.preprocessing import quantile_transform\n\nimport pickle\n\nDEVICE = torch.device(\"cuda:0\")\nDATA_SOURCE = os.path.join(\"..\",\"input\",\"aptos2019-blindness-detection\")\nMODEL_SOURCE = os.path.join(\"..\",\"input\",\"torchvisionmodelspartial1\")\nMODEL_SIZE = 224","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Pre-processing"},{"metadata":{"trusted":true},"cell_type":"code","source":"def crop_image(img,tol=7):\n    w, h = img.shape[1],img.shape[0]\n    gray_img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)\n    gray_img = cv2.blur(gray_img,(5,5))\n    shape = gray_img.shape \n    gray_img = gray_img.reshape(-1,1)\n    quant = quantile_transform(gray_img, n_quantiles=256, random_state=0, copy=True)\n    quant = (quant*256).astype(int)\n    gray_img = quant.reshape(shape)\n    xp = (gray_img.mean(axis=0)>tol)\n    yp = (gray_img.mean(axis=1)>tol)\n    x1, x2 = np.argmax(xp), w-np.argmax(np.flip(xp))\n    y1, y2 = np.argmax(yp), h-np.argmax(np.flip(yp))\n    if x1 >= x2 or y1 >= y2 : # something wrong with the crop\n        return img # return original image\n    else:\n        img1=img[y1:y2,x1:x2,0]\n        img2=img[y1:y2,x1:x2,1]\n        img3=img[y1:y2,x1:x2,2]\n        img = np.stack([img1,img2,img3],axis=-1)\n    return img\n\ndef process_image(image, size=512):\n    image = cv2.resize(image, (size,int(size*image.shape[0]/image.shape[1])))\n    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    try:\n        image = crop_image(image, tol=15)\n    except Exception as e:\n        image = image\n        print( str(e) )\n    return image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"*inspired by: https://www.kaggle.com/ratthachat/aptos-updated-preprocessing-ben-s-cropping*"},{"metadata":{},"cell_type":"markdown","source":"# PyTorch's style data loader defintion"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"class RetinopathyDataset(Dataset):\n\n    def __init__(self, transform, is_test=False):\n        self.transform = transform\n        self.base_transform = transforms.Resize((MODEL_SIZE, MODEL_SIZE))\n        self.is_test = is_test \n        if not os.path.exists(\"cache\"): os.mkdir(\"cache\")\n        if is_test : file = \"test.csv\"\n        else : file = \"train.csv\"\n        csv_file = os.path.join(DATA_SOURCE, file)\n        df = pd.read_csv(csv_file)\n        self.data = df.reset_index(drop=True)\n            \n    def __len__(self):\n        return len(self.data)\n\n    def __getitem__(self, idx):\n        if self.is_test : archive = \"test_images\"\n        else : archive = \"train_images\"\n        folder = os.path.join(DATA_SOURCE, archive)\n        code = str(self.data.loc[idx, 'id_code'])\n        file = code + \".png\"\n        cache_path = os.path.join(\"cache\",code+\".png\")\n        cached = os.path.exists(cache_path)\n        if not cached : \n            path = os.path.join(folder, file)\n            image = cv2.imread(path)\n            image = process_image(image)\n            imgpil = Image.fromarray(image)\n            imgpil = self.base_transform(imgpil)\n            imgpil.save(cache_path,\"PNG\")\n        imgpil = Image.open(cache_path)\n        img_tensor = self.transform(imgpil)\n        if self.is_test : return {'image': img_tensor} \n        else : \n            label = self.data.loc[idx, \"diagnosis\"]\n            return {'image': img_tensor, 'label': label}\n        \n\n    def get_df(self):\n        return self.data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"*adapted from: https://www.kaggle.com/abhishek/very-simple-pytorch-training-0-59*"},{"metadata":{},"cell_type":"markdown","source":"# Re-train the pre-trained model"},{"metadata":{"trusted":true},"cell_type":"code","source":"# first we will prepare the dataset and create the folds\nNUM_FOLDS = 5\n\ndata_augmentation = transforms.Compose([\n    transforms.RandomRotation((-15, 15)),\n    transforms.Resize(224),\n    transforms.RandomHorizontalFlip(),\n    transforms.RandomVerticalFlip(),\n    transforms.ToTensor(),\n    transforms.Normalize([0.485, 0.456, 0.406], \n                         [0.229, 0.224, 0.225])\n])\n\nDATA = RetinopathyDataset(data_augmentation)\ndf = DATA.get_df()\nskf = StratifiedKFold(n_splits=NUM_FOLDS)\nfolds_generator = skf.split(df.index.values, df.diagnosis.values)\ndata_train, data_eval = [], [] \nfor t, e in folds_generator:\n    data_train.append(t)\n    data_eval.append(e)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_dataloader_for_fold(n, data, train_data, eval_data, batch_size):\n    \"\"\" return the train and eval dataloader for a fold\n    \"\"\"    \n    train_sampler = SubsetRandomSampler(train_data[n])\n    valid_sampler = SubsetRandomSampler(eval_data[n])\n\n    data_loader_train = torch.utils.data.DataLoader(data, \n                    batch_size=batch_size, drop_last=False, \n                    sampler=train_sampler)\n    data_loader_eval = torch.utils.data.DataLoader(data, \n                    batch_size=batch_size, drop_last=False, \n                    sampler=valid_sampler)\n    \n    return data_loader_train, data_loader_eval","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Classificator(nn.Module):\n    \"\"\" classifier layer used to retrain the CNN\n    \"\"\"    \n    def __init__(self, size=128):\n        super(Classificator, self).__init__()\n        self.size = size\n        self.network = nn.Sequential(\n              nn.BatchNorm1d(size),\n              nn.Dropout(p=0.3),\n              nn.Linear(in_features=size, out_features=5, bias=True),\n        )        \n    def forward(self, x):\n        ## Define forward behavior\n        return self.network(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_base_model():\n    \"\"\" get the pretrained model\n    \"\"\"    \n    model = torchvision.models.densenet161(pretrained=False)\n    model_path = os.path.join(MODEL_SOURCE, \"densenet161.pth\")\n    model.load_state_dict(torch.load(model_path))\n    in_features = model.classifier.in_features\n    model.classifier = Classificator(in_features)\n    model = model.to(DEVICE)\n    model.eval()\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_model(model, optimizer, scheduler, train_data_loader, eval_data_loader, \n                file_name, num_epochs = 50, patience = 7, prev_loss = 1000.00):\n    \"\"\" train the model\n    arguments : model, optimizer, scheduler, train_data_loader, eval_data_loader\n        file_name: name of the file to save the best model \n        num_epochs: maximum number of epochs\n        patience: number of epochs to wait if no improvements\n        prev_loss: previous loss achieved, to surpass to have the model saved\n    return: \n        best loss achieved (previous loss if not surpassed)\n    \"\"\"    \n    criterion = nn.CrossEntropyLoss()\n    countdown = patience\n    best_loss = 1000.00\n    since = time.time()\n    for epoch in range(num_epochs):\n        running_loss = 0.0\n        counter = 0\n        for bi, d in enumerate(train_data_loader):\n            inputs = d[\"image\"].to(DEVICE, dtype=torch.float)\n            labels = d[\"label\"].to(DEVICE, dtype=torch.long)\n            # batch norm layers needs more than 1 set of data\n            # this is to skip the last batch if it's only 1 image\n            if inputs.shape[0] > 1 :\n                counter += inputs.size(0)\n                model.to(DEVICE)\n                model.train()\n                optimizer.zero_grad()\n                outputs = model(inputs) \n                loss = criterion(outputs, labels)\n                loss.backward()\n                optimizer.step()\n                running_loss += loss.item() * inputs.size(0)\n                loss_val = running_loss / counter\n                print(\"{:7} {:.4f} {:.4f}\".format(counter, loss.item()*1, loss_val), end=\"\\r\")\n        epoch_loss = running_loss / ( len(train_data_loader) * train_data_loader.batch_size)\n        time_elapsed = time.time() - since\n        print(\" T{:3}/{:3} loss: {:.4f} ({:3.0f}m {:2.0f}s)\".format( \n            epoch, num_epochs - 1, epoch_loss,time_elapsed // 60, time_elapsed % 60))\n        running_loss = 0.0\n        counter = 0\n        for bi, d in enumerate(eval_data_loader):\n            inputs = d[\"image\"].to(DEVICE, dtype=torch.float)\n            counter += inputs.size(0)\n            labels = d[\"label\"].to(DEVICE, dtype=torch.long)\n            model.to(DEVICE)\n            model.eval()\n            with torch.no_grad():\n                outputs = model(inputs)\n            loss = criterion(outputs, labels)\n            running_loss += loss.item() * inputs.size(0)\n            loss_val = running_loss / counter\n            print(\"{:7} {:.4f} {:.4f}\".format(counter, loss.item()*1, loss_val), end=\"\\r\")\n        epoch_loss = running_loss / ( len(eval_data_loader) * eval_data_loader.batch_size)\n        if epoch_loss < best_loss : \n            best_loss = epoch_loss\n            if epoch_loss < prev_loss:\n                torch.save(model.state_dict(), file_name)\n                prev_loss = epoch_loss\n                print(\"*\", end=\"\")\n            else:\n                print(\".\", end=\"\")\n            countdown = patience\n        else:\n            print(\"{:1}\".format(countdown), end=\"\")\n            countdown -= 1\n        time_elapsed = time.time() - since\n        print(\"E{:3}/{:3} loss: {:.4f} ({:3.0f}m {:2.0f}s)\".format( \n            epoch, num_epochs - 1, epoch_loss,time_elapsed // 60, time_elapsed % 60 ))\n        scheduler.step() #epoch_loss\n\n        if countdown <= 0 : break\n\n    return prev_loss\n    print(\"done.\")\n# Model training","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train the model num_round_per_fold times for each fold\n# and the save the best model for each fold\nbatch_size = 56\nnum_round_per_fold = 3\nfor no in range(NUM_FOLDS):\n    print(\"-\"*22, \"fold\",no)\n    bst_loss = 10000.00\n    for r in range(num_round_per_fold):\n        print(\"-\"*11,\"round\",r)\n        data_loader_train, data_loader_eval = get_dataloader_for_fold(no, \n                                    DATA, data_train, data_eval, batch_size)\n        model = get_base_model()\n        plist = [{\"params\": model.features.denseblock3.parameters(), \"lr\":0.0001},\n                 {\"params\": model.features.denseblock4.parameters(), \"lr\":0.0001},\n                 {\"params\": model.classifier.parameters()}]\n        optimizer = optim.Adam(plist, lr=0.001, amsgrad=True)\n        scheduler = optim.lr_scheduler.MultiStepLR(optimizer, [1,5], gamma=0.1, last_epoch=-1)        \n        bst_loss = train_model(model, optimizer, scheduler, \n                               data_loader_train, data_loader_eval, \n                               \"tmp\"+str(no)+\".pth\", prev_loss=bst_loss, \n                               num_epochs=9, patience=3)\n    print(\"-\"*22, \"best loss\", bst_loss)\n    print(\"\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Extract train features from CNN"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_trained_model(no): \n    \"\"\" reload and return the retrained model for the given fold \n    \"\"\"    \n    extractor = torchvision.models.densenet161(pretrained=False)\n    in_features = extractor.classifier.in_features\n    extractor.classifier = Classificator(in_features)\n    model_path = os.path.join(\"tmp\"+str(no)+\".pth\") #no\n    extractor.load_state_dict(torch.load(model_path))\n    extractor = extractor.to(DEVICE)\n    extractor.eval()\n    return extractor","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_extractor_model(no):\n    \"\"\" reload and return the retrained model for the given fold\n        and make last layer identity\n    \"\"\"    \n    extractor = get_trained_model(no)\n    extractor.classifier = nn.Identity()\n    extractor = extractor.to(DEVICE)\n    extractor.eval()\n    return extractor","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_train_features(data_loader, extractor):\n    \"\"\" return 2 arrays of features extracted, and targets  \n    \"\"\"    \n    for bi, d in enumerate(data_loader):\n        print(\".\", end=\"\")\n        img_tensor = d[\"image\"].to(DEVICE)\n        target = d[\"label\"].numpy()\n        with torch.no_grad(): feature = extractor(img_tensor)\n        feature = feature.cpu().detach().squeeze(0).numpy()\n        if bi == 0 :\n            features = feature \n            targets = target \n        else :\n            features = np.concatenate([features, feature], axis=0)\n            targets = np.concatenate([targets, target], axis=0)\n\n    return features, targets","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fit the XGBoost model"},{"metadata":{"trusted":true},"cell_type":"code","source":"XGBOOST_PARAM = {\n    \"random_state\"      : 42,\n    \"n_estimators\"      : 200,\n    \"objective\"         : \"multi:softmax\",\n    \"num_class\"         : 5,\n    \"eval_metric\"       : \"mlogloss\",\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# for each fold, get the data loader, extractor, \n# extract the features (loaded and process each time, but we can have data aug.) \n# calcul the weights table, create and fit the XGB model\nbatch_size = 64\neval_set = []\nfor no in range(NUM_FOLDS):\n    print(\"-\"*22, \"fold\",no)\n    data_loader_train, data_loader_eval = get_dataloader_for_fold(no, \n                                DATA, data_train, data_eval, batch_size)\n    extractor = get_extractor_model(no)\n\n    print(\"...........|.............................................|\")\n    features_eval, targets_eval = get_train_features(data_loader_eval,\n                                                     extractor)\n    features_train, targets_train = get_train_features(data_loader_train,\n                                                       extractor)\n    print(\"\")\n\n    xgb_model = xgb.XGBClassifier(**XGBOOST_PARAM)\n    xgb_model = xgb_model.fit(features_train,targets_train.reshape(-1),\n                        eval_set=[(features_eval, targets_eval.reshape(-1))],\n                        early_stopping_rounds=20,\n                        verbose=False)\n    print(\"score\",xgb_model.evals_result()[\"validation_0\"][\"mlogloss\"][-1])\n    pickle.dump(xgb_model, open(\"xgb_model_\"+str(no), \"wb\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Cross-Validation\n## .using the full CNN models"},{"metadata":{"trusted":true},"cell_type":"code","source":"# change the data augmentation of the dataset object\nbase_transform = transforms.Compose([\n    transforms.Resize(224),\n    transforms.ToTensor(),\n    transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])\n])\nDATA.transform = base_transform","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# first we make predictions using the CNN models for each fold \n# and calculate the score of all predictions\npredictions, targets = np.zeros(len(DATA)), np.zeros(len(DATA))\nbatch_slice = (0, 0)\nsoftmax = nn.Softmax(dim=1)\nfor no in range(NUM_FOLDS):\n    _, data_loader_eval = get_dataloader_for_fold(no, \n                DATA, data_train, data_eval, 64)\n    model = get_trained_model(no) #no\n    for bi, d in enumerate(data_loader_eval):\n        inputs = d[\"image\"].to(DEVICE, dtype=torch.float)\n        batch_slice = (batch_slice[1], batch_slice[1]+inputs.size(0))\n        with torch.no_grad():\n            outputs = model(inputs)\n            outputs = softmax(outputs)\n        predictions[batch_slice[0]:batch_slice[1]] = \\\n                outputs.cpu().detach().squeeze(0).numpy().argmax(axis=1)\n        targets[batch_slice[0]:batch_slice[1]] = d[\"label\"]\n        \nprint(\"Cohen Kappa quadratic score\", \n      cohen_kappa_score(targets, predictions, weights=\"quadratic\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## .using the XGB models"},{"metadata":{"trusted":true},"cell_type":"code","source":"# now let's do it with the XGB models \npredictions, targets = np.zeros(0), np.zeros(0)\nfor no in range(NUM_FOLDS):\n    _, data_loader_eval = get_dataloader_for_fold(no, \n                                DATA, data_train, data_eval, batch_size)\n    features_eval, targets_eval = get_train_features(data_loader_eval,\n                                                     extractor)\n    print(\"\")\n    xgb_model = xgb.XGBClassifier()\n    model_path = os.path.join(\"xgb_model_\"+str(no))\n    xgb_model = pickle.load(open(model_path, \"rb\"))\n    prediction = xgb_model.predict(features_eval)\n    predictions = np.concatenate([predictions, prediction], axis=0)\n    targets = np.concatenate([targets, targets_eval], axis=0)\n\nprint(\"Cohen Kappa quadratic score\", \n      cohen_kappa_score(targets, predictions, weights=\"quadratic\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# cleaning\nif os.path.exists(\"cache\"):\n    for e in os.listdir(\"cache\"):\n        os.remove(os.path.join(\"cache\", e))\n    os.rmdir(\"cache\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Extract test features using CNN model"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_augmentation = transforms.Compose([\n    transforms.Resize((MODEL_SIZE, MODEL_SIZE)),\n    transforms.RandomHorizontalFlip(),\n    transforms.RandomVerticalFlip(),\n    transforms.ToTensor(),\n    transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])\n])\n\ndata_test = RetinopathyDataset(data_augmentation, is_test=True)\ndata_loader = torch.utils.data.DataLoader(data_test, \n                            batch_size=16, shuffle=False, \n                            num_workers=0, drop_last=False)\n\ndef get_test_features(data_loader, extractor):\n    \"\"\" return an array of features extracted  \n    \"\"\"    \n    for bi, d in enumerate(data_loader):\n        if bi % 8 == 0 : print(\".\", end=\"\")\n        img_tensor = d[\"image\"].to(DEVICE)\n        with torch.no_grad(): feature = extractor(img_tensor)\n        feature = feature.cpu().detach().numpy() #.squeeze(0) for batch_size > 1\n        if bi == 0 :\n            features = feature \n        else :\n            features = np.concatenate([features, feature], axis=0)\n    return features","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prediction using XGB"},{"metadata":{"trusted":true},"cell_type":"code","source":"# adding prediction for each model \n# we can loop several time to perform data augmentation (tta) \n# (note: a bit risky as we cache the image, should be done bucket by bucket\n#  and clean after each bucket to avoid filling all the disk space)\nprint(\"................................ v\")\npredictions = np.zeros((len(data_test),5))\nfor tta in range(2):\n    print(\"............tta\"+str(tta)+\"................\")\n    for no in range(NUM_FOLDS):\n        extractor = get_extractor_model(no)\n        features = get_test_features(data_loader, extractor)\n        print(\"\",no)\n        xgb_model = xgb.XGBClassifier()\n        model_path = os.path.join(\"xgb_model_\"+str(no))\n        xgb_model = pickle.load(open(model_path, \"rb\"))\n        prediction = xgb_model.predict_proba(features)\n        predictions = predictions + prediction ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# voting\nprediction_final = predictions.argmax(axis=1)\ncsv_file = os.path.join(DATA_SOURCE, \"sample_submission.csv\")\ndf = pd.read_csv(csv_file)\ndf[\"diagnosis\"] = prediction_final\ndf.to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# cache cleaning\nif os.path.exists(\"cache\"):\n    for e in os.listdir(\"cache\"):\n        os.remove(os.path.join(\"cache\", e))\n    os.rmdir(\"cache\")","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}