main-fusion-append.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 51,
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import csv\n",
    "import time\n",
    "import glob\n",
    "import socket\n",
    "import logging\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from google.cloud import bigquery\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "import torch.nn.functional as F\n",
    "import torchvision\n",
    "import torchvision.transforms as transforms\n",
    "from torch.utils.data import Dataset, DataLoader\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split, KFold, StratifiedKFold, GridSearchCV\n",
    "from sklearn.feature_selection import RFE, VarianceThreshold\n",
    "from sklearn.metrics import average_precision_score, roc_auc_score, roc_curve, accuracy_score, confusion_matrix\n",
    "from functools import partial"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "outputs": [],
   "source": [
    "def load_excel_data(path, sheet=0):\n",
    "    filename = os.path.basename(path).strip()\n",
    "    if isinstance(sheet, str):\n",
    "        logging.info(f'Loading {filename}, Sheet: {sheet}...')\n",
    "    else:\n",
    "        logging.info('Loading ' + filename + '...')\n",
    "    df_data = pd.read_excel(path, sheet)\n",
    "    logging.info(\"Done loading.\")\n",
    "    return df_data\n",
    "\n",
    "\n",
    "def load_image_data(path, patients, limit=False):\n",
    "    data_images = {}\n",
    "    for root, dirs, files in os.walk(path):\n",
    "        dirs.sort(key=int)\n",
    "        dirs = list(map(int, dirs))\n",
    "        dirs = [patient for patient in dirs if patient in patients]\n",
    "        if limit:\n",
    "            dirs = dirs[:limit]\n",
    "        for d in dirs:\n",
    "            logging.info(f\"Loading Patient {d}...\")\n",
    "            np_filenames = glob.glob(f\"{os.path.join(root, f'{d}')}/*/*.npy\")\n",
    "            data_images[d] = [np.load(np_filenames[0]), np.load(np_filenames[1])]\n",
    "        break\n",
    "    return data_images, dirs\n",
    "\n",
    "\n",
    "def random_seed(seed_value, use_cuda):\n",
    "    np.random.seed(seed_value)  # set np random seed\n",
    "    torch.manual_seed(seed_value)  # set torch seed\n",
    "    random.seed(seed_value)  # set python random seed\n",
    "    if use_cuda:\n",
    "        torch.cuda.manual_seed(seed_value)\n",
    "        torch.cuda.manual_seed_all(seed_value)\n",
    "        # reproducibility\n",
    "        torch.use_deterministic_algorithms(True)\n",
    "        torch.backends.cudnn.benchmark = False"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "# Plot Training Curve\n",
    "def plot_training_curve(path):\n",
    "    \"\"\" Plots the training curve for a model run, given the csv files\n",
    "    containing the train/validation error/loss.\n",
    "\n",
    "    Args:\n",
    "        path: The base path of the csv files produced during training\n",
    "    \"\"\"\n",
    "    train_err = np.loadtxt(\"{}_train_err.csv\".format(path))\n",
    "    val_err = np.loadtxt(\"{}_val_err.csv\".format(path))\n",
    "    train_loss = np.loadtxt(\"{}_train_loss.csv\".format(path))\n",
    "    val_loss = np.loadtxt(\"{}_val_loss.csv\".format(path))\n",
    "    plt.title(\"Train vs Validation Error\")\n",
    "    n = len(train_err)  # number of epochs\n",
    "    plt.plot(range(1, n + 1), train_err, label=\"Train\")\n",
    "    plt.plot(range(1, n + 1), val_err, label=\"Validation\")\n",
    "    plt.xlabel(\"Epoch\")\n",
    "    plt.ylabel(\"Error\")\n",
    "    plt.legend(loc='best')\n",
    "    plt.show()\n",
    "    plt.title(\"Train vs Validation Loss\")\n",
    "    plt.plot(range(1, n + 1), train_loss, label=\"Train\")\n",
    "    plt.plot(range(1, n + 1), val_loss, label=\"Validation\")\n",
    "    plt.xlabel(\"Epoch\")\n",
    "    plt.ylabel(\"Loss\")\n",
    "    plt.legend(loc='best')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_roc(labels, preds):\n",
    "    fpr, tpr, _ = roc_curve(labels, preds)\n",
    "    auc = roc_auc_score(labels, preds)\n",
    "\n",
    "    plt.plot(fpr, tpr, label=\"data 1, auc=\" + str(auc))\n",
    "    plt.legend(loc=4)\n",
    "    plt.show()\n",
    "    return\n",
    "\n",
    "###############################################################################"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "outputs": [],
   "source": [
    "# Model Classes\n",
    "class CNNDataset(Dataset):\n",
    "    def __init__(self, data, patient_ids):\n",
    "        self.data = data\n",
    "        self.patient_ids = patient_ids\n",
    "\n",
    "    def __len__(self):\n",
    "        return len(self.patient_ids)\n",
    "\n",
    "    def __getitem__(self, idx):\n",
    "        return self.data[self.patient_ids[idx]][\"input\"], self.data[self.patient_ids[idx]][\"label\"], self.data[self.patient_ids[idx]][\"features\"]\n",
    "\n",
    "\n",
    "def conv3x3x3(in_planes, out_planes, stride=1):\n",
    "    return nn.Conv3d(in_planes,\n",
    "                     out_planes,\n",
    "                     kernel_size=3,\n",
    "                     stride=stride,\n",
    "                     padding=1,\n",
    "                     bias=False)\n",
    "\n",
    "\n",
    "def conv1x1x1(in_planes, out_planes, stride=1):\n",
    "    return nn.Conv3d(in_planes,\n",
    "                     out_planes,\n",
    "                     kernel_size=1,\n",
    "                     stride=stride,\n",
    "                     bias=False)\n",
    "\n",
    "\n",
    "class BasicBlock(nn.Module):\n",
    "    expansion = 1\n",
    "\n",
    "    def __init__(self, in_planes, planes, stride=1, downsample=None):\n",
    "        super().__init__()\n",
    "\n",
    "        self.name = \"BasicBlock\"\n",
    "        self.conv1 = conv3x3x3(in_planes, planes, stride)\n",
    "        self.bn1 = nn.BatchNorm3d(planes)\n",
    "        self.relu = nn.ReLU(inplace=True)\n",
    "        self.conv2 = conv3x3x3(planes, planes)\n",
    "        self.bn2 = nn.BatchNorm3d(planes)\n",
    "        self.downsample = downsample\n",
    "        self.stride = stride\n",
    "\n",
    "    def forward(self, x):\n",
    "        residual = x\n",
    "\n",
    "        out = self.conv1(x)\n",
    "        out = self.bn1(out)\n",
    "        out = self.relu(out)\n",
    "\n",
    "        out = self.conv2(out)\n",
    "        out = self.bn2(out)\n",
    "\n",
    "        if self.downsample is not None:\n",
    "            residual = self.downsample(x)\n",
    "\n",
    "        out += residual\n",
    "        out = self.relu(out)\n",
    "\n",
    "        return out\n",
    "\n",
    "\n",
    "class Bottleneck(nn.Module):\n",
    "    expansion = 4\n",
    "\n",
    "    def __init__(self, in_planes, planes, stride=1, downsample=None):\n",
    "        super().__init__()\n",
    "\n",
    "        self.name = \"Bottleneck\"\n",
    "        self.conv1 = conv1x1x1(in_planes, planes)\n",
    "        self.bn1 = nn.BatchNorm3d(planes)\n",
    "        self.conv2 = conv3x3x3(planes, planes, stride)\n",
    "        self.bn2 = nn.BatchNorm3d(planes)\n",
    "        self.conv3 = conv1x1x1(planes, planes * self.expansion)\n",
    "        self.bn3 = nn.BatchNorm3d(planes * self.expansion)\n",
    "        self.relu = nn.ReLU(inplace=True)\n",
    "        self.downsample = downsample\n",
    "        self.stride = stride\n",
    "\n",
    "    def forward(self, x):\n",
    "        residual = x\n",
    "\n",
    "        out = self.conv1(x)\n",
    "        out = self.bn1(out)\n",
    "        out = self.relu(out)\n",
    "\n",
    "        out = self.conv2(out)\n",
    "        out = self.bn2(out)\n",
    "        out = self.relu(out)\n",
    "\n",
    "        out = self.conv3(out)\n",
    "        out = self.bn3(out)\n",
    "\n",
    "        if self.downsample is not None:\n",
    "            residual = self.downsample(x)\n",
    "\n",
    "        out += residual\n",
    "        out = self.relu(out)\n",
    "\n",
    "        return out\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "outputs": [],
   "source": [
    "class ResNet(nn.Module):\n",
    "\n",
    "    def __init__(self,\n",
    "                 block,\n",
    "                 layers,\n",
    "                 block_inplanes,\n",
    "                 model_depth,\n",
    "                 n_input_channels=3,\n",
    "                 conv1_t_size=7,\n",
    "                 conv1_t_stride=1,\n",
    "                 no_max_pool=False,\n",
    "                 shortcut_type='B',\n",
    "                 widen_factor=1.0,\n",
    "                 n_classes=400,\n",
    "                 num_additional_features=851):\n",
    "        super().__init__()\n",
    "\n",
    "        block_inplanes = [int(x * widen_factor) for x in block_inplanes]\n",
    "\n",
    "        self.name = f\"ResNet_pLGG_Classifer_depth{model_depth}\"\n",
    "\n",
    "        self.in_planes = block_inplanes[0]\n",
    "        self.no_max_pool = no_max_pool\n",
    "\n",
    "        self.conv1 = nn.Conv3d(n_input_channels,\n",
    "                               self.in_planes,\n",
    "                               kernel_size=(conv1_t_size, 7, 7),\n",
    "                               stride=(conv1_t_stride, 2, 2),\n",
    "                               padding=(conv1_t_size // 2, 3, 3),\n",
    "                               bias=False)\n",
    "        self.bn1 = nn.BatchNorm3d(self.in_planes)\n",
    "        self.relu = nn.ReLU(inplace=True)\n",
    "        self.maxpool = nn.MaxPool3d(kernel_size=3, stride=2, padding=1)\n",
    "        self.layer1 = self._make_layer(block, block_inplanes[0], layers[0],\n",
    "                                       shortcut_type)\n",
    "        self.layer2 = self._make_layer(block,\n",
    "                                       block_inplanes[1],\n",
    "                                       layers[1],\n",
    "                                       shortcut_type,\n",
    "                                       stride=2)\n",
    "        self.layer3 = self._make_layer(block,\n",
    "                                       block_inplanes[2],\n",
    "                                       layers[2],\n",
    "                                       shortcut_type,\n",
    "                                       stride=2)\n",
    "        self.layer4 = self._make_layer(block,\n",
    "                                       block_inplanes[3],\n",
    "                                       layers[3],\n",
    "                                       shortcut_type,\n",
    "                                       stride=2)\n",
    "\n",
    "        self.avgpool = nn.AdaptiveAvgPool3d((1, 1, 1))\n",
    "        self.fc_features = nn.Linear(512 + 851, 512)\n",
    "        self.fc = nn.Linear(block_inplanes[3] * block.expansion, n_classes)\n",
    "        self.dropout = nn.Dropout(dropout_rate)\n",
    "\n",
    "        for m in self.modules():\n",
    "            if isinstance(m, nn.Conv3d):\n",
    "                nn.init.kaiming_normal_(m.weight,\n",
    "                                        mode='fan_out',\n",
    "                                        nonlinearity='relu')\n",
    "            elif isinstance(m, nn.BatchNorm3d):\n",
    "                nn.init.constant_(m.weight, 1)\n",
    "                nn.init.constant_(m.bias, 0)\n",
    "\n",
    "    def _downsample_basic_block(self, x, planes, stride):\n",
    "        out = F.avg_pool3d(x, kernel_size=1, stride=stride)\n",
    "        zero_pads = torch.zeros(out.size(0), planes - out.size(1), out.size(2),\n",
    "                                out.size(3), out.size(4))\n",
    "        if isinstance(out.data, torch.cuda.FloatTensor):\n",
    "            zero_pads = zero_pads.cuda()\n",
    "\n",
    "        out = torch.cat([out.data, zero_pads], dim=1)\n",
    "\n",
    "        return out\n",
    "\n",
    "    def _make_layer(self, block, planes, blocks, shortcut_type, stride=1):\n",
    "        downsample = None\n",
    "        if stride != 1 or self.in_planes != planes * block.expansion:\n",
    "            if shortcut_type == 'A':\n",
    "                downsample = partial(self._downsample_basic_block,\n",
    "                                     planes=planes * block.expansion,\n",
    "                                     stride=stride)\n",
    "            else:\n",
    "                downsample = nn.Sequential(\n",
    "                    conv1x1x1(self.in_planes, planes * block.expansion, stride),\n",
    "                    nn.BatchNorm3d(planes * block.expansion))\n",
    "\n",
    "        layers = []\n",
    "        layers.append(\n",
    "            block(in_planes=self.in_planes,\n",
    "                  planes=planes,\n",
    "                  stride=stride,\n",
    "                  downsample=downsample))\n",
    "        self.in_planes = planes * block.expansion\n",
    "        for i in range(1, blocks):\n",
    "            layers.append(block(self.in_planes, planes))\n",
    "\n",
    "        return nn.Sequential(*layers)\n",
    "\n",
    "    def forward(self, x, additional_features):\n",
    "        x = self.dropout(x)\n",
    "        x = self.conv1(x)\n",
    "        x = self.bn1(x)\n",
    "        x = self.relu(x)\n",
    "        if not self.no_max_pool:\n",
    "            x = self.maxpool(x)\n",
    "\n",
    "        x = self.layer1(x)\n",
    "        x = self.dropout(x)\n",
    "        x = self.layer2(x)\n",
    "        x = self.dropout(x)\n",
    "        x = self.layer3(x)\n",
    "        x = self.dropout(x)\n",
    "        x = self.layer4(x)\n",
    "        x = self.dropout(x)\n",
    "\n",
    "        x = self.avgpool(x)\n",
    "\n",
    "        x = x.view(x.size(0), -1)\n",
    "        additional_features = additional_features.view(x.size(0), -1)\n",
    "        x = torch.cat((x, additional_features), dim=1)\n",
    "        x = F.relu(self.fc_features(x))\n",
    "        x = self.dropout(x)\n",
    "        x = self.fc(x)\n",
    "        x = torch.sigmoid(x)\n",
    "\n",
    "        return x\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "outputs": [],
   "source": [
    "\n",
    "def generate_model(model_depth, inplanes, **kwargs):\n",
    "    assert model_depth in [10, 18, 34, 50, 101, 152, 200]\n",
    "    model = None\n",
    "    if model_depth == 10:\n",
    "        model = ResNet(BasicBlock, [1, 1, 1, 1], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "    elif model_depth == 18:\n",
    "        model = ResNet(BasicBlock, [2, 2, 2, 2], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "    elif model_depth == 34:\n",
    "        model = ResNet(BasicBlock, [3, 4, 6, 3], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "    elif model_depth == 50:\n",
    "        model = ResNet(Bottleneck, [3, 4, 6, 3], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "    elif model_depth == 101:\n",
    "        model = ResNet(Bottleneck, [3, 4, 23, 3], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "    elif model_depth == 152:\n",
    "        model = ResNet(Bottleneck, [3, 8, 36, 3], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "    elif model_depth == 200:\n",
    "        model = ResNet(Bottleneck, [3, 24, 36, 3], inplanes, model_depth, widen_factor=widen_factor, **kwargs)\n",
    "\n",
    "    return model\n",
    "\n",
    "\n",
    "def get_model_name(trial, name, batch_size, learning_rate, dropout_rate, epoch):\n",
    "    \"\"\" Generate a name for the model consisting of all the hyperparameter values\n",
    "\n",
    "    Args:\n",
    "        config: Configuration object containing the hyperparameters\n",
    "    Returns:\n",
    "        path: A string with the hyperparameter name and value concatenated\n",
    "    \"\"\"\n",
    "    path = \"trial_{0}_model_{1}_bs{2}_lr{3}_dr{4}_epoch{5}\".format(trial,\n",
    "                                                                   name,\n",
    "                                                                   batch_size,\n",
    "                                                                   learning_rate,\n",
    "                                                                   dropout_rate,\n",
    "                                                                   epoch)\n",
    "    return path\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "outputs": [],
   "source": [
    "\n",
    "########################################################################\n",
    "# Other functions\n",
    "def create_label(mutation, fusion):\n",
    "    if mutation == 1:\n",
    "        return 1\n",
    "    elif fusion == 1:\n",
    "        return 0\n",
    "    else:\n",
    "        return None\n",
    "\n",
    "\n",
    "def process_excel(df_data, exclusions, custom_features=None):\n",
    "    nanmask = np.isnan(df_data[\"code\"])\n",
    "    data_data_new = df_data[~nanmask]\n",
    "    data_data_new = data_data_new.reindex()\n",
    "\n",
    "    # Remove exluded patients\n",
    "    data_data_new = data_data_new[~data_data_new[\"code\"].isin(exclusions)]\n",
    "    data_data_new = data_data_new.reindex()\n",
    "\n",
    "    # Remove data that we don't need for this analysis\n",
    "    data_data_new = data_data_new.drop(columns=['WT', 'NF1',\n",
    "                                                'CDKN2A (0=balanced, 1=Del, 2=Undetermined)', 'FGFR 1', 'FGFR 2',\n",
    "                                                'FGFR 4',\n",
    "                                                'Further gen info', 'Notes', 'Pathology Dx_Original', 'Pathology Coded',\n",
    "                                                'Location_1', 'Location_2', 'Location_Original', 'Gender', 'Age Dx'])\n",
    "\n",
    "    data_data_new['label'] = data_data_new.apply(lambda x: create_label(x['BRAF V600E final'], x['BRAF fusion final']),\n",
    "                                                 axis=1)\n",
    "    data_data_new = data_data_new.drop(columns=[\"BRAF V600E final\", \"BRAF fusion final\"])\n",
    "\n",
    "    # Drop rows where the outcome is not mutation or fusion\n",
    "    nanmask = np.isnan(data_data_new[\"label\"])\n",
    "    data_data_new = data_data_new[~nanmask]\n",
    "    data_data_new = data_data_new.reindex()\n",
    "    patient_codes = [int(x) for x in list(data_data_new[\"code\"].values)]\n",
    "\n",
    "    training_labels = dict(zip(patient_codes, list(data_data_new[\"label\"].values)))\n",
    "    data_data_new = data_data_new.drop(columns=[\"label\"])\n",
    "\n",
    "    # Organize the radiomic features into a dictionary with patient codes and corresponding patient features\n",
    "    if custom_features is not None:\n",
    "        radiomic_features = {}\n",
    "        for index, row in custom_features.iterrows():\n",
    "            radiomic_features[custom_features['id'][index]] = row.values[1:]\n",
    "    else:\n",
    "        data_data_new.set_index(\"code\", inplace=True)\n",
    "        radiomic_features = {}\n",
    "        for index, row in data_data_new.iterrows():\n",
    "            radiomic_features[index] = row.values\n",
    "    return radiomic_features, training_labels"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "outputs": [],
   "source": [
    "\n",
    "###############################################################################\n",
    "# Model Training\n",
    "def evaluate(net, loader, criterion=nn.BCEWithLogitsLoss()):\n",
    "    \"\"\" Evaluate the network on the validation set.\n",
    "\n",
    "     Args:\n",
    "         net: PyTorch neural network object\n",
    "         loader: PyTorch data loader for the validation set\n",
    "         criterion: The loss function\n",
    "     Returns:\n",
    "         err: A scalar for the avg classification error over the validation set\n",
    "         loss: A scalar for the average loss function over the validation set\n",
    "     \"\"\"\n",
    "    net.eval()\n",
    "    total_loss = 0.0\n",
    "    total_err = 0.0\n",
    "    total_epoch = 0\n",
    "    true = []\n",
    "    estimated = []\n",
    "    for i, data in enumerate(loader, 0):\n",
    "        inputs, labels, features = data\n",
    "        # print(labels.float())\n",
    "        # print(data)\n",
    "        # print(i)\n",
    "        use_cuda = True\n",
    "        if use_cuda and torch.cuda.is_available():\n",
    "            inputs = inputs.cuda()\n",
    "            labels = labels.cuda()\n",
    "            features = features.cuda()\n",
    "            net = net.cuda()\n",
    "        outputs = net(inputs, features)\n",
    "        # print(outputs.float())\n",
    "        loss_func = nn.BCEWithLogitsLoss()\n",
    "        loss = loss_func(outputs, labels.float())\n",
    "        corr = (outputs > 0.0).squeeze().long() != labels\n",
    "        total_err += int(corr.sum())\n",
    "        total_loss += loss.item()\n",
    "        total_epoch += len(labels)\n",
    "\n",
    "        for i in range(len(labels.tolist())):\n",
    "            true.append(labels.tolist()[i][0])\n",
    "            estimated.append(outputs.tolist()[i][0])\n",
    "    err = float(total_err) / total_epoch\n",
    "    loss = float(total_loss) / (i + 1)\n",
    "    return err, loss, true, estimated\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "outputs": [],
   "source": [
    "def train_net(train_dataloader, val_dataloader, test_dataloader, trial, net, optimizer, criterion, batch_size=64,\n",
    "              learning_rate=0.01, num_epochs=30, checkpoint=False,\n",
    "              save_folder=os.getcwd()):\n",
    "    # total_train_err = np.zeros(num_epochs)\n",
    "    total_train_loss = np.zeros(num_epochs)\n",
    "    total_train_auc = np.zeros(num_epochs)\n",
    "    # total_val_err = np.zeros(num_epochs)\n",
    "    total_val_loss = np.zeros(num_epochs)\n",
    "    total_val_auc = np.zeros(num_epochs)\n",
    "\n",
    "    total_test_loss = np.zeros(num_epochs)\n",
    "    total_test_auc = np.zeros(num_epochs)\n",
    "\n",
    "    total_train_roc = []\n",
    "    total_val_roc = []\n",
    "\n",
    "    training_start_time = time.time()\n",
    "\n",
    "    for epoch in range(num_epochs):\n",
    "        if len(total_val_loss) > 3:\n",
    "            if total_val_loss[-1] > total_val_loss[-2] > total_val_loss[-3]:\n",
    "                break\n",
    "\n",
    "        epoch_start_time = time.time()\n",
    "\n",
    "        # Training\n",
    "        # train_err = 0.0\n",
    "        net.train()\n",
    "        train_loss = 0.0\n",
    "        total_epoch = 0\n",
    "        training_true = []\n",
    "        training_estimated = []\n",
    "        n = 0\n",
    "        for inputs, labels, features in train_dataloader:\n",
    "            inputs, labels, features = inputs.to(device), labels.to(device), features.to(device)\n",
    "\n",
    "            # Add noise to images\n",
    "            noise = torch.randn_like(inputs, device=device) * 0.1\n",
    "            inputs = inputs + noise\n",
    "\n",
    "            # Forward + Backward + Optimize\n",
    "            optimizer.zero_grad()\n",
    "            outputs = net(inputs, features)\n",
    "            loss = criterion(outputs, labels.float())\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "\n",
    "            if use_scheduler:\n",
    "                scheduler.step()\n",
    "\n",
    "            # corr = (outputs > 0.0).squeeze().long() != labels\n",
    "            # train_err += int(corr.sum())\n",
    "            # Keep track of loss through the entire epoch\n",
    "            train_loss += loss.item()\n",
    "            total_epoch += len(labels)\n",
    "            n = n + 1\n",
    "\n",
    "            for i in range(len(labels.tolist())):\n",
    "                training_true.append(labels.tolist()[i][0])\n",
    "                training_estimated.append(outputs.tolist()[i][0])\n",
    "\n",
    "        # Calculate average over epoch\n",
    "        # total_train_err[epoch] = float(train_err) / total_epoch\n",
    "        total_train_loss[epoch] = float(train_loss) / n\n",
    "\n",
    "        # Validation\n",
    "        eval_err_val, eval_loss_val, validation_true, validation_estimated = evaluate(net, val_dataloader)\n",
    "        total_val_loss[epoch] = eval_loss_val\n",
    "\n",
    "        eval_err_test, eval_loss_test, test_true, test_estimated = evaluate(net, test_dataloader)\n",
    "        total_test_loss[epoch] = eval_loss_test\n",
    "\n",
    "        # Calculate the AUC for the different models\n",
    "        train_auc = roc_auc_score(training_true, training_estimated)\n",
    "        val_auc = roc_auc_score(validation_true, validation_estimated)\n",
    "        test_auc = roc_auc_score(test_true, test_estimated)\n",
    "\n",
    "        total_train_auc[epoch] = train_auc\n",
    "        total_val_auc[epoch] = val_auc\n",
    "        total_test_auc[epoch] = test_auc\n",
    "\n",
    "        train_fpr, train_tpr, _ = roc_curve(training_true, training_estimated)\n",
    "        val_fpr, val_tpr, _ = roc_curve(validation_true, validation_estimated)\n",
    "        test_fpr, test_tpr, _ = roc_curve(test_true, test_estimated)\n",
    "\n",
    "        # total_train_roc.append((train_fpr, train_tpr))\n",
    "        # total_val_roc.append((val_fpr, val_tpr))\n",
    "\n",
    "        # logging.info(\n",
    "        #     \"Epoch {}: Train err: {}, Train loss: {}, Train AUC: {} | Val err: {}, Val loss: {}, Val AUC: {} \".format(\n",
    "        #         epoch + 1,\n",
    "        #         total_train_err[epoch],\n",
    "        #         total_train_loss[epoch],\n",
    "        #         total_train_auc[epoch],\n",
    "        #         total_val_err[epoch],\n",
    "        #         total_val_loss[epoch],\n",
    "        #         total_val_auc[epoch]))\n",
    "\n",
    "        logging.info(\n",
    "            \"Epoch {}: Train loss: {}, Train AUC: {} | Val loss: {}, Val AUC: {} | Test loss: {}, Test AUC: {}\".format(\n",
    "                epoch + 1,\n",
    "                total_train_loss[epoch],\n",
    "                total_train_auc[epoch],\n",
    "                total_val_loss[epoch],\n",
    "                total_val_auc[epoch],\n",
    "                total_test_loss[epoch],\n",
    "                total_test_auc[epoch]))\n",
    "\n",
    "        model_path = get_model_name(trial=trial, name=net.name, batch_size=batch_size, learning_rate=learning_rate,\n",
    "                                    dropout_rate=dropout_rate, epoch=epoch + 1)\n",
    "\n",
    "        if checkpoint:\n",
    "            torch.save(net.state_dict(), model_path)\n",
    "\n",
    "    # np.savetxt(os.path.join(save_folder, \"{}_train_err.csv\".format(model_path)), total_train_err)\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_train_loss.csv\".format(model_path)), total_train_loss)\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_train_auc.csv\".format(model_path)), total_train_auc)\n",
    "    # np.savetxt(os.path.join(save_folder, \"{}_val_err.csv\".format(model_path)), total_val_err)\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_val_loss.csv\".format(model_path)), total_val_loss)\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_val_auc.csv\".format(model_path)), total_val_auc)\n",
    "\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_test_loss.csv\".format(model_path)), total_test_loss)\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_test_auc.csv\".format(model_path)), total_test_auc)\n",
    "\n",
    "    # with open(os.path.join(save_folder, \"{}_train_roc.csv\".format(model_path)), 'wb') as csvfile:\n",
    "    #     fwriter = csv.writer(csvfile)\n",
    "    #     for x in total_train_roc:\n",
    "    #         fwriter.writerow(x)\n",
    "    #\n",
    "    # with open(os.path.join(save_folder, \"{}_val_roc.csv\".format(model_path)), 'wb') as csvfile:\n",
    "    #     fwriter = csv.writer(csvfile)\n",
    "    #     for x in total_val_roc:\n",
    "    #         fwriter.writerow(x)\n",
    "\n",
    "    logging.info('Finished training.')\n",
    "    logging.info(f'Time elapsed: {round(time.time() - training_start_time, 3)} seconds.')\n",
    "\n",
    "    # return total_train_err, total_train_loss, total_train_auc, total_val_err, total_val_loss, total_val_auc\n",
    "    return total_train_loss, total_train_auc, total_val_loss, total_val_auc, total_test_loss, total_test_auc\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "outputs": [],
   "source": [
    "def remove_correlated_features(X=None, threshold=0.95):\n",
    "    if X is not None:\n",
    "        correlation_matrix = X.corr().abs()\n",
    "        upper_triangle = correlation_matrix.where(np.triu(np.ones(correlation_matrix.shape), k=1).astype(bool))\n",
    "        to_drop = [column for column in upper_triangle.columns if any(upper_triangle[column] > threshold)]\n",
    "        X = X.drop(to_drop, axis='columns')\n",
    "    return X\n",
    "\n",
    "\n",
    "def variance_threshold(X):\n",
    "    selector = VarianceThreshold(threshold=0.05)\n",
    "    ids = X['id']\n",
    "    X = selector.fit_transform(X.loc[:, X.columns != 'id'])\n",
    "    X = pd.DataFrame(X)\n",
    "    X['id'] = ids\n",
    "    return X"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "########################################################################\n",
    "# Currently only SK images are used, so \"test_dataloader\" represents\n",
    "# internal testing data. Once Stanford images are processed, they will\n",
    "# also be a \"testing set\" but rather \"external testing data\".\n",
    "########################################################################\n",
    "# run\n",
    "if __name__ == '__main__':\n",
    "    start_up_time = time.time()\n",
    "\n",
    "    save_folder = os.path.join(os.getcwd(), f\"CNN_results_{time.strftime('%Y_%m_%d-%H_%M_%S')}\")\n",
    "    os.mkdir(save_folder)\n",
    "    targets = logging.StreamHandler(sys.stdout), logging.FileHandler(os.path.join(save_folder, 'output_log.log'))\n",
    "    logging.basicConfig(format='%(message)s', level=logging.INFO, handlers=targets)\n",
    "\n",
    "    random_seed(random.randint(1, 100), True)\n",
    "    pd.set_option('display.max_rows', None)\n",
    "\n",
    "    # use numpy files instead of .nii\n",
    "    # no need to normalize images between [0,1] as input images are already preprocessed\n",
    "    # https://github.com/kenshohara/3D-ResNets-PyTorch\n",
    "\n",
    "    radiomics_directory = r'C:\\Users\\Justin\\Documents\\Data'\n",
    "    image_directory = r'K:\\Projects\\SickKids_Brain_Preprocessing\\preprocessed_all_seq_kk_july_2022'\n",
    "\n",
    "    # Parameters\n",
    "    load_model = False\n",
    "    use_scheduler = False\n",
    "    limit = 5\n",
    "\n",
    "    num_trials = 2\n",
    "    num_epochs = 2\n",
    "    batch_size = 8\n",
    "    learning_rate = 0.01\n",
    "    dropout_rate = 0.5  # default\n",
    "    widen_factor = 0.5  # speeds up training without affect performance too much\n",
    "    inplanes = [64, 128, 256, 512]\n",
    "\n",
    "    excluded_patients = [2, 3, 4, 6, 7, 9, 11, 12, 13, 16, 21, 23, 25, 28, 29, 30, 33, 35, 36, 37, 38, 44, 45, 49, 50,\n",
    "                         52, 53, 54, 55, 58, 59, 61, 63, 66, 73, 74, 75, 77, 78, 80, 84, 85, 86, 92, 95, 96, 98, 100,\n",
    "                         102, 103, 105, 107, 108, 110, 113, 117, 121, 122, 123, 125, 128, 130, 131, 132, 136, 137, 138,\n",
    "                         139, 140, 141, 142, 143, 147, 148, 150, 152, 156, 158, 159, 165, 166, 171, 173, 174, 176, 182,\n",
    "                         183, 184, 187, 190, 191, 192, 194, 196, 199, 203, 204, 209, 210, 213, 221, 222, 224, 226, 227,\n",
    "                         228, 232, 233, 234, 235, 237, 240, 242, 243, 245, 246, 250, 254, 255, 256, 258, 259, 260, 261,\n",
    "                         263, 264, 266, 270, 272, 274, 277, 278, 283, 284, 285, 288, 293, 298, 299, 303, 306, 309, 311,\n",
    "                         312, 317, 318, 321, 322, 324, 325, 327, 328, 330, 332, 333, 334, 336, 337, 341, 343, 347,\n",
    "                         349, 350, 351, 352, 354, 356, 359, 364, 367, 370, 371, 374, 376, 377, 378, 380, 383, 386, 387,\n",
    "                         388, 392, 396, 243, 255, 261, 264, 288, 299, 309, 327, 351, 387]\n",
    "\n",
    "    if torch.cuda.is_available():\n",
    "        device = \"cuda\"\n",
    "    else:\n",
    "        device = \"cpu\"\n",
    "\n",
    "    # Load data\n",
    "    df_sickkids = load_excel_data(os.path.join(radiomics_directory, 'Nomogram_study_LGG_data_Nov.27.xlsx'), sheet='SK')\n",
    "    df_features = pd.read_csv(r'C:\\Users\\Justin\\Documents\\Data\\radiomics_features_normalized_08-15-22_filtered_851.csv')\n",
    "\n",
    "    df_features = remove_correlated_features(df_features)\n",
    "    df_features = variance_threshold(df_features)\n",
    "\n",
    "    sickkids_radiomics_features, sickkids_labels = process_excel(df_data=df_sickkids, exclusions=excluded_patients,\n",
    "                                                                 custom_features=df_features)\n",
    "\n",
    "    # Prepare CNN data\n",
    "    radiomics_patients_list = set(sickkids_labels.keys())\n",
    "    patients_with_FLAIR = []\n",
    "    for each_patient in os.listdir(image_directory):\n",
    "        try:\n",
    "            patients_with_FLAIR.append(int(each_patient))\n",
    "        except:\n",
    "            logging.info(f'Patient {each_patient} FLAIR not found.')\n",
    "    patients_with_FLAIR.sort(key=int)\n",
    "    patients_list = list(radiomics_patients_list.intersection(patients_with_FLAIR))\n",
    "    logging.info(f\"Total number of patients: {len(patients_list)}.\")\n",
    "    logging.info(f\"Start-up time: {round(time.time() - start_up_time, 3)} seconds.\\n\")\n",
    "\n",
    "    load_image_time = time.time()\n",
    "    images, patients_used = load_image_data(image_directory, patients=patients_list, limit=limit)\n",
    "    data_pre_norm = {}\n",
    "    for each_patient in patients_used:\n",
    "        image = images[each_patient][0]\n",
    "        seg = images[each_patient][1]\n",
    "        data_pre_norm[each_patient] = (image, seg)\n",
    "    data_post_norm = {}\n",
    "    for each_patient in data_pre_norm.keys():\n",
    "        image_pre_norm = data_pre_norm[each_patient][0]\n",
    "        seg = data_pre_norm[each_patient][1]\n",
    "        image_post_norm = (image_pre_norm - np.min(image_pre_norm)) / (np.max(image_pre_norm) - np.min(image_pre_norm))\n",
    "        data_post_norm[each_patient] = (image_post_norm, seg)\n",
    "    data = {}\n",
    "    for each_patient in patients_used:\n",
    "        input = torch.tensor(\n",
    "            np.multiply(data_post_norm[each_patient][0], data_post_norm[each_patient][1])).float().unsqueeze(0)\n",
    "        label = sickkids_labels[each_patient]\n",
    "        label = torch.tensor(label).float().unsqueeze(0)\n",
    "        features = torch.tensor(sickkids_radiomics_features[each_patient].tolist()).float().unsqueeze(0)\n",
    "        patient = {\n",
    "            \"input\": input,\n",
    "            \"label\": label,\n",
    "            \"features\": features\n",
    "        }\n",
    "        data[each_patient] = patient\n",
    "\n",
    "    logging.info(\"Done loading images.\")\n",
    "    logging.info(f\"Number of patients included: {len(patients_used)}.\")\n",
    "    logging.info(f\"Image loading time: {round(time.time() - load_image_time, 3)} seconds.\\n\")\n",
    "\n",
    "    # if load_model:\n",
    "    #     try:\n",
    "    #         net = generate_model(model_depth=18, inplanes=inplanes, n_classes=1039)\n",
    "    #         model_path = get_model_name(trial=1, name=net.name, batch_size=batch_size, learning_rate=learning_rate,\n",
    "    #                                     dropout_rate=dropout_rate, epoch=num_epochs)\n",
    "    #         state = torch.load(model_path)\n",
    "    #         net.load_state_dict(state)\n",
    "    #     except FileNotFoundError:\n",
    "    #         logging.info('Model not found.')\n",
    "    #     else:\n",
    "    #         logging.info(\"Insert code...\")\n",
    "    #     sys.exit()\n",
    "\n",
    "    training_aucs = []\n",
    "    validation_aucs = []\n",
    "    test_aucs = []\n",
    "    best_epochs = []\n",
    "    trial_times = []\n",
    "\n",
    "    for t in range(num_trials):\n",
    "        logging.info(f\"Beginning trial {t + 1} of {num_trials}...\")\n",
    "        begin_trial_time = time.time()\n",
    "\n",
    "        # Set the seed for this iteration\n",
    "        if t == 0:\n",
    "            random_seed(1, True)\n",
    "            next_seed = random.randint(0, 1000)\n",
    "        else:\n",
    "            random_seed(next_seed, True)\n",
    "            next_seed = random.randint(0, 1000)\n",
    "\n",
    "        dataset = CNNDataset(data, patients_used)\n",
    "        train_size = int(0.6 * len(dataset))\n",
    "        validation_size = int(0.2 * len(dataset))\n",
    "        test_size = len(dataset) - train_size - validation_size\n",
    "        train_dataset, validation_dataset, test_dataset = torch.utils.data.random_split(dataset, [train_size,\n",
    "                                                                                                  validation_size,\n",
    "                                                                                                  test_size])\n",
    "        train_dataloader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True, num_workers=0)\n",
    "        validation_dataloader = DataLoader(validation_dataset, batch_size=batch_size, shuffle=True, num_workers=0)\n",
    "        test_dataloader = DataLoader(test_dataset, batch_size=batch_size, shuffle=True, num_workers=0)\n",
    "        logging.info(f\"Datasplit -> Training: {train_size}, Validation: {validation_size}, Testing: {test_size}.\")\n",
    "\n",
    "        net = generate_model(model_depth=18, inplanes=inplanes, n_classes=1039)\n",
    "\n",
    "        net.conv1 = nn.Conv3d(1, int(64 * widen_factor), kernel_size=(7, 7, 7), stride=(1, 2, 2), padding=(3, 3, 3),\n",
    "                              bias=False)\n",
    "        net.fc = net.fc = nn.Linear(512 * widen_factor, 1)\n",
    "\n",
    "        net.to(device)\n",
    "\n",
    "        criterion = nn.BCEWithLogitsLoss()\n",
    "        optimizer = optim.Adam(net.parameters(), lr=learning_rate)\n",
    "\n",
    "        if use_scheduler:\n",
    "            scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=75, gamma=0.1)\n",
    "\n",
    "        results = train_net(train_dataloader, validation_dataloader, test_dataloader,\n",
    "                            trial=t + 1,\n",
    "                            net=net,\n",
    "                            optimizer=optimizer,\n",
    "                            criterion=criterion,\n",
    "                            batch_size=batch_size,\n",
    "                            learning_rate=learning_rate,\n",
    "                            num_epochs=num_epochs, save_folder=save_folder)\n",
    "\n",
    "        # train_err, train_loss, train_auc, val_err, val_loss, val_auc = results\n",
    "        train_loss, train_auc, val_loss, val_auc, test_loss, test_auc = results\n",
    "\n",
    "        epoch = np.where(val_loss == min(val_loss))[0]\n",
    "        best_epochs.append(epoch[0])\n",
    "\n",
    "        best_train_auc = train_auc[np.where(val_loss == min(val_loss))]\n",
    "        training_aucs.append(best_train_auc[0])\n",
    "\n",
    "        best_val_auc = val_auc[np.where(val_loss == min(val_loss))]\n",
    "        validation_aucs.append(best_val_auc[0])\n",
    "\n",
    "        best_test_auc = test_auc[np.where(val_loss == min(val_loss))]\n",
    "        test_aucs.append(best_test_auc[0])\n",
    "\n",
    "        # logging.info(f\"Best epoch (lowest validation loss): {epoch[0]}, \"\n",
    "        #              f\"Lowest training error {round(min(train_err), 3)}, \"\n",
    "        #              f\"Lowest training loss {round(min(train_loss), 3)}, \"\n",
    "        #              f\"Training AUC corresponding to loss: {round(best_train_auc[0], 3)}, \"\n",
    "        #              f\"Lowest validation error {round(min(train_err), 3)}, \"\n",
    "        #              f\"Lowest validation loss {round(min(val_loss), 3)}, \"\n",
    "        #              f\"Validation AUC corresponding to loss: {round(best_val_auc[0], 3)}\")\n",
    "\n",
    "        logging.info(f\"Best epoch (lowest validation loss): {epoch[0]}, \"\n",
    "                     f\"Lowest training loss {round(min(train_loss), 3)}, \"\n",
    "                     f\"Lowest validation loss {round(min(val_loss), 3)}, \"\n",
    "                     f\"Training AUC corresponding to lowest validation loss: {round(best_train_auc[0], 3)}, \"\n",
    "                     f\"Validation AUC corresponding to lowest validation loss: {round(best_val_auc[0], 3)}, \"\n",
    "                     f\"Testing AUC corresponding to lowest validation loss: {round(best_test_auc[0], 3)}\")\n",
    "\n",
    "        trial_duration = time.time() - begin_trial_time\n",
    "        trial_times.append(round(trial_duration, 3))\n",
    "        logging.info(f\"Trial {t + 1} ended. Duration: {round(trial_duration, 3)} seconds.\\n\")\n",
    "    logging.info(f'Experiment done. Time elapsed: {round(time.time() - start_up_time, 3)} seconds.')\n",
    "logging.info('---------------------')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}