{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Normalized Voxels: Align Planes, Adjust Contrast, and Crop\n\nAs shown in several notebooks, MRI plane type (Axial, Coronal, and Sagittal) is not consistent among patients or MRI scan types (FLAIR, T1w, T1wCE, T2w).\nWhile augmentations might alleviate this inconsistency, it is better to train models using MRI voxels that are consistent in terms of plane type.\nThis notebook shows we can obtain normalized voxels by appropriately rotating MRI voxels.\nI found that simply rotating MRI voxels is not enough because the order of planes is also inconsistent in some cases.\nFor example, even with the same Sagittal type, some of scans were in left-to-right order, while others were the other way around.\nAs Instance Number does not help, Image Position (Patient) is used in this notebook to reorder stacked images.\nAfter normalized voxels with respect to planes, contrast is adjusted and then voxels are cropped.\nFinally, the voxel is resized to arbitrary fixed size.\n\n## Normalized Voxel Datasets\n\nThe normalized voxels created the above procedure were stored as a dataset. Please refer to the second half of this notebook.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"# Because of the different scanning distances, i think 3dconv is not suitable.so use conv2d for voxel;\nfor mcr2 loss,label so little.\nmcr2: https://github.com/Ma-Lab-Berkeley/ReduNet\nbut when submit commit error: Notebook Threw Exception, seven time;","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport cv2\nimport pydicom\nimport matplotlib.pyplot as plt\n\nDATASET = 'train'\nscan_types = ['FLAIR','T1w','T1wCE','T2w']\ndata_root = Path(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification\")","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:56.98277Z","iopub.execute_input":"2021-09-13T12:42:56.983149Z","iopub.status.idle":"2021-09-13T12:42:56.990767Z","shell.execute_reply.started":"2021-09-13T12:42:56.983114Z","shell.execute_reply":"2021-09-13T12:42:56.987542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:56.992533Z","iopub.execute_input":"2021-09-13T12:42:56.993274Z","iopub.status.idle":"2021-09-13T12:42:57.000513Z","shell.execute_reply.started":"2021-09-13T12:42:56.993235Z","shell.execute_reply":"2021-09-13T12:42:56.999615Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_image_plane(data):\n    cords = [round(j) for j in data.ImageOrientationPatient]\n    return cords\n\ndef get_voxel(study_id, scan_type, split=\"train\"):\n    imgs = []\n    dcm_dir = data_root.joinpath( split , study_id, scan_type)\n    dcm_paths = sorted(dcm_dir.glob(\"*.dcm\"), key=lambda x: int(x.stem.split(\"-\")[-1]))\n    positions = []\n\n    for dcm_path in dcm_paths:\n        img = pydicom.dcmread(str(dcm_path))\n        imgs.append(img.pixel_array)\n        positions.append(img.ImagePositionPatient)\n\n    plane = get_image_plane(img)        \n    voxel = np.stack(imgs)\n\n    rotDir = []\n    rotDir.append(positions[-1][0]-positions[0][0])\n    rotDir.append(positions[-1][1]-positions[0][1])\n    rotDir.append(positions[-1][2]-positions[0][2])\n\n    rotDir = np.array(rotDir)\n    rotDir = rotDir / np.max(np.absolute(rotDir))\n    rotDir = np.around(rotDir)\n\n    rotVec = []\n\n    rotVec.append(np.arctan2(rotDir[0], rotDir[1]))\n    rotVec.append(np.arctan2(rotDir[1], rotDir[2]))\n    rotVec.append(np.arctan2(rotDir[2], rotDir[0]))\n\n    rotVec = np.array(rotVec)\n    rotVec = rotVec / np.max(np.absolute(rotVec))\n    rotVec = np.around(rotVec)\n\n    voxel = np.rot90(voxel, plane[1], (2, 0))\n    voxel = np.rot90(voxel, plane[2], (0, 1))\n    voxel = np.rot90(voxel, plane[3], (1, 2))\n    voxel = np.rot90(voxel, plane[5], (0, 1))\n\n    if plane[0] == 0:\n        voxel = np.flip(voxel, 1)\n    if plane[4] == 0:\n        voxel = np.flip(voxel, 0)\n    if rotDir[1] == 1:\n        voxel = np.flip(voxel, 1)\n    if rotDir[2] == -1:\n        voxel = np.flip(voxel, 0)\n\n    return voxel, plane","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.04053Z","iopub.execute_input":"2021-09-13T12:42:57.040764Z","iopub.status.idle":"2021-09-13T12:42:57.051397Z","shell.execute_reply.started":"2021-09-13T12:42:57.040742Z","shell.execute_reply":"2021-09-13T12:42:57.050472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_contrast(voxel):\n    if voxel.sum() == 0:\n        return voxel\n    voxel = voxel - np.min(voxel)\n    voxel = voxel / np.max(voxel)\n    voxel = (voxel * 255).astype(np.uint8)\n    return voxel","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.05288Z","iopub.execute_input":"2021-09-13T12:42:57.053403Z","iopub.status.idle":"2021-09-13T12:42:57.060673Z","shell.execute_reply.started":"2021-09-13T12:42:57.053366Z","shell.execute_reply":"2021-09-13T12:42:57.059782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def crop_voxel(voxel):\n    if voxel.sum() == 0:\n        return voxel\n    keep = (voxel.mean(axis=(0, 1)) > 0)\n    voxel = voxel[:, :, keep]\n    keep = (voxel.mean(axis=(0, 2)) > 0)\n    voxel = voxel[:, keep]\n    keep = (voxel.mean(axis=(1, 2)) > 0)\n    voxel = voxel[keep]\n    return voxel","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.063553Z","iopub.execute_input":"2021-09-13T12:42:57.063821Z","iopub.status.idle":"2021-09-13T12:42:57.070413Z","shell.execute_reply.started":"2021-09-13T12:42:57.063783Z","shell.execute_reply":"2021-09-13T12:42:57.069409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Sample planes along the longest axis and resize the sampled planes.\nBy sampling along the longest axis, the degradation due to sampling is minimized.\nThe best way is to resize twice (e.g. (x, y) axis then (y, z) axis) but it is computationally expensive.","metadata":{}},{"cell_type":"code","source":"def resize_orlay_voxel_forconv2d(voxel, sz=256, NUM_IMAGES=16):\n#     my_log.info(\"voxel.shape  no  0  resize is {}:\".format(voxel.shape))\n    NUM_IMAGES = sz\n    if np.argmin(voxel.shape) == 0:\n        if voxel.shape[0] < NUM_IMAGES:\n            output = np.zeros((NUM_IMAGES, sz, sz), dtype=np.uint8)\n            for i,s  in enumerate(np.linspace(0, voxel.shape[0] - 1, voxel.shape[0])):\n                b = int(abs(int(voxel.shape[0] / 2) - NUM_IMAGES / 2) + i)\n                output[b] = cv2.resize(voxel[i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 11  resize is {}:\".format(output.shape))\n        else:\n            output = np.zeros((voxel.shape[0], sz, sz), dtype=np.uint8)\n            for i, s in enumerate(np.linspace(0, voxel.shape[0] - 1, voxel.shape[0])):\n                output[i] = cv2.resize(voxel[i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 12  resize is {}:\".format(output.shape))\n            voxel = output\n            output = np.zeros((sz, sz, sz), dtype=np.uint8)\n            for i, s in enumerate(np.linspace(0, sz - 1, sz)):\n                output[:, i] = cv2.resize(voxel[:, i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 13  resize is {}:\".format(output.shape))\n        output = torch.tensor(output).permute(0, 1, 2)\n    elif np.argmin(voxel.shape) == 1:\n        if voxel.shape[1] < NUM_IMAGES:\n            output = np.zeros((sz, NUM_IMAGES, sz), dtype=np.uint8)\n            for i,s  in enumerate(np.linspace(0, voxel.shape[1] - 1, voxel.shape[1])):\n                b = int(abs(int(voxel.shape[1] / 2) - NUM_IMAGES / 2) + i)\n                output[:, b] = cv2.resize(voxel[:, i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 21  resize is {}:\".format(output.shape))\n        else:\n            output = np.zeros((sz, voxel.shape[1], sz), dtype=np.uint8)\n            for i, s in enumerate(np.linspace(0, voxel.shape[1] - 1, voxel.shape[1])):\n                output[:, i] = cv2.resize(voxel[:, i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 22  resize is {}:\".format(output.shape))\n            voxel = output\n            output = np.zeros((sz, sz, sz), dtype=np.uint8)\n            for i, s in enumerate(np.linspace(0, sz - 1, sz)):\n                output[i] = cv2.resize(voxel[i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 23  resize is {}:\".format(output.shape))\n        output = torch.tensor(output).permute(1, 0, 2)\n    elif np.argmin(voxel.shape) == 2:\n        if voxel.shape[2] < NUM_IMAGES:\n            output = np.zeros((sz, sz, NUM_IMAGES), dtype=np.uint8)\n            for i,s  in enumerate(np.linspace(0, voxel.shape[2] - 1, voxel.shape[2])):\n                b = int(abs(int(voxel.shape[2] / 2) - NUM_IMAGES / 2) + i)\n                output[:, :, b] = cv2.resize(voxel[:, :, i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 31  resize is {}:\".format(output.shape))\n        else:\n            output = np.zeros((sz, sz, voxel.shape[2]), dtype=np.uint8)\n            for i, s in enumerate(np.linspace(0, voxel.shape[2] - 1, voxel.shape[2])):\n                output[:, :, i] = cv2.resize(voxel[:, :, i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 32  resize is {}:\".format(output.shape))\n            voxel = output\n            output = np.zeros((sz, sz, sz), dtype=np.uint8)\n            for i, s in enumerate(np.linspace(0, sz - 1, sz)):\n                output[:, i] = cv2.resize(voxel[:, i], (sz, sz))\n#             my_log.info(\"voxel.shape  do 33  resize is {}:\".format(output.shape))\n        output = torch.tensor(output).permute(2, 0, 1)\n    return output\n","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.100452Z","iopub.execute_input":"2021-09-13T12:42:57.1007Z","iopub.status.idle":"2021-09-13T12:42:57.122545Z","shell.execute_reply.started":"2021-09-13T12:42:57.100678Z","shell.execute_reply":"2021-09-13T12:42:57.121674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resize_orlay_voxel_forconv2d000(voxel, sz=256, NUM_IMAGES=16):\n    output = np.zeros((sz, sz, sz), dtype=np.uint8)\n\n    if np.argmax(voxel.shape) == 0:\n        for i, s in enumerate(np.linspace(0, voxel.shape[0] - 1, sz)):\n            output[i] = cv2.resize(voxel[int(s)], (sz, sz))\n    elif np.argmax(voxel.shape) == 1:\n        for i, s in enumerate(np.linspace(0, voxel.shape[1] - 1, sz)):\n            output[:, i] = cv2.resize(voxel[:, int(s)], (sz, sz))\n    elif np.argmax(voxel.shape) == 2:\n        for i, s in enumerate(np.linspace(0, voxel.shape[2] - 1, sz)):\n            output[:, :, i] = cv2.resize(voxel[:, :, int(s)], (sz, sz))\n\n    return output","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.123988Z","iopub.execute_input":"2021-09-13T12:42:57.124602Z","iopub.status.idle":"2021-09-13T12:42:57.136148Z","shell.execute_reply.started":"2021-09-13T12:42:57.124567Z","shell.execute_reply":"2021-09-13T12:42:57.135266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import glob\nfrom torch.utils import data as torch_data\nfrom torch.utils.data import Dataset\nimport os\n\n#scan_types = ['FLAIR'] #,'T1w','T1wCE','T2w']\n\n#DataRetriever7voxelconv\nclass DataRetriever(torch_data.Dataset):\n    def __init__(self, paths, targets,  split='train', vsz = 96):\n        self.paths = paths\n        self.targets = targets\n        self.vsz = vsz\n        self.split = split\n        self.NUM_IMAGES = vsz\n        if not os.path.exists(os.path.join('./', self.split)):\n            os.mkdir(os.path.join('./', self.split))\n        os.system(\"ls -la ./train/\")  \n        #os.mkdir(os.path.join('./', test))\n\n    def __len__(self):\n        return len(self.paths)\n\n    def __getitem__(self, index):\n        _id = self.paths[index]\n        # _id = 388\n        #print(_id)\n        case = str(_id).zfill(5)\n        save_path = os.path.join('../input/newvox96', self.split, case + \".size.\" + str( self.vsz) + \".voxel4type.npy\")\n        #save_path = os.path.join('./', self.split, case + \".size.\" + str( self.vsz) + \".voxel4type.npy\")\n        if Path(save_path).exists():\n            voxels =  np.load(save_path)\n        else:\n            for i, scan_type in enumerate(scan_types):\n                voxel, plane = get_voxel(case, scan_type, self.split)\n                voxel = normalize_contrast(voxel)\n#                 print(\"voxel.shape  no  crop is :\", voxel.shape)\n                voxel = crop_voxel(voxel)\n#                 print(\"voxel.shape  do  crop is :\", voxel.shape)\n                voxel = resize_orlay_voxel_forconv2d(voxel, self.vsz, self.NUM_IMAGES)\n                if i > 0:\n                    voxel = np.concatenate((voxels, voxel), axis = 0)\n#                     print(voxel.shape)\n                voxels = voxel\n            np.save(save_path, voxels)\n        y = torch.tensor(self.targets[index], dtype=torch.int8)\n        #time.sleep(0.1)\n        return torch.tensor(voxels).float(), y\n    \nclass DataRetrievertest(torch_data.Dataset):\n    def __init__(self, paths, targets,  split='test', vsz = 96):\n        self.paths = paths\n        self.targets = targets\n        self.vsz = vsz\n        self.split = split\n        self.NUM_IMAGES = vsz\n        if not os.path.exists(os.path.join('./', self.split)):\n            os.mkdir(os.path.join('./', self.split))\n        #os.system(\"ls -la ./test/\")  checkpoints\n        \n\n\n    def __len__(self):\n        return len(self.paths)\n\n    def __getitem__(self, index):\n        _id = self.paths[index]\n        # _id = 388\n        #print(_id)\n        case = str(_id).zfill(5)\n        save_path = os.path.join('./', self.split, case + \".size.\" + str( self.vsz) + \".voxel4type.npy\")\n        #save_path = os.path.join('../input/voxel96', self.split, case + \".size.\" + str( self.vsz) + \".voxel4type.npy\")\n        if Path(save_path).exists():\n            voxels =  np.load(save_path)\n        else:\n            for i, scan_type in enumerate(scan_types):\n                voxel, plane = get_voxel(case, scan_type, self.split)\n                voxel = normalize_contrast(voxel)\n#                 print(\"voxel.shape  no  crop is :\", voxel.shape)\n                voxel = crop_voxel(voxel)\n#                 print(\"voxel.shape  do  crop is :\", voxel.shape)\n                voxel = resize_orlay_voxel_forconv2d(voxel, self.vsz, self.NUM_IMAGES)\n                if i > 0:\n                    voxel = np.concatenate((voxels, voxel), axis = 0)\n#                     print(voxel.shape)\n                voxels = voxel\n            np.save(save_path, voxels)\n#         y = torch.tensor(self.targets[index], dtype=torch.int8)\n        #time.sleep(0.1)\n        return torch.tensor(voxels).float(), _id","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.159769Z","iopub.execute_input":"2021-09-13T12:42:57.160049Z","iopub.status.idle":"2021-09-13T12:42:57.178797Z","shell.execute_reply.started":"2021-09-13T12:42:57.160025Z","shell.execute_reply":"2021-09-13T12:42:57.177924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Normalized Voxel Datasets\nThe normalized voxels created the above procedure were stored as a dataset:\n\n- [64x64x64 voxel](https://www.kaggle.com/ren4yu/rsna-miccai-voxel-64-dataset)\n- [128x128x128 voxel](https://www.kaggle.com/ren4yu/rsna-miccai-voxel-128-dataset)\n- [256x256x256 voxel](https://www.kaggle.com/ren4yu/rsna-miccai-voxel-256-dataset)\n\nThe directory structure is as follows:\n\n```\nvoxel\n├── train\n│   ├── 00000\n│   │   ├── FLAIR.npy\n│   │   ├── T1w.npy\n│   │   ├── T1wCE.npy\n│   │   └── T2w.npy\n│   ├── 00002\n│   │   ├── FLAIR.npy\n...\n├── test\n│   ├── 00001\n│   │   ├── FLAIR.npy\n│   │   ├── T1w.npy\n│   │   ├── T1wCE.npy\n│   │   └── T2w.npy\n│   ├── 00013\n│   │   ├── FLAIR.npy\n```\n\nSome voxels do not exist because all images belonging to these scans are completely black:\n\n- ('train', '00109', 'FLAIR.npy')\n- ('train', '00123', 'T1w.npy')\n- ('train', '00123', 'T2w.npy')\n- ('train', '00709', 'FLAIR.npy')\n\nLet's get one voxel and visualize it.","metadata":{}},{"cell_type":"code","source":"! ls ../input","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.180487Z","iopub.execute_input":"2021-09-13T12:42:57.180868Z","iopub.status.idle":"2021-09-13T12:42:57.854052Z","shell.execute_reply.started":"2021-09-13T12:42:57.180819Z","shell.execute_reply":"2021-09-13T12:42:57.853016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_voxel(study_id, scan_type=\"FLAIR\", split=\"train\", sz=256):\n    assert sz in (64, 128, 256)\n    data_root = Path(f\"../input/rsna-miccai-voxel-{sz}-dataset\")\n    npy_path = Path(data_root).joinpath(\"voxel\", split, study_id, f\"{scan_type}.npy\")\n    voxel = np.load(str(npy_path))\n    return voxel","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.856478Z","iopub.execute_input":"2021-09-13T12:42:57.856907Z","iopub.status.idle":"2021-09-13T12:42:57.864252Z","shell.execute_reply.started":"2021-09-13T12:42:57.856863Z","shell.execute_reply":"2021-09-13T12:42:57.863264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport logging\nimport json\nimport numpy as np\nimport torch\n\n\n\n\ndef save_params(model_dir, params):\n    \"\"\"Save params to a .json file. Params is a dictionary of parameters.\"\"\"\n    path = os.path.join(model_dir, 'params.json')\n    with open(path, 'w') as f:\n        json.dump(params, f, indent=2, sort_keys=True)\n\ndef update_params(model_dir, pretrain_dir):\n    \"\"\"Updates architecture and feature dimension from pretrain directory \n    to new directoy. \"\"\"\n    params = load_params(model_dir)\n    old_params = load_params(pretrain_dir)\n    params['arch'] = old_params[\"arch\"]\n    params['fd'] = old_params['fd']\n    save_params(model_dir, params)\n\ndef load_params(model_dir):\n    \"\"\"Load params.json file in model directory and return dictionary.\"\"\"\n    _path = os.path.join(model_dir, \"params.json\")\n    with open(_path, 'r') as f:\n        _dict = json.load(f)\n    return _dict\n\ndef save_state(model_dir, *entries, filename='losses.csv'):\n    \"\"\"Save entries to csv. Entries is list of numbers. \"\"\"\n    csv_path = os.path.join(model_dir, filename)\n    assert os.path.exists(csv_path), 'CSV file is missing in project directory.'\n    with open(csv_path, 'a') as f:\n        f.write('\\n'+','.join(map(str, entries)))\n\ndef save_ckpt(model_dir, net, epoch):\n    \"\"\"Save PyTorch checkpoint to ./checkpoints/ directory in model directory. \"\"\"\n    torch.save(net.state_dict(), os.path.join(model_dir, 'checkpoints', \n        'model-epoch{}.pt'.format(epoch)))\n\ndef save_labels(model_dir, labels, epoch):\n    \"\"\"Save labels of a certain epoch to directory. \"\"\"\n    path = os.path.join(model_dir, 'plabels', f'epoch{epoch}.npy')\n    np.save(path, labels)\n\ndef compute_accuracy0(y_pred, y_true):\n    \"\"\"Compute accuracy by counting correct classification. \"\"\"\n    assert y_pred.shape == y_true.shape\n    return 1 - np.count_nonzero(y_pred - y_true) / y_true.size\n\ndef compute_accuracy(y_pred, y_true):\n    \"\"\"Compute accuracy by counting correct classification. \"\"\"\n    assert y_pred.shape == y_true.shape\n    if type(y_pred) == torch.Tensor:\n        n_wrong = torch.count_nonzero(y_pred - y_true).item()\n    elif type(y_pred) == np.ndarray:\n        n_wrong = np.count_nonzero(y_pred - y_true)\n    else:\n        raise TypeError(\"Not Tensor nor Array type.\")\n    n_samples = len(y_pred)\n    return 1 - n_wrong / n_samples\n\n\ndef clustering_accuracy(labels_true, labels_pred):\n    \"\"\"Compute clustering accuracy.\"\"\"\n    from sklearn.metrics.cluster import supervised\n    from scipy.optimize import linear_sum_assignment\n    labels_true, labels_pred = supervised.check_clusterings(labels_true, labels_pred)\n    value = supervised.contingency_matrix(labels_true, labels_pred)\n    [r, c] = linear_sum_assignment(-value)\n    return value[r, c].sum() / len(labels_true)","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.865716Z","iopub.execute_input":"2021-09-13T12:42:57.866042Z","iopub.status.idle":"2021-09-13T12:42:57.883414Z","shell.execute_reply.started":"2021-09-13T12:42:57.866016Z","shell.execute_reply":"2021-09-13T12:42:57.882578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.linear_model import SGDClassifier\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.svm import LinearSVC,SVC\nimport time\nimport os\nfrom sklearn.calibration import CalibratedClassifierCV as CalibratedClassifierCV\n#from simple_log import get_log, add_file\n\n# my_log = get_log('acc_log')\n# timelog= time.time()\n# print(\"timelog is : \",timelog)\n# #add_file('acc_log', os.path.join(time.time() +  'acc_log.log'))\n# add_file('acc_log', os.path.join(str(time.time()) +  'acc_log.log'))\n# # my_log.info(\"acc is {},{}:\".format(,))\n\nfrom sklearn.datasets import load_iris\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.naive_bayes import GaussianNB\n\n\ndef bys( train_features, train_labels, test_features, test_labels,maxtest,subtest,test='train'):\n    gnb = GaussianNB()\n    gnb.fit(train_features, train_labels)\n    acc_train = gnb.score(train_features, train_labels)\n    acc_test = gnb.score(test_features, test_labels)\n    pred = []\n    testp = test\n    print(testp)\n    if testp == 'test':\n        pred = gnb.predict_proba(subtest)\n    else:\n        if acc_test > maxtest :\n            maxtest = acc_test\n        if acc_test > 0.51:\n            pred = gnb.predict_proba(subtest)\n    print(\"gnb acc is {},{},{}:\".format(maxtest,acc_train,acc_test))\n    return acc_train, acc_test ,maxtest, pred\n\n\n\n\n\ndef svm( train_features, train_labels, test_features, test_labels,maxtest,subtest,test='train'):\n    svm = SVC(verbose=0, random_state=10, max_iter=2100, probability=True)\n    svm.fit(train_features, train_labels)\n    acc_train = svm.score(train_features, train_labels)\n    acc_test = svm.score(test_features, test_labels)\n    pred = []\n    testp = test\n    print(testp)\n    if testp == 'test':\n        pred = svm.predict_proba(subtest)\n    else:\n        if acc_test > maxtest :\n            maxtest = acc_test\n        if acc_test > 0.51:\n            pred = svm.predict_proba(subtest)\n    print(\"svm acc is {},{},{}:\".format(maxtest,acc_train,acc_test))\n    return acc_train, acc_test ,maxtest, pred\n\ndef sgd( train_features, train_labels, test_features, test_labels,maxtest,subtest,test='train'):\n    clf = CalibratedClassifierCV(base_estimator=SGDClassifier(max_iter=5000, tol=1e-3, early_stopping=True, validation_fraction=0.8, n_iter_no_change=120 ))\n    clf.fit(train_features, train_labels)\n    acc_train = clf.score(train_features, train_labels)\n    acc_test = clf.score(test_features, test_labels)\n    pred = []\n    testp = test\n    print(testp)\n    if testp == 'test':\n        pred = clf.predict_proba(subtest)\n    else:\n        if acc_test > maxtest :\n            maxtest = acc_test\n        if acc_test > 0.51:\n            pred = clf.predict_proba(subtest)\n    print(\"sgd acc is {},{},{}:\".format(maxtest,acc_train,acc_test))\n    return acc_train, acc_test ,maxtest, pred","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.884727Z","iopub.execute_input":"2021-09-13T12:42:57.88512Z","iopub.status.idle":"2021-09-13T12:42:57.901971Z","shell.execute_reply.started":"2021-09-13T12:42:57.885085Z","shell.execute_reply":"2021-09-13T12:42:57.901003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\n\nclass BasicBlock(nn.Module):\n    expansion = 1\n\n    def __init__(self, in_planes, planes, stride=1):\n        super(BasicBlock, self).__init__()\n        self.conv1 = nn.Conv2d(in_planes, planes, kernel_size=3,\n                               stride=stride, padding=1, bias=False)\n        self.bn1 = nn.BatchNorm2d(planes)\n        self.conv2 = nn.Conv2d(planes, planes, kernel_size=3,\n                               stride=1, padding=1, bias=False)\n        self.bn2 = nn.BatchNorm2d(planes)\n\n        self.shortcut = nn.Sequential()\n        if stride != 1 or in_planes != self.expansion * planes:\n            self.shortcut = nn.Sequential(\n                nn.Conv2d(in_planes, self.expansion * planes,\n                          kernel_size=1, stride=stride, bias=False),\n                nn.BatchNorm2d(self.expansion * planes)\n            )\n\n    def forward(self, x):\n        out = F.relu(self.bn1(self.conv1(x)))\n        out = self.bn2(self.conv2(out))\n        out += self.shortcut(x)\n        out = F.relu(out)\n        return out\n\n\nclass Bottleneck(nn.Module):\n    expansion = 4\n\n    def __init__(self, in_planes, planes, stride=1):\n        super(Bottleneck, self).__init__()\n        self.conv1 = nn.Conv2d(in_planes, planes, kernel_size=1, bias=False)\n        self.bn1 = nn.BatchNorm2d(planes)\n        self.conv2 = nn.Conv2d(planes, planes, kernel_size=3,\n                               stride=stride, padding=1, bias=False)\n        self.bn2 = nn.BatchNorm2d(planes)\n        self.conv3 = nn.Conv2d(planes, self.expansion * planes,\n                               kernel_size=1, bias=False)\n        self.bn3 = nn.BatchNorm2d(self.expansion * planes)\n\n        self.shortcut = nn.Sequential()\n        if stride != 1 or in_planes != self.expansion * planes:\n            self.shortcut = nn.Sequential(\n                nn.Conv2d(in_planes, self.expansion * planes, kernel_size=1,\n                          stride=stride, bias=False),\n                nn.BatchNorm2d(self.expansion * planes)\n            )\n\n    def forward(self, x):\n        out = F.relu(self.bn1(self.conv1(x)))\n        out = F.relu(self.bn2(self.conv2(out)))\n        out = self.bn3(self.conv3(out))\n        out += self.shortcut(x)\n        out = F.relu(out)\n        return out\n\n\n\nclass ResNetControl(nn.Module):\n    def __init__(self, block, num_blocks, feature_dim=512):\n        super(ResNetControl, self).__init__()\n        self.in_planes = 256\n        self.feature_dim = feature_dim\n        self.conv1 = nn.Conv2d(384, 256, kernel_size=2, stride=1,\n                               padding=1, bias=False)\n        self.bn1 = nn.BatchNorm2d(256)\n        self.layer1 = self._make_layer(block, 256, num_blocks[0], stride=1)\n        self.layer2 = self._make_layer(block, 256, num_blocks[1], stride=2)\n        self.layer3 = self._make_layer(block, 256, num_blocks[2], stride=2)\n        self.layer4 = self._make_layer(block, feature_dim, num_blocks[3], stride=2)\n\n    def _make_layer(self, block, planes, num_blocks, stride):\n        strides = [stride] + [1] * (num_blocks - 1)\n        layers = []\n        for stride in strides:\n            layers.append(block(self.in_planes, planes, stride))\n            self.in_planes = planes * block.expansion\n        return nn.Sequential(*layers)\n\n    def forward(self, x):\n        out = F.relu(self.bn1(self.conv1(x)))\n        out = self.layer1(out)\n        out = self.layer2(out)\n        out = self.layer3(out)\n        out = self.layer4(out)\n        out = F.avg_pool2d(out, 4)\n        out = out.view(out.size(0), -1)\n        return F.normalize(out)\n\n\ndef ResNet18(feature_dim=512):\n    return ResNet(BasicBlock, [2, 2, 2, 2], feature_dim)\n\ndef ResNet18Control(feature_dim=512):\n    return ResNetControl(BasicBlock, [2, 2, 2, 2], feature_dim)\n","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.903505Z","iopub.execute_input":"2021-09-13T12:42:57.903866Z","iopub.status.idle":"2021-09-13T12:42:57.929251Z","shell.execute_reply.started":"2021-09-13T12:42:57.903816Z","shell.execute_reply":"2021-09-13T12:42:57.928398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nfrom tqdm import tqdm\n\nimport cv2\nimport numpy as np\nimport torch\nimport torch.nn\nimport torchvision\nimport torchvision.transforms as transforms\nfrom torch.utils.data import DataLoader\n\n# from cluster import ElasticNetSubspaceClustering, clustering_accuracy\n# import utils\n\n\n\n\n\ndef load_checkpoint(model_dir, epoch=None, eval_=False):\n    \"\"\"Load checkpoint from model directory. Checkpoints should be stored in \n    `model_dir/checkpoints/model-epochX.ckpt`, where `X` is the epoch number.\n    \n    Parameters:\n        model_dir (str): path to model directory\n        epoch (int): epoch number; set to None for last available epoch\n        eval_ (bool): PyTorch evaluation mode. set to True for testing\n        \n    Returns:\n        net (torch.nn.Module): PyTorch checkpoint at `epoch`\n        epoch (int): epoch number\n    \n    \"\"\"\n    if epoch is None: # get last epoch\n        ckpt_dir = os.path.join(model_dir, 'checkpoints')\n        epochs = [int(e[11:-3]) for e in os.listdir(ckpt_dir) if e[-3:] == \".pt\"]\n        epoch = np.sort(epochs)[-1]\n    ckpt_path = os.path.join(model_dir, 'checkpoints', 'model-epoch{}.pt'.format(epoch))\n    params = load_params(model_dir)\n    print('Loading checkpoint: {}'.format(ckpt_path))\n    state_dict = torch.load(ckpt_path)\n    net = load_architectures(params['arch'], params['fd'])\n    net.load_state_dict(state_dict)\n    del state_dict\n    if eval_:\n        net.eval()\n    return net, epoch\n\n    \ndef get_features(net, trainloader, verbose=True):\n    '''Extract all features out into one single batch. \n    \n    Parameters:\n        net (torch.nn.Module): get features using this model\n        trainloader (torchvision.dataloader): dataloader for loading data\n        verbose (bool): shows loading staus bar\n\n    Returns:\n        features (torch.tensor): with dimension (num_samples, feature_dimension)\n        labels (torch.tensor): with dimension (num_samples, )\n    '''\n    features = []\n    labels = []\n    if verbose:\n        train_bar = tqdm(trainloader, desc=\"extracting all features from dataset\")\n    else:\n        train_bar = trainloader\n    for step, (batch_imgs, batch_lbls) in enumerate(train_bar):\n        batch_features = net(batch_imgs.cuda())\n        features.append(batch_features.cpu().detach())\n        labels.append(batch_lbls)\n    return torch.cat(features), torch.cat(labels)\n    \n\ndef corrupt_labels(mode=\"default\"):\n    \"\"\"Returns higher corder function\"\"\"\n    if mode == \"default\":\n        from corrupt import default_corrupt\n        return default_corrupt\n    elif mode == \"asymmetric_noise\":\n        from corrupt import asymmetric_noise\n        return asymmetric_noise\n    elif mode == \"noisify_pairflip\":\n        from corrupt import noisify_pairflip\n        return noisify_pairflip\n    elif mode == \"noisify_multiclass_symmetric\":\n        from corrupt import noisify_multiclass_symmetric\n        return noisify_multiclass_symmetric\n\n\n\ndef label_to_membership(targets, num_classes=None):\n    \"\"\"Generate a true membership matrix, and assign value to current Pi.\n\n    Parameters:\n        targets (np.ndarray): matrix with one hot labels\n\n    Return:\n        Pi: membership matirx, shape (num_classes, num_samples, num_samples)\n\n    \"\"\"\n    targets = one_hot(targets, num_classes)\n    num_samples, num_classes = targets.shape\n    Pi = np.zeros(shape=(num_classes, num_samples, num_samples))\n    for j in range(len(targets)):\n        k = np.argmax(targets[j])\n        Pi[k, j, j] = 1.\n    return Pi\n\n\ndef membership_to_label(membership):\n    \"\"\"Turn a membership matrix into a list of labels.\"\"\"\n    _, num_classes, num_samples, _ = membership.shape\n    labels = np.zeros(num_samples)\n    for i in range(num_samples):\n        labels[i] = np.argmax(membership[:, i, i])\n    return labels\n\ndef one_hot(labels_int, n_classes):\n    \"\"\"Turn labels into one hot vector of K classes. \"\"\"\n    labels_onehot = torch.zeros(size=(len(labels_int), n_classes)).float()\n    for i, y in enumerate(labels_int):\n        labels_onehot[i, y] = 1.\n    return labels_onehot\n\n\n## Additional Augmentations\nclass GaussianBlur():\n    # Implements Gaussian blur as described in the SimCLR paper\n    def __init__(self, kernel_size, min=0.1, max=2.0):\n        self.min = min\n        self.max = max\n        # kernel size is set to be 10% of the image height/width\n        self.kernel_size = kernel_size\n\n    def __call__(self, sample):\n        sample = np.array(sample)\n\n        # blur the image with a 50% chance\n        prob = np.random.random_sample()\n\n        if prob < 0.5:\n            sigma = (self.max - self.min) * np.random.random_sample() + self.min\n            sample = cv2.GaussianBlur(sample, (self.kernel_size, self.kernel_size), sigma)\n\n        return sample\n\ndef sparse2coarse(targets):\n    \"\"\"CIFAR100 Coarse Labels. \"\"\"\n    coarse_targets = [ 4,  1, 14,  8,  0,  6,  7,  7, 18,  3,  3, 14,  9, 18,  7, 11,  3,\n                       9,  7, 11,  6, 11,  5, 10,  7,  6, 13, 15,  3, 15,  0, 11,  1, 10,\n                      12, 14, 16,  9, 11,  5,  5, 19,  8,  8, 15, 13, 14, 17, 18, 10, 16,\n                       4, 17,  4,  2,  0, 17,  4, 18, 17, 10,  3,  2, 12, 12, 16, 12,  1,\n                       9, 19,  2, 10,  0,  1, 16, 12,  9, 13, 15, 13, 16, 19,  2,  4,  6,\n                      19,  5,  5,  8, 19, 18,  1,  2, 15,  6,  0, 17,  8, 14, 13]\n    return np.array(coarse_targets)[targets]\n","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.932217Z","iopub.execute_input":"2021-09-13T12:42:57.932456Z","iopub.status.idle":"2021-09-13T12:42:57.956802Z","shell.execute_reply.started":"2021-09-13T12:42:57.932434Z","shell.execute_reply":"2021-09-13T12:42:57.956065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport torch\nfrom itertools import combinations\n\n\nclass MaximalCodingRateReduction(torch.nn.Module):\n    def __init__(self, gam1=1.0, gam2=1.0, eps=0.01):\n        super(MaximalCodingRateReduction, self).__init__()\n        self.gam1 = gam1\n        self.gam2 = gam2\n        self.eps = eps\n\n    def compute_discrimn_loss_empirical(self, W):\n        \"\"\"Empirical Discriminative Loss.\"\"\"\n        p, m = W.shape\n        I = torch.eye(p).cuda()\n        scalar = p / (m * self.eps)\n        logdet = torch.logdet(I + self.gam1 * scalar * W.matmul(W.T))\n        return logdet / 2.\n\n    def compute_compress_loss_empirical(self, W, Pi):\n        \"\"\"Empirical Compressive Loss.\"\"\"\n        p, m = W.shape\n        k, _, _ = Pi.shape\n        I = torch.eye(p).cuda()\n        compress_loss = 0.\n        for j in range(k):\n            trPi = torch.trace(Pi[j]) + 1e-8\n            scalar = p / (trPi * self.eps)\n            log_det = torch.logdet(I + scalar * W.matmul(Pi[j]).matmul(W.T))\n            compress_loss += log_det * trPi / m\n        return compress_loss / 2.\n\n    def compute_discrimn_loss_theoretical(self, W):\n        \"\"\"Theoretical Discriminative Loss.\"\"\"\n        p, m = W.shape\n        I = torch.eye(p).cuda()\n        scalar = p / (m * self.eps)\n        logdet = torch.logdet(I + scalar * W.matmul(W.T))\n        return logdet / 2.\n\n    def compute_compress_loss_theoretical(self, W, Pi):\n        \"\"\"Theoretical Compressive Loss.\"\"\"\n        p, m = W.shape\n        k, _, _ = Pi.shape\n        I = torch.eye(p).cuda()\n        compress_loss = 0.\n        for j in range(k):\n            trPi = torch.trace(Pi[j]) + 1e-8\n            scalar = p / (trPi * self.eps)\n            log_det = torch.logdet(I + scalar * W.matmul(Pi[j]).matmul(W.T))\n            compress_loss += trPi / (2 * m) * log_det\n        return compress_loss\n\n    def forward(self, X, Y, num_classes=None):\n        if num_classes is None:\n            num_classes = Y.max() + 1\n        W = X.T\n        Pi = label_to_membership(Y.numpy(), num_classes)\n        Pi = torch.tensor(Pi, dtype=torch.float32).cuda()\n\n        discrimn_loss_empi = self.compute_discrimn_loss_empirical(W)\n        compress_loss_empi = self.compute_compress_loss_empirical(W, Pi)\n        discrimn_loss_theo = self.compute_discrimn_loss_theoretical(W)\n        compress_loss_theo = self.compute_compress_loss_theoretical(W, Pi)\n \n        total_loss_empi = self.gam2 * -discrimn_loss_empi + compress_loss_empi\n        return (total_loss_empi,\n                [discrimn_loss_empi.item(), compress_loss_empi.item()],\n                [discrimn_loss_theo.item(), compress_loss_theo.item()])","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.958142Z","iopub.execute_input":"2021-09-13T12:42:57.958523Z","iopub.status.idle":"2021-09-13T12:42:57.976379Z","shell.execute_reply.started":"2021-09-13T12:42:57.958488Z","shell.execute_reply":"2021-09-13T12:42:57.975626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch.optim.lr_scheduler as lr_scheduler\nimport torch.optim as optim\nimport torch.optim.lr_scheduler as lr_scheduler\n\n\n\nnet = ResNet18Control(512)\n\n\ncriterion = MaximalCodingRateReduction(gam1=1, gam2=1, eps=0.001)\n#optimizer = optim.SGD(net.parameters(), lr=0.00001, momentum=0.9, weight_decay=5e-4)\noptimizer = torch.optim.Adam(net.parameters(),lr=0.000001)\nscheduler = lr_scheduler.MultiStepLR(optimizer, [200, 400, 600], gamma=0.1)\n# utils.save_params(model_dir, vars(args))\n\nfrom sklearn import model_selection as sk_model_selection\nfrom torch.utils import data as torch_data\nimport torch\nimport pandas as pd\n\n# import DataRetriever7voxelconv as DataRetriever\ndf = pd.read_csv(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv\")\n\nto_exclude = [109, 123, 123, 709] # [514, 11, 658, 537, 544, 688, 692, 694, 703, 709, 710, 711, 713, 719, 720, 722, 631, 109, 750, 621, 501,\n             #         118, 503, 630, 121, 123, 125, 126, 127, 551, 605, 539, 505]\ndf = df[~df['BraTS21ID'].isin(to_exclude)]\n\ndf_train, df_valid = sk_model_selection.train_test_split(df,  test_size=0.15 ,     random_state=42,     stratify=df[\"MGMT_value\"],  )\n# df_train = df_train.head(10)\n# df_valid = df_valid.head(10)\n\ntrain_data_retriever = DataRetriever(    df_train[\"BraTS21ID\"].values,       df_train[\"MGMT_value\"].values,  )\nvalid_data_retriever = DataRetriever(     df_valid[\"BraTS21ID\"].values,     df_valid[\"MGMT_value\"].values, )\ntrain_loader = torch_data.DataLoader(\n    train_data_retriever,\n    batch_size=16, #args.samples,  #100\n    shuffle=True,\n    num_workers=2,\n)\ntrain_loader2 = torch_data.DataLoader(\n    train_data_retriever,\n    batch_size=16, #200, #args.samples,  #100\n    shuffle=False,\n    num_workers=1,\n)\nvalid_loader = torch_data.DataLoader(\n    valid_data_retriever,\n    batch_size=8,\n    shuffle=False,\n    num_workers=1,\n)\n\nsubmission = pd.read_csv(f\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv\")\nsub_data_retriever = DataRetrievertest(    submission[\"BraTS21ID\"].values,       submission[\"MGMT_value\"].values, split='test',  )\nsub_loader = torch_data.DataLoader(\n    sub_data_retriever,\n    batch_size=8, #args.bs, #200, #args.samples,  #100\n    shuffle=False,\n    num_workers=2,\n)\n\ndf_train.head()","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:57.979087Z","iopub.execute_input":"2021-09-13T12:42:57.979508Z","iopub.status.idle":"2021-09-13T12:42:58.224912Z","shell.execute_reply.started":"2021-09-13T12:42:57.97947Z","shell.execute_reply":"2021-09-13T12:42:58.224047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"os.system('ls -la ./train/')\ntorch.cuda.empty_cache() \n# os.system('killall -9 python')\n# os.system('kill -9 32512')\n \n\nimport os\nfor dirname, _, filenames in os.walk('./checkpoints/'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:58.226601Z","iopub.execute_input":"2021-09-13T12:42:58.226964Z","iopub.status.idle":"2021-09-13T12:42:58.275822Z","shell.execute_reply.started":"2021-09-13T12:42:58.226925Z","shell.execute_reply":"2021-09-13T12:42:58.274876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# X_train, y_train = next(iter(train_loader))\n# print(\"X_train.shape, y_train.shape is :\", X_train.shape, y_train.shape)\n# X_train, y_train = X_train.to(device), y_train.to(device)\n\ndevice = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')\n      \nif os.path.exists(os.path.join('../', 'input/checkpt')):\n    maxtest=0\n    maxtestsave=0\n    ## Training\n    allsubpre=[]\n    \n    ckpt_dir = os.path.join('../', 'input/checkpt')\n    epochs = [int(e[11:-3]) for e in os.listdir(ckpt_dir) if e[-3:] == \".pt\"]\n    epoch = np.sort(epochs)[-1]\n    ckpt_path = os.path.join('../', 'input/checkpt', 'model-epoch{}.pt'.format(epoch))\n    print('Loading checkpoint: {}'.format(ckpt_path))\n    state_dict = torch.load(ckpt_path)\n    net = ResNet18Control(512)\n    net.to(device)\n    net.load_state_dict(state_dict)\n    del state_dict\n    net.eval()\n    \n    \n    train_features, train_labels = get_features(net, train_loader2)\n    test_features, test_labels = get_features(net, valid_loader)\n    sub_features, _sub_labels = get_features(net, sub_loader)\n    _,_,maxtest,pred = svm( train_features, train_labels, test_features, test_labels,maxtest,sub_features,'test')\n    if len(pred) >0  :\n        print()\n        allsubpre.append(pred)\n    _,_,maxtest,pred = sgd( train_features, train_labels, test_features, test_labels,maxtest,sub_features,'test')\n    if len(pred) >0  :\n        print()\n        allsubpre.append(pred)\n    _,_,maxtest,pred = bys( train_features, train_labels, test_features, test_labels,maxtest,sub_features,'test')\n    if len(pred) >0  :\n        print()\n        allsubpre.append(pred)\n        #print(allsubpre)\n    print(np.mean(allsubpre, axis=0))\n\n    print(submission.head())\n    submission = pd.DataFrame({\"BraTS21ID\": _sub_labels, \"MGMT_value\": np.mean(allsubpre, axis=0)[:,0]})\n    print(submission.head())\n    submission.to_csv(\"submission.csv\", index=False)\n    \n    \nelse:\n    if not os.path.exists(os.path.join('./', 'checkpoints')):\n        os.mkdir(os.path.join('./', 'checkpoints'))\n    #net = ResNet18Control(512)\n    net.to(device)\n    maxtest=0\n    maxtestsave=0\n    ## Training\n    allsubpre=[]\n    for epoch in range(25):\n        for step, (batch_imgs, batch_lbls) in enumerate(train_loader):\n            #print(batch_imgs.shape)\n            # features = net(batch_imgs)\n            features = net(batch_imgs.cuda())\n            loss, loss_empi, loss_theo = criterion(features, batch_lbls, num_classes=2) #trainset.num_classes)\n            if ( 'nan' in str(loss_empi[0]) or 'nan' in str(loss_empi[1]) )  or  ( 'nan' in str(loss_theo[0]) or 'nan' in str(loss_theo[1]) ) :\n                print(epoch, step, loss.item(), *loss_empi, *loss_theo)\n                print( \"because nan ,so  exit\")\n                exit()\n            optimizer.zero_grad()\n            loss.backward()\n            optimizer.step()\n            if step % 6 == 0:\n                print(epoch, step, loss.item(), *loss_empi, *loss_theo)\n        scheduler.step()\n\n        #print('train_features, train_labels start ')\n        if epoch % 1 == 0:\n\n            train_features, train_labels = get_features(net, train_loader2)\n            #print('train to cpu end ,start test feature  start ')\n            test_features, test_labels = get_features(net, valid_loader)\n    #        knn(args, train_features, train_labels, test_features, test_labels)\n            sub_features, _sub_labels = get_features(net, sub_loader)\n            _,_,maxtest,pred = svm( train_features, train_labels, test_features, test_labels,maxtest,sub_features,'train')\n            if len(pred) >0  :\n                print()\n                allsubpre.append(pred)\n                #utils.save_ckpt(model_dir, net, epoch)\n                #print(allsubpre)\n            _,_,maxtest,pred = sgd( train_features, train_labels, test_features, test_labels,maxtest,sub_features,'train')\n            if len(pred) >0  :\n                print()\n                allsubpre.append(pred)\n                #utils.save_ckpt(model_dir, net, epoch)\n                #print(allsubpre)\n            _,_,maxtest,pred = bys( train_features, train_labels, test_features, test_labels,maxtest,sub_features,'train')\n            if len(pred) >0  :\n                print()\n                allsubpre.append(pred)\n                #print(allsubpre)\n            if maxtest > maxtestsave:  # maxtestsave:\n                maxtestsave = maxtest\n                print('save ck ....')\n                save_ckpt('./', net, epoch)\n    print(np.mean(allsubpre, axis=0))\n\n    print(submission.head())\n    submission = pd.DataFrame({\"BraTS21ID\": _sub_labels, \"MGMT_value\": np.mean(allsubpre, axis=0)[:,0]})\n    print(submission.head())\n    submission.to_csv(\"submission.csv\", index=False)\n\n\n","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:42:58.277529Z","iopub.execute_input":"2021-09-13T12:42:58.277914Z","iopub.status.idle":"2021-09-13T12:43:40.200515Z","shell.execute_reply.started":"2021-09-13T12:42:58.277871Z","shell.execute_reply":"2021-09-13T12:43:40.199241Z"},"trusted":true},"execution_count":null,"outputs":[]}]}