{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"display_name":"pytorch_gpu","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.12"},"papermill":{"default_parameters":{},"duration":17947.044951,"end_time":"2024-07-21T19:59:57.358595","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-07-21T15:00:50.313644","version":"2.5.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport pydicom\nimport numpy as np\nimport os\nimport glob\nfrom tqdm import tqdm\nimport gc\n\nimport torchvision\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset\nfrom fastai.vision.all import *\nimport segmentation_models_pytorch as smp\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","metadata":{"execution":{"iopub.execute_input":"2024-07-21T15:00:53.042579Z","iopub.status.busy":"2024-07-21T15:00:53.042231Z","iopub.status.idle":"2024-07-21T15:01:24.461173Z","shell.execute_reply":"2024-07-21T15:01:24.460287Z"},"papermill":{"duration":31.433451,"end_time":"2024-07-21T15:01:24.463598","exception":false,"start_time":"2024-07-21T15:00:53.030147","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"SEED = 777\nFOLDS = [1,2,3,4,5]\nPATH = 'C:/Users/Angel/kaggle/'# Main path\nTRAIN_PATH = 'C:/Users/Angel/kaggle/train/'# Training images folder\nPATCH_SIZE = 128\nBS = 32\nEPOCHS = 10\nLmax = 15","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta_f = pd.read_csv(PATH + 'train_series_descriptions.csv')\ndf_meta_f.tail()","metadata":{"execution":{"iopub.execute_input":"2024-07-21T15:01:24.675045Z","iopub.status.busy":"2024-07-21T15:01:24.674395Z","iopub.status.idle":"2024-07-21T15:01:24.692944Z","shell.execute_reply":"2024-07-21T15:01:24.692058Z"},"papermill":{"duration":0.035373,"end_time":"2024-07-21T15:01:24.695008","exception":false,"start_time":"2024-07-21T15:01:24.659635","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_coor = pd.read_csv(PATH +'train_label_coordinates.csv')\ndf_coor.tail()","metadata":{"execution":{"iopub.execute_input":"2024-07-21T15:01:24.724580Z","iopub.status.busy":"2024-07-21T15:01:24.724234Z","iopub.status.idle":"2024-07-21T15:01:24.850435Z","shell.execute_reply":"2024-07-21T15:01:24.849512Z"},"papermill":{"duration":0.143427,"end_time":"2024-07-21T15:01:24.852417","exception":false,"start_time":"2024-07-21T15:01:24.708990","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"left_df_coor = df_coor[df_coor.condition == 'Left Neural Foraminal Narrowing'][['series_id','level','instance_number']].groupby(['series_id','level']).mean().sort_values(['series_id','level'])\nleft_df_coor","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"right_df_coor = df_coor[df_coor.condition == 'Right Neural Foraminal Narrowing'][['series_id','level','instance_number']].groupby(['series_id','level']).mean().sort_values(['series_id','level'])\nright_df_coor","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flipped_foraminal = left_df_coor.merge(right_df_coor,left_on=left_df_coor.index.to_numpy(),right_on=right_df_coor.index.to_numpy())\nflipped_foraminal.tail()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flipped_foraminal['flipped'] = flipped_foraminal.instance_number_x > flipped_foraminal.instance_number_y\nflipped_foraminal.tail()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flipped_foraminal['level'] = flipped_foraminal['key_0'].apply(lambda v: v[1])\nflipped_foraminal['series_id'] = flipped_foraminal['key_0'].apply(lambda v: v[0])\nflipped_foraminal_counts = flipped_foraminal[['series_id','flipped']].groupby('series_id').count().reset_index()\nprint(flipped_foraminal_counts.tail())\nflipped_foraminal = flipped_foraminal[['series_id','flipped']].groupby('series_id').sum().reset_index()\nprint(flipped_foraminal.tail())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flipped_foraminal['flipped'] = flipped_foraminal['flipped']/flipped_foraminal_counts['flipped']\nflipped_foraminal.tail()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flipped_foraminal.groupby('flipped').count()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flipped_foraminal = flipped_foraminal[flipped_foraminal['flipped'].isin([0.,1.])].reset_index(drop=True)\nflipped_foraminal['flipped'] = flipped_foraminal['flipped'] > 0\nflipped_foraminal.tail()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(PATH + 'train_split.csv')\ntrain_sagittal_direction = df_meta_f.merge(\n    flipped_foraminal,\n    left_on='series_id',\n    right_on='series_id'\n).merge(\n    train[['study_id','fold']],\n    left_on='study_id',\n    right_on='study_id'\n)\ntrain_sagittal_direction['flip'] = False\nftrain_sagittal_direction = train_sagittal_direction.copy()\nftrain_sagittal_direction[['flipped','flip']] = ~ftrain_sagittal_direction[['flipped','flip']]\ntrain_sagittal_direction = pd.concat([train_sagittal_direction,ftrain_sagittal_direction]).reset_index(drop=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sagittal_direction.groupby('fold').count()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sagittal_direction.groupby('series_description').count()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Sagittal_T1_direction_Dataset(Dataset):\n    def __init__(self, df, VALID=False, alpha=0):\n        self.data = df\n        self.VALID = VALID\n        self.alpha = alpha\n        self.resize = torchvision.transforms.Resize((PATCH_SIZE,PATCH_SIZE),antialias=True)\n\n    def __len__(self):\n        return len(self.data)\n\n    def __getitem__(self, index):\n        row = self.data.iloc[index]\n        \n        sample = TRAIN_PATH + str(row['study_id']) + '/' + str(row['series_id'])\n\n        images = [x.replace('\\\\','/') for x in glob.glob(sample+'/*.dcm')]\n        \n        images.sort(reverse=False, key=lambda x: int(x.split('/')[-1].replace('.dcm', '')))\n        \n        image = np.stack([pydicom.dcmread(x).pixel_array for x in images])\n        D,H,W = image.shape\n\n        slices_mask = torch.ones(Lmax,dtype=bool).to(device)\n        if D > Lmax:\n            slices = np.rint(torch.arange(Lmax)*D/Lmax).long()\n            slices_mask[:] = False\n            D = Lmax\n        else:\n            slices = np.rint(torch.arange(D)*Lmax/D).long()\n            slices_mask[slices] = False\n            slices = torch.arange(D)\n#       By plane resizing I've been distorting the proportions\n        if H > W:\n            d = W\n            if not self.VALID:\n                h = int((H - d)*(.5 + self.alpha*(.5 - np.random.rand())))\n            else:\n                h = (H - d)//2\n            image = image[slices,h:h+d]\n            H = W\n        elif H < W:\n            d = H\n            if not self.VALID:\n                w = int((W - d)*(.5 + self.alpha*(.5 - np.random.rand())))\n            else:\n                w = (W - d)//2\n            image = image[slices,:,w:w+d]\n            W = H\n        else:\n            image = image[slices]\n\n        image = self.resize(torch.as_tensor((image/np.max(image)).astype(np.float32)))\n        if row['flip']: image = image.flip(0)\n        \n        X = torch.zeros(Lmax,PATCH_SIZE,PATCH_SIZE).to(device).float()\n        X[~slices_mask] = image.to(device)\n\n        return [X,slices_mask],torch.as_tensor(row['flipped']).long().to(device)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CosineAnnealingAlpha\ndef nt(nmin,nmax,tcur,tmax):\n    return (nmax - .5*(nmax-nmin)*(1+np.cos(tcur*np.pi/tmax))).astype(np.float32)\n\nplt.plot(nt(.25,1,np.arange(EPOCHS),EPOCHS))\nplt.show()\n\n# callback to update alpha during training\ndef cb(self):\n    alpha = torch.as_tensor(nt(.25,1,learn.train_iter,EPOCHS*n_iter))\n    learn.dls.train_ds.alpha = alpha\nalpha_cb = Callback(before_batch=cb)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SinusoidalPosEmb(nn.Module):\n    def __init__(self, dim=16, M=10000):\n        super().__init__()\n        self.dim = dim\n        self.M = M\n\n    def forward(self, x):\n        device = x.device\n        half_dim = self.dim // 2\n        emb = math.log(self.M) / half_dim\n        emb = torch.exp(torch.arange(half_dim, device=device) * (-emb))\n        emb = x[...,None] * emb[None,...]\n        emb = torch.cat((emb.sin(), emb.cos()), dim=-1)\n        return emb\n\nclass Sagittal_T1_direction_ViT(nn.Module):\n    def __init__(self, dim=512, depth=24, head_size=64, **kwargs):\n        super().__init__()\n        CNN = torchvision.models.resnet18(weights='DEFAULT')\n        W = nn.Parameter(CNN.conv1.weight.sum(1, keepdim=True))\n        CNN.conv1 = nn.Conv2d(1, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)\n        CNN.conv1.weight = W\n        CNN.fc = nn.Identity()\n        self.emb = CNN.to(device)\n        self.slices_enc = nn.Parameter(SinusoidalPosEmb(dim)(torch.arange(Lmax, device=device).unsqueeze(0)))\n        self.slices_transformer = nn.TransformerEncoder(\n                nn.TransformerEncoderLayer(d_model=dim, nhead=dim//head_size, dim_feedforward=4*dim,\n                dropout=0.1, activation=nn.GELU(), batch_first=True, norm_first=True, device=device), depth)\n        self.proj_out = nn.Linear(dim,2).to(device)\n    \n    def forward(self, x):\n        x,slices_mask = x\n        slices_mask = slices_mask.view(-1,15)\n        x = self.emb(x.view(-1,1,PATCH_SIZE,PATCH_SIZE))\n        x = x.view(-1,15,512)\n        x = x + self.slices_enc\n        x = self.slices_transformer(x,src_key_padding_mask=slices_mask)\n        x[slices_mask] = 0\n        x = x.sum(1)/(~slices_mask).sum(1).unsqueeze(-1)\n        x = self.proj_out(x.view(-1,512))\n        return x","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def seed_everything(seed):\n    random.seed(seed)\n    np.random.seed(seed)\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = False","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = train_sagittal_direction\nfor f in FOLDS:\n    seed_everything(SEED)\n    model = Sagittal_T1_direction_ViT()\n    \n    tdf = df[df.fold != f]\n    vdf = df[df.fold == f]\n    tds = Sagittal_T1_direction_Dataset(tdf)\n    vds = Sagittal_T1_direction_Dataset(vdf,VALID=True)\n    tdl = torch.utils.data.DataLoader(tds, batch_size=BS, shuffle=True, drop_last=True)\n    vdl = torch.utils.data.DataLoader(vds, batch_size=BS, shuffle=False)\n    dls = DataLoaders(tdl,vdl)\n\n    n_iter = len(tds)//BS\n\n    learn = Learner(\n        dls,\n        model,\n        loss_func=nn.CrossEntropyLoss(),\n        cbs=[\n            ShowGraphCallback(),\n            GradientClip(3.0),\n            alpha_cb\n        ]\n    )\n    learn.fit_one_cycle(EPOCHS, lr_max=5e-4, wd=0.05, pct_start=0.02)\n    torch.save(model,'Sagittal_T1_direction_'+str(f))\n    del model,tdf,vdf,tds,vds,tdl,vdl,dls,learn\n    gc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = train_sagittal_direction\nCV_pred = []\nCV_true = []\nfor f in FOLDS:\n    model = torch.load('Sagittal_T1_direction_'+str(f))\n\n    vdf = df[df.fold == f]\n    vds = Sagittal_T1_direction_Dataset(vdf,VALID=True)\n    vdl = torch.utils.data.DataLoader(vds, batch_size=BS, shuffle=False)\n    y_pred = []\n    y_true = []\n    with torch.no_grad():\n        for X,Y in tqdm(vdl):\n            y_pred = y_pred + model(X).tolist()\n            y_true = y_true + Y.tolist()\n    print(nn.CrossEntropyLoss()(torch.as_tensor(y_pred),torch.as_tensor(y_true)))\n    CV_pred = CV_pred + y_pred\n    CV_true = CV_true + y_true\n    del model,vdf,vds,vdl,y_pred,y_true\n    gc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sklearn\n\nCV_pred = torch.as_tensor(CV_pred)\nCV_true = torch.as_tensor(CV_true)\nprint('CV Loss: ',nn.CrossEntropyLoss()(CV_pred,CV_true))\nprint('CV Confusion Matrix:')\nsklearn.metrics.confusion_matrix(CV_pred[:,1]>CV_pred[:,0], CV_true)","metadata":{},"execution_count":null,"outputs":[]}]}