main-fusion-risk_score.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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",
    "from torch.utils.data import SubsetRandomSampler\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": 2,
   "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": 3,
   "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": 4,
   "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\"],\n",
    "        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": 5,
   "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",
    "        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 = 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):\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",
    "        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": 6,
   "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": 7,
   "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": 8,
   "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",
    "    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",
    "            net = net.cuda()\n",
    "        outputs = net(inputs)\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": 9,
   "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": 12,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading Nomogram_study_LGG_data_Nov.27.xlsx, Sheet: SK...\n",
      "Done loading.\n",
      "Total number of patients: 214.\n",
      "Start-up time: 16.208 seconds.\n",
      "\n",
      "Loading Patient 1...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Justin\\anaconda3\\envs\\pLGG_Predict_Status\\lib\\site-packages\\sklearn\\feature_selection\\_variance_threshold.py:104: RuntimeWarning: Degrees of freedom <= 0 for slice.\n",
      "  self.variances_ = np.nanvar(X, axis=0)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading Patient 5...\n",
      "Loading Patient 8...\n",
      "Loading Patient 10...\n",
      "Loading Patient 14...\n",
      "Loading Patient 15...\n",
      "Loading Patient 17...\n",
      "Loading Patient 18...\n",
      "Loading Patient 19...\n",
      "Loading Patient 20...\n",
      "Loading Patient 22...\n",
      "Loading Patient 24...\n",
      "Loading Patient 26...\n",
      "Loading Patient 27...\n",
      "Loading Patient 31...\n",
      "Loading Patient 32...\n",
      "Loading Patient 34...\n",
      "Loading Patient 39...\n",
      "Loading Patient 40...\n",
      "Loading Patient 41...\n",
      "Done loading images.\n",
      "Number of patients included: 20.\n",
      "Image loading time: 36.242 seconds.\n",
      "\n",
      "Beginning trial 1 of 2...\n",
      "Datasplit -> Development: 15, Testing: 5.\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[1;31mKeyboardInterrupt\u001B[0m                         Traceback (most recent call last)",
      "Input \u001B[1;32mIn [14]\u001B[0m, in \u001B[0;36m<cell line: 210>\u001B[1;34m()\u001B[0m\n\u001B[0;32m    381\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m use_scheduler:\n\u001B[0;32m    382\u001B[0m     scheduler \u001B[38;5;241m=\u001B[39m optim\u001B[38;5;241m.\u001B[39mlr_scheduler\u001B[38;5;241m.\u001B[39mStepLR(optimizer, step_size\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m75\u001B[39m, gamma\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.1\u001B[39m)\n\u001B[1;32m--> 384\u001B[0m results \u001B[38;5;241m=\u001B[39m \u001B[43mtrain_net\u001B[49m\u001B[43m(\u001B[49m\u001B[43mdevelopment_dataset\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mtest_dataset\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    385\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mtrial\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mt\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;241;43m+\u001B[39;49m\u001B[43m \u001B[49m\u001B[38;5;241;43m1\u001B[39;49m\u001B[43m,\u001B[49m\n\u001B[0;32m    386\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mnet\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mnet\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    387\u001B[0m \u001B[43m                    \u001B[49m\u001B[43moptimizer\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43moptimizer\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    388\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mcriterion\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mcriterion\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    389\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mbatch_size\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mbatch_size\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    390\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mlearning_rate\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mlearning_rate\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    391\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mnum_epochs\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mnum_epochs\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43msave_folder\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43msave_folder\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    392\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mradiomics_grid_parameters\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mradiomics_grid_parameters\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m    393\u001B[0m \u001B[43m                    \u001B[49m\u001B[43mfusion_grid_parameters\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mfusion_grid_parameters\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m    395\u001B[0m dev_fusion_auc, test_fusion_auc \u001B[38;5;241m=\u001B[39m results\n\u001B[0;32m    397\u001B[0m df_trial \u001B[38;5;241m=\u001B[39m pd\u001B[38;5;241m.\u001B[39mDataFrame(\u001B[38;5;28mlist\u001B[39m(\u001B[38;5;28mzip\u001B[39m([dev_fusion_auc], [test_fusion_auc])),\n\u001B[0;32m    398\u001B[0m                         columns\u001B[38;5;241m=\u001B[39m[\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mDevelopment\u001B[39m\u001B[38;5;124m'\u001B[39m, \u001B[38;5;124m'\u001B[39m\u001B[38;5;124mTesting\u001B[39m\u001B[38;5;124m'\u001B[39m])\n",
      "Input \u001B[1;32mIn [14]\u001B[0m, in \u001B[0;36mtrain_net\u001B[1;34m(dev_dataset, test_dataset, trial, net, optimizer, criterion, batch_size, learning_rate, num_epochs, checkpoint, save_folder, radiomics_grid_parameters, fusion_grid_parameters)\u001B[0m\n\u001B[0;32m      8\u001B[0m \u001B[38;5;66;03m# fit radiomics model on dev data\u001B[39;00m\n\u001B[0;32m      9\u001B[0m rfc \u001B[38;5;241m=\u001B[39m RandomForestClassifier(n_jobs\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m-\u001B[39m\u001B[38;5;241m1\u001B[39m)  \u001B[38;5;66;03m# all CPUs used\u001B[39;00m\n\u001B[1;32m---> 10\u001B[0m rfc_cv \u001B[38;5;241m=\u001B[39m GridSearchCV(estimator\u001B[38;5;241m=\u001B[39mrfc, param_grid\u001B[38;5;241m=\u001B[39mradiomics_grid_parameters, verbose\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0\u001B[39m, scoring\u001B[38;5;241m=\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mroc_auc\u001B[39m\u001B[38;5;124m\"\u001B[39m,\n\u001B[0;32m     11\u001B[0m                       n_jobs\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m-\u001B[39m\u001B[38;5;241m1\u001B[39m, return_train_score\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m, refit\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m)\n\u001B[0;32m     13\u001B[0m dev_patient_data \u001B[38;5;241m=\u001B[39m \u001B[38;5;28mlist\u001B[39m(dev_dataset[\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mdata\u001B[39m\u001B[38;5;124m'\u001B[39m]\u001B[38;5;241m.\u001B[39mvalues())\n\u001B[0;32m     14\u001B[0m Y_dev \u001B[38;5;241m=\u001B[39m np\u001B[38;5;241m.\u001B[39mconcatenate([d[\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mlabel\u001B[39m\u001B[38;5;124m'\u001B[39m]\u001B[38;5;241m.\u001B[39mnumpy() \u001B[38;5;28;01mfor\u001B[39;00m d \u001B[38;5;129;01min\u001B[39;00m dev_patient_data])\n",
      "Input \u001B[1;32mIn [14]\u001B[0m, in \u001B[0;36mtrain_net\u001B[1;34m(dev_dataset, test_dataset, trial, net, optimizer, criterion, batch_size, learning_rate, num_epochs, checkpoint, save_folder, radiomics_grid_parameters, fusion_grid_parameters)\u001B[0m\n\u001B[0;32m      8\u001B[0m \u001B[38;5;66;03m# fit radiomics model on dev data\u001B[39;00m\n\u001B[0;32m      9\u001B[0m rfc \u001B[38;5;241m=\u001B[39m RandomForestClassifier(n_jobs\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m-\u001B[39m\u001B[38;5;241m1\u001B[39m)  \u001B[38;5;66;03m# all CPUs used\u001B[39;00m\n\u001B[1;32m---> 10\u001B[0m rfc_cv \u001B[38;5;241m=\u001B[39m GridSearchCV(estimator\u001B[38;5;241m=\u001B[39mrfc, param_grid\u001B[38;5;241m=\u001B[39mradiomics_grid_parameters, verbose\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0\u001B[39m, scoring\u001B[38;5;241m=\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mroc_auc\u001B[39m\u001B[38;5;124m\"\u001B[39m,\n\u001B[0;32m     11\u001B[0m                       n_jobs\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m-\u001B[39m\u001B[38;5;241m1\u001B[39m, return_train_score\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m, refit\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m)\n\u001B[0;32m     13\u001B[0m dev_patient_data \u001B[38;5;241m=\u001B[39m \u001B[38;5;28mlist\u001B[39m(dev_dataset[\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mdata\u001B[39m\u001B[38;5;124m'\u001B[39m]\u001B[38;5;241m.\u001B[39mvalues())\n\u001B[0;32m     14\u001B[0m Y_dev \u001B[38;5;241m=\u001B[39m np\u001B[38;5;241m.\u001B[39mconcatenate([d[\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mlabel\u001B[39m\u001B[38;5;124m'\u001B[39m]\u001B[38;5;241m.\u001B[39mnumpy() \u001B[38;5;28;01mfor\u001B[39;00m d \u001B[38;5;129;01min\u001B[39;00m dev_patient_data])\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_39_64.pyx:1179\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_39_64.SafeCallWrapper.__call__\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_39_64.pyx:620\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_39_64.PyDBFrame.trace_dispatch\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_39_64.pyx:1095\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_39_64.PyDBFrame.trace_dispatch\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_39_64.pyx:1057\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_39_64.PyDBFrame.trace_dispatch\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_39_64.pyx:317\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_39_64.PyDBFrame.do_wait_suspend\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32mC:\\Program Files\\JetBrains\\PyCharm 2022.2\\plugins\\python\\helpers\\pydev\\pydevd.py:1160\u001B[0m, in \u001B[0;36mPyDB.do_wait_suspend\u001B[1;34m(self, thread, frame, event, arg, send_suspend_message, is_unhandled_exception)\u001B[0m\n\u001B[0;32m   1157\u001B[0m         from_this_thread\u001B[38;5;241m.\u001B[39mappend(frame_id)\n\u001B[0;32m   1159\u001B[0m \u001B[38;5;28;01mwith\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_threads_suspended_single_notification\u001B[38;5;241m.\u001B[39mnotify_thread_suspended(thread_id, stop_reason):\n\u001B[1;32m-> 1160\u001B[0m     \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_do_wait_suspend\u001B[49m\u001B[43m(\u001B[49m\u001B[43mthread\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mframe\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mevent\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43marg\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43msuspend_type\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mfrom_this_thread\u001B[49m\u001B[43m)\u001B[49m\n",
      "File \u001B[1;32mC:\\Program Files\\JetBrains\\PyCharm 2022.2\\plugins\\python\\helpers\\pydev\\pydevd.py:1175\u001B[0m, in \u001B[0;36mPyDB._do_wait_suspend\u001B[1;34m(self, thread, frame, event, arg, suspend_type, from_this_thread)\u001B[0m\n\u001B[0;32m   1172\u001B[0m             \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_call_mpl_hook()\n\u001B[0;32m   1174\u001B[0m         \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mprocess_internal_commands()\n\u001B[1;32m-> 1175\u001B[0m         \u001B[43mtime\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43msleep\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;241;43m0.01\u001B[39;49m\u001B[43m)\u001B[49m\n\u001B[0;32m   1177\u001B[0m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mcancel_async_evaluation(get_current_thread_id(thread), \u001B[38;5;28mstr\u001B[39m(\u001B[38;5;28mid\u001B[39m(frame)))\n\u001B[0;32m   1179\u001B[0m \u001B[38;5;66;03m# process any stepping instructions\u001B[39;00m\n",
      "\u001B[1;31mKeyboardInterrupt\u001B[0m: "
     ]
    }
   ],
   "source": [
    "def train_net(dev_dataset, test_dataset, trial, net, optimizer, criterion, batch_size=64,\n",
    "              learning_rate=0.01, num_epochs=30, checkpoint=False,\n",
    "              save_folder=os.getcwd(), radiomics_grid_parameters=dict, fusion_grid_parameters=dict):\n",
    "    training_start_time = time.time()\n",
    "\n",
    "    kfold = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)\n",
    "\n",
    "    # fit radiomics model on dev data\n",
    "    rfc = RandomForestClassifier(n_jobs=-1)  # all CPUs used\n",
    "    rfc_cv = GridSearchCV(estimator=rfc, param_grid=radiomics_grid_parameters, verbose=0, scoring=\"roc_auc\",\n",
    "                          n_jobs=-1, return_train_score=True, refit=True)\n",
    "\n",
    "    dev_patient_data = list(dev_dataset['data'].values())\n",
    "    Y_dev = np.concatenate([d['label'].numpy() for d in dev_patient_data])\n",
    "\n",
    "    test_patient_data = list(test_dataset['data'].values())\n",
    "    Y_test = np.concatenate([d['label'].numpy() for d in test_patient_data])\n",
    "\n",
    "    test_loader = DataLoader(test_dataset, batch_size=batch_size, shuffle=True, num_workers=0)\n",
    "\n",
    "    all_dev_predictions_from_radiomics = {}\n",
    "    all_test_predictions_from_radiomics = {}\n",
    "    all_dev_predictions_from_resnet = {}\n",
    "    all_test_predictions_from_resnet = {}\n",
    "\n",
    "    all_dev_true_split_by_fold = {}\n",
    "    for fold, (train_idx, val_idx) in enumerate(kfold.split(np.arange(len(dev_dataset)))):\n",
    "        print('Fold {}'.format(fold + 1))\n",
    "\n",
    "        train_sampler = SubsetRandomSampler(train_idx)\n",
    "        val_sampler = SubsetRandomSampler(val_idx)\n",
    "\n",
    "        #RADIOMICS\n",
    "        train_patient_data = list(train_idx['data'].values())\n",
    "        X_train = pd.DataFrame([d['features'].numpy() for d in train_patient_data])\n",
    "        Y_train = np.concatenate([d['label'].numpy() for d in train_patient_data])\n",
    "        fit_model = rfc_cv.fit(X_train, Y_train)\n",
    "        best_model_radiomics = fit_model.best_estimator_\n",
    "\n",
    "        # get val predictions for best radiomics model\n",
    "        val_patient_data = list(val_idx['data'].values())\n",
    "        X_val = pd.DataFrame([d['features'].numpy() for d in val_patient_data])\n",
    "        Y_val = np.concatenate([d['label'].numpy() for d in val_patient_data])\n",
    "        all_dev_predictions_from_radiomics[fold] = best_model_radiomics.predict_proba(X_val)[:, 1]\n",
    "\n",
    "        # get test predictions for radiomics model\n",
    "        test_patient_data = list(test_dataset['data'].values())\n",
    "        X_test = pd.DataFrame([d['features'].numpy() for d in test_patient_data])\n",
    "        Y_test = np.concatenate([d['label'].numpy() for d in test_patient_data])\n",
    "        all_test_predictions_from_radiomics[fold] = best_model_radiomics.predict_proba(X_test)[:, 1]\n",
    "\n",
    "        #CNN\n",
    "        train_loader = DataLoader(dev_dataset, batch_size=batch_size, sampler=train_sampler)\n",
    "        val_loader = DataLoader(dev_dataset, batch_size=batch_size, sampler=val_sampler)\n",
    "\n",
    "        train_history = {'train_auc': np.zeros(num_epochs), 'train_loss': np.zeros(num_epochs),\n",
    "                         'val_auc': np.zeros(num_epochs), 'val_loss': np.zeros(num_epochs),\n",
    "                         'test_auc': np.zeros(num_epochs), 'test_loss': np.zeros(num_epochs)}\n",
    "\n",
    "        dev_predictions_from_resnet_epochs = {}\n",
    "        test_predictions_from_resnet_epochs = {}\n",
    "\n",
    "        model_results_name = get_model_name(trial=trial, name=net.name, batch_size=batch_size,\n",
    "                                            learning_rate=learning_rate, dropout_rate=dropout_rate, epoch=num_epochs)\n",
    "        for epoch in range(num_epochs):\n",
    "            if len(train_history['val_loss']) > 3:\n",
    "                if train_history['val_loss'][-1] > train_history['val_loss'][-2] > train_history['val_loss'][-3]:\n",
    "                    break\n",
    "\n",
    "            epoch_start_time = time.time()\n",
    "\n",
    "            # Training\n",
    "            # train_err = 0.0\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_loader:\n",
    "                inputs, labels = inputs.to(device), labels.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)\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",
    "            train_history['train_loss'][epoch] = float(train_loss) / n\n",
    "\n",
    "            # Validation\n",
    "            eval_err_val, eval_loss_val, validation_true, validation_estimated = evaluate(net, val_loader)\n",
    "            train_history['val_loss'][epoch] = eval_loss_val\n",
    "\n",
    "            eval_err_test, eval_loss_test, test_true, test_estimated = evaluate(net, test_loader)\n",
    "            train_history['test_loss'][epoch] = eval_loss_test\n",
    "\n",
    "            dev_predictions_from_resnet_epochs[epoch] = validation_estimated\n",
    "            test_predictions_from_resnet_epochs[epoch] = test_estimated\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",
    "            train_history['train_auc'][epoch] = train_auc\n",
    "            train_history['val_auc'][epoch] = val_auc\n",
    "            train_history['test_auc'][epoch] = test_auc\n",
    "\n",
    "            logging.info(\n",
    "                \"Epoch {}: Train loss: {}, Train AUC: {} | Val loss: {}, Val AUC: {} | Test loss: {}, Test AUC: {}\".format(\n",
    "                    epoch + 1,\n",
    "                    train_history['train_loss'][epoch],\n",
    "                    train_history['train_auc'][epoch],\n",
    "                    train_history['val_loss'][epoch],\n",
    "                    train_history['val_auc'][epoch],\n",
    "                    train_history['test_loss'][epoch],\n",
    "                    train_history['test_auc'][epoch]))\n",
    "\n",
    "            if checkpoint:\n",
    "                torch.save(net.state_dict(), model_results_name)\n",
    "\n",
    "        all_epochs_val_loss = train_history['val_loss']\n",
    "        best_epoch = np.where(all_epochs_val_loss == min(all_epochs_val_loss))[0]\n",
    "\n",
    "        all_dev_predictions_from_resnet[fold] = dev_predictions_from_resnet_epochs[best_epoch]\n",
    "        all_test_predictions_from_resnet[fold] = test_predictions_from_resnet_epochs[best_epoch]\n",
    "\n",
    "        all_dev_true_split_by_fold[fold] = Y_val\n",
    "\n",
    "    rfc_fusion = RandomForestClassifier(n_jobs=-1)  # all CPUs used\n",
    "    rfc_cv_fusion = GridSearchCV(estimator=rfc_fusion, param_grid=fusion_grid_parameters, cv=kfold, verbose=0,\n",
    "                                 scoring=\"roc_auc\", n_jobs=-1, return_train_score=True, refit=True)\n",
    "\n",
    "    val_aucs_by_fold_radiomics = []\n",
    "    for fold in all_dev_predictions_from_radiomics.keys():\n",
    "        val_aucs_by_fold_radiomics.append((fold, roc_auc_score(all_dev_true_split_by_fold[fold],\n",
    "                                                               all_dev_predictions_from_radiomics[fold])))\n",
    "    best_fold_radiomics = min(val_aucs_by_fold_radiomics, key=lambda x: x[1])[0]\n",
    "    test_estimated_radiomics_best_fold = all_test_predictions_from_radiomics[best_fold_radiomics]\n",
    "\n",
    "    val_aucs_by_fold_resnet = []\n",
    "    for fold in all_dev_predictions_from_resnet.keys():\n",
    "        val_aucs_by_fold_resnet.append((fold, roc_auc_score(all_dev_true_split_by_fold[fold],\n",
    "                                                            all_dev_predictions_from_resnet[fold])))\n",
    "    best_fold_resnet = min(val_aucs_by_fold_resnet, key=lambda x: x[1])[0]\n",
    "    test_estimated_resnet_best_fold = all_test_predictions_from_resnet[best_fold_resnet]\n",
    "\n",
    "    dev_predictions_radiomics = []\n",
    "    for fold in all_dev_predictions_from_radiomics.keys():\n",
    "        dev_predictions_radiomics = dev_predictions_radiomics + all_dev_predictions_from_radiomics[fold]\n",
    "\n",
    "    dev_predictions_resnet = []\n",
    "    for fold in all_dev_predictions_from_resnet.keys():\n",
    "        dev_predictions_resnet = dev_predictions_resnet + all_dev_predictions_from_resnet[fold]\n",
    "\n",
    "    X_fusion_dev = pd.DataFrame(\n",
    "        list(zip(dev_predictions_resnet, dev_predictions_resnet)),\n",
    "        columns=['radiomics', 'resnet'])\n",
    "    X_fusion_test = pd.DataFrame(\n",
    "        list(zip(test_estimated_radiomics_best_fold, test_estimated_resnet_best_fold)),\n",
    "        columns=['radiomics', 'resnet'])\n",
    "\n",
    "    fit_model_fusion = rfc_cv_fusion.fit(X_fusion_dev, Y_dev)\n",
    "    best_model_radiomics_fusion = fit_model_fusion.best_estimator_\n",
    "\n",
    "    # get predictions for fusion RFC model\n",
    "    X_fusion_test_predicted = best_model_radiomics_fusion.predict_proba(X_fusion_test)\n",
    "\n",
    "    dev_auc = fit_model_fusion.cv_results_['mean_train_score'][fit_model_fusion.best_index_]\n",
    "    test_auc = roc_auc_score(Y_test, X_fusion_test_predicted[:, 1])\n",
    "\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_dev_fusion_auc.csv\".format(model_results_name)), dev_auc)\n",
    "    np.savetxt(os.path.join(save_folder, \"{}_test_fusion_auc.csv\".format(model_results_name)), test_auc)\n",
    "\n",
    "    logging.info('Finished training.')\n",
    "    logging.info(f'Time elapsed: {round(time.time() - training_start_time, 3)} seconds.')\n",
    "\n",
    "    return dev_auc, test_auc\n",
    "\n",
    "\n",
    "import math\n",
    "\n",
    "########################################################################\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 = 10\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",
    "    k = 5\n",
    "\n",
    "    fusion_grid_parameters = {\n",
    "        'n_estimators': [5, 10, 25],\n",
    "        'random_state': [42, 100],\n",
    "        'criterion': [\"entropy\"],\n",
    "        'min_samples_leaf': [1, 2, 4],\n",
    "        'max_depth': [1, 2],\n",
    "        'max_features': [\"auto\", None],\n",
    "        'max_samples': [0.5, 0.75, 1]\n",
    "    }\n",
    "\n",
    "    radiomics_grid_parameters = {\n",
    "        'n_estimators': [25, 50, 100],\n",
    "        'random_state': [42, 100],\n",
    "        'criterion': [\"entropy\"],\n",
    "        'min_samples_leaf': [2, 4, 8],\n",
    "        'max_depth': [1, 2, 4],\n",
    "        'max_features': [\"auto\", None],\n",
    "        'max_samples': [0.5, 0.75, 1]\n",
    "    }\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",
    "    df_result = pd.DataFrame(columns=['Development AUC', 'Testing AUC'])\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",
    "        development_size = int(0.75 * len(dataset))\n",
    "        test_size = int(0.25 * len(dataset))\n",
    "        development_dataset, test_dataset = torch.utils.data.random_split(dataset, [development_size, test_size])\n",
    "        logging.info(f\"Datasplit -> Development: {development_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(math.floor(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(development_dataset, test_dataset,\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",
    "                            radiomics_grid_parameters=radiomics_grid_parameters,\n",
    "                            fusion_grid_parameters=fusion_grid_parameters)\n",
    "\n",
    "        dev_fusion_auc, test_fusion_auc = results\n",
    "\n",
    "        df_trial = pd.DataFrame(list(zip([dev_fusion_auc], [test_fusion_auc])),\n",
    "                                columns=['Development', 'Testing'])\n",
    "        df_result = pd.concat([df_result, df_trial], axis=0, ignore_index=True)\n",
    "\n",
    "        logging.info(\"\\r\",\n",
    "                     f\"Trial: {t + 1}/{num_trials}...Completed. Fusion AUC on development set: {dev_fusion_auc}, fusion AUC on test set: {test_fusion_auc}.\\n\")\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",
    "\n",
    "    curr_directory = os.path.dirname(os.path.realpath(__file__))\n",
    "    sub_directory = \"RFC_fusion_results\"\n",
    "    file_name = \"RFC_fusion_results_\" + time.strftime(\"%Y_%m_%d-%H_%M_%S\") + \".csv\"\n",
    "    filepath = os.path.join(curr_directory, sub_directory, file_name)\n",
    "    df_result.to_csv(filepath, index=False)\n",
    "\n",
    "logging.info('---------------------')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "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
}