{"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":"# RSNA MICCAI Brain Tumor Radiogenomic Classification\n\n[](http://google.com.au/)\n\nIn this notebook I will try to classify the images using differente EfficientNet models. To deal with 3D data I will try several method:\n* Using the maximum as its last layer\n* Adding and attention layer\n* Concatenating the neural network with a gaussian process\n\n## Importing necessary libraries","metadata":{}},{"cell_type":"code","source":"!pip install tensorflow-gpu==2.6 &> /dev/null\n!pip install gpflow &> /dev/null\n# Little hack so that it works on GPU\n!echo \"__version__ = '2.3.0'\" > /opt/conda/lib/python3.7/site-packages/gpflow/versions.py\n# Same but for CPU\n# !cp -r /opt/conda/lib/python3.7/site-packages/tensorboard-2.7.0.dist-info/ /opt/conda/lib/python3.7/site-packages/tensorboard-2.4.1.dist-info/\n!pip install gpflux &> /dev/null\n!echo \"__version__ = '0.2.4'\" > /opt/conda/lib/python3.7/site-packages/gpflux/version.py","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:38:32.261151Z","iopub.execute_input":"2021-11-08T12:38:32.261449Z","iopub.status.idle":"2021-11-08T12:41:25.597401Z","shell.execute_reply.started":"2021-11-08T12:38:32.261387Z","shell.execute_reply":"2021-11-08T12:41:25.596351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Basic libraries\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n# Neural network libraries\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras import layers\nfrom tensorflow.keras.callbacks import EarlyStopping\n\n# Reading images and creating video libraries\nimport cv2\nfrom IPython.display import HTML\nfrom base64 import b64encode\nimport matplotlib.animation as animation\nimport os\n\nimport SimpleITK as sitk","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-11-08T12:43:54.438521Z","iopub.execute_input":"2021-11-08T12:43:54.43889Z","iopub.status.idle":"2021-11-08T12:43:54.908787Z","shell.execute_reply.started":"2021-11-08T12:43:54.438857Z","shell.execute_reply":"2021-11-08T12:43:54.907926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Utility functions to visualize the images\n\nI display a video with a collection of the images of each folder.","metadata":{}},{"cell_type":"code","source":"def play(filename):\n html = ''\n video = open(filename,'rb').read()\n src = 'data:video/mp4;base64,' + b64encode(video).decode()\n html += '' % src \n return HTML(html)\n\ndef create_video(imgs, output='/kaggle/working/predicted.mp4', duration=30, subplot=True, \n frame_delay=200):\n fig, ax = plt.subplots(figsize=(15, 10))\n ims = []\n if not subplot:\n shape = imgs.shape[0]\n for i in range(duration):\n im = ax.imshow(imgs[i % shape], animated=True)\n ims.append([im])\n plt.close(fig)\n else:\n shapes = [imgs[views[0]].shape[0], imgs[views[1]].shape[0], \n imgs[views[2]].shape[0], imgs[views[3]].shape[0]]\n fig, ax = plt.subplots(2,2, figsize=(10,10))\n for k in range(duration):\n im_ = []\n for i in range(2):\n for j in range(2):\n im = ax[i,j].imshow(imgs[views[2*i+j]][k % shapes[2*i+j]], animated=True)\n im_.append(im)\n ax[i,j].set_title(views[2*i+j])\n plt.close()\n ims.append(im_)\n\n ani = animation.ArtistAnimation(fig, ims, interval=frame_delay, blit=True, repeat_delay=1000)\n\n ani.save(output)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:52.547797Z","iopub.execute_input":"2021-11-08T12:52:52.548133Z","iopub.status.idle":"2021-11-08T12:52:52.559335Z","shell.execute_reply.started":"2021-11-08T12:52:52.548102Z","shell.execute_reply":"2021-11-08T12:52:52.558488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Labels","metadata":{}},{"cell_type":"code","source":"target = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv')\npreds = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:53.061888Z","iopub.execute_input":"2021-11-08T12:52:53.062322Z","iopub.status.idle":"2021-11-08T12:52:53.093106Z","shell.execute_reply.started":"2021-11-08T12:52:53.062272Z","shell.execute_reply":"2021-11-08T12:52:53.092364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Read images utility function","metadata":{}},{"cell_type":"code","source":"# specify your image path\nviews = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\ndef load_imgs(idx, ignore_zeros=True, train=True):\n imgs = {}\n for view in views:\n save_ds = []\n if train:\n dir_path = os.walk(os.path.join(\n '../input/rsna-miccai-png/train/', idx, view\n ))\n else:\n dir_path = os.walk(os.path.join(\n '../input/rsna-miccai-png/test/', idx, view\n ))\n for path, subdirs, files in dir_path:\n for name in files:\n image_path = os.path.join(path, name) \n ds = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)\n save_ds.append(np.array(ds))\n if len(save_ds) == 0:\n save_ds = np.zeros((1,256,256))\n imgs[view] = np.array(save_ds)\n return imgs","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:54.165039Z","iopub.execute_input":"2021-11-08T12:52:54.165399Z","iopub.status.idle":"2021-11-08T12:52:54.172901Z","shell.execute_reply.started":"2021-11-08T12:52:54.165368Z","shell.execute_reply":"2021-11-08T12:52:54.172033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we try loading 32 images to see how much it takes. This will be the base to set the batch size later on so that each iteration is less expensive in time.","metadata":{}},{"cell_type":"code","source":"# %%time\n# for i in range(32):\n# idx = str(target.BraTS21ID[i]).zfill(5)\n# imgs = load_imgs(idx)","metadata":{"execution":{"iopub.status.busy":"2021-10-28T14:17:18.168139Z","iopub.execute_input":"2021-10-28T14:17:18.168598Z","iopub.status.idle":"2021-10-28T14:19:31.771205Z","shell.execute_reply.started":"2021-10-28T14:17:18.168553Z","shell.execute_reply":"2021-10-28T14:19:31.770292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Also, there are some folders without images. For those we simply define a zero-valued image so that the models work fine.","metadata":{}},{"cell_type":"code","source":"# Pathological one\nidx = str(109).zfill(5)\nimgs = load_imgs(idx)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Example of Image Visualization","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(2,2, figsize=(10,10))\nfor i in range(2):\n for j in range(2):\n m = ax[i,j].imshow(imgs[views[2*i+j]].mean(axis=0))\n ax[i,j].set_title(views[2*i+j])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-28T14:20:18.461128Z","iopub.execute_input":"2021-10-28T14:20:18.461449Z","iopub.status.idle":"2021-10-28T14:20:19.168249Z","shell.execute_reply.started":"2021-10-28T14:20:18.461421Z","shell.execute_reply":"2021-10-28T14:20:19.166628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Example of Video Visualization","metadata":{}},{"cell_type":"code","source":"create_video(imgs, duration=60, subplot=True, frame_delay=300)\nplay('predicted.mp4')","metadata":{"execution":{"iopub.status.busy":"2021-10-28T14:20:26.918198Z","iopub.execute_input":"2021-10-28T14:20:26.918551Z","iopub.status.idle":"2021-10-28T14:20:46.763541Z","shell.execute_reply.started":"2021-10-28T14:20:26.918522Z","shell.execute_reply":"2021-10-28T14:20:46.762653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### EfficientNet","metadata":{}},{"cell_type":"code","source":"!pip install efficientnet &> /dev/null","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:59.077823Z","iopub.execute_input":"2021-11-08T12:52:59.078144Z","iopub.status.idle":"2021-11-08T12:53:06.163896Z","shell.execute_reply.started":"2021-11-08T12:52:59.078112Z","shell.execute_reply":"2021-11-08T12:53:06.16283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## DataGenerator3D\n\nSince the data is massive we need to use a data generator. In the preprocessing we fetch all the images to the same dimensions and reduce bias using N4 Bias Field Correction. Apart from that the values of the images are zero mean and unit variance.\n\nNote: Since we need 3 channels I am just erasing the last one.","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\n\nclass DataGenerator3D(keras.utils.Sequence):\n 'Generates data for Keras'\n def __init__(self, list_IDs, labels=None, batch_size=256, dim=(512,512,512), n_channels=1,\n n_classes=2, shuffle=True, is_train=True):\n 'Initialization'\n self.dim = dim\n self.batch_size = batch_size\n self.labels = labels\n self.is_train = (labels is not None)\n self.list_IDs = list_IDs\n self.n_channels = n_channels\n self.n_classes = n_classes\n self.shuffle = shuffle\n self.on_epoch_end()\n\n def __len__(self):\n 'Denotes the number of batches per epoch'\n return int(np.floor(len(self.list_IDs) / self.batch_size))\n\n def __getitem__(self, index):\n 'Generate one batch of data'\n # Generate indexes of the batch\n list_IDs_temp = self.list_IDs[index*self.batch_size:(index+1)*self.batch_size]\n\n X = self.__data_generation(list_IDs_temp)\n # Generate data\n if self.is_train:\n y = self.labels[index*self.batch_size:(index+1)*self.batch_size]\n return np.array(X), np.array(y, dtype='float64')\n else:\n return np.array(X)\n\n def on_epoch_end(self):\n 'Updates indexes after each epoch'\n self.indexes = np.arange(len(self.list_IDs))\n if self.shuffle == True:\n np.random.shuffle(self.indexes)\n\n def __data_generation(self, list_IDs_temp):\n 'Generates data containing batch_size samples' # X : (n_samples, *dim, n_channels)\n # Initialization\n X = np.empty((self.batch_size, *self.dim, self.n_channels))\n\n # Generate data\n for i, ID in enumerate(list_IDs_temp):\n # Store sample\n idx = str(ID).zfill(5)\n imgs = load_imgs(idx, ignore_zeros=False, train=self.is_train)\n new_imgs = []\n for ii in range(2):\n for jj in range(2):\n if (ii == 1 and jj == 1):\n img_ = imgs[views[2*ii+jj]]\n img_ = np.array([cv2.resize(img_[k], dsize=(self.dim[1],self.dim[0]), interpolation=cv2.INTER_LINEAR) for k in range(img_.shape[0])])\n img_ = np.array([cv2.resize(img_.transpose(1,2,0)[k], dsize=(self.dim[2],self.dim[1]), interpolation=cv2.INTER_LINEAR) for k in range(self.dim[0])])\n\n # Removing radiofrequency inhomogeneity using N4 Bias Field Correction \n for p in range(len(img_)):\n inputImage = sitk.GetImageFromArray(img_[p])\n maskImage = sitk.GetImageFromArray((img_[p] >0.1) * 1)\n inputImage = sitk.Cast(inputImage, sitk.sitkFloat32)\n maskImage = sitk.Cast(maskImage, sitk.sitkUInt8)\n corrector = sitk.N4BiasFieldCorrectionImageFilter()\n numberFittingLevels = 4\n maxIter = 100\n if maxIter is not None:\n corrector.SetMaximumNumberOfIterations([maxIter]\n * numberFittingLevels)\n corrected_image = corrector.Execute(inputImage, maskImage)\n img_[p] = sitk.GetArrayFromImage(corrected_image)\n\n # Normalization\n sc = StandardScaler()\n img_ = np.array([sc.fit_transform(img_[i]) for i in range(img_.shape[0])])\n\n new_imgs.append(img_)\n new_imgs = np.concatenate(new_imgs).transpose(1,2,0).reshape((*self.dim,-1))\n X[i,] = new_imgs\n \n return X","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.16693Z","iopub.execute_input":"2021-11-08T12:53:06.167317Z","iopub.status.idle":"2021-11-08T12:53:06.69007Z","shell.execute_reply.started":"2021-11-08T12:53:06.167277Z","shell.execute_reply":"2021-11-08T12:53:06.689249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Usual train-validation split","metadata":{}},{"cell_type":"code","source":"target_red = target[288:(288+288)]","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.692028Z","iopub.execute_input":"2021-11-08T12:53:06.692431Z","iopub.status.idle":"2021-11-08T12:53:06.697145Z","shell.execute_reply.started":"2021-11-08T12:53:06.692392Z","shell.execute_reply":"2021-11-08T12:53:06.696137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train, X_val, y_train, y_val = train_test_split(target_red.BraTS21ID, target_red.MGMT_value,\n test_size=0.3, random_state=0,\n stratify=target_red.MGMT_value)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.699011Z","iopub.execute_input":"2021-11-08T12:53:06.6997Z","iopub.status.idle":"2021-11-08T12:53:06.760795Z","shell.execute_reply.started":"2021-11-08T12:53:06.699657Z","shell.execute_reply":"2021-11-08T12:53:06.760047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dim=(32,32,34)\nbatch_size=24\ntrain_dataset = DataGenerator3D(X_train, y_train, batch_size=batch_size, dim=dim)\nval_dataset = DataGenerator3D(X_val, y_val, batch_size=batch_size, dim=dim)\ntest_dataset = DataGenerator3D(preds.BraTS21ID, batch_size=batch_size, dim=dim)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.762015Z","iopub.execute_input":"2021-11-08T12:53:06.762373Z","iopub.status.idle":"2021-11-08T12:53:06.768914Z","shell.execute_reply.started":"2021-11-08T12:53:06.762338Z","shell.execute_reply":"2021-11-08T12:53:06.767856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #1\n\nApply efficientNet to each band of the image, then take the maximum.","metadata":{}},{"cell_type":"code","source":"import efficientnet.tfkeras as efn\n\nwith tf.device('/gpu:0'):\n def effNet(inp):\n inp = tf.reshape(tf.transpose(inp, perm=(0,3,1,2,4)), shape=(-1,dim[0], dim[1], 3))\n return efn.EfficientNetB0(include_top=False, pooling='avg')(inp)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.770586Z","iopub.execute_input":"2021-11-08T12:53:06.771049Z","iopub.status.idle":"2021-11-08T12:53:07.158719Z","shell.execute_reply.started":"2021-11-08T12:53:06.77101Z","shell.execute_reply":"2021-11-08T12:53:07.157838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n \n def build_model():\n inp = keras.Input(shape=(*dim,3))\n intermediate = effNet(inp)\n out = layers.Dense(1, activation='sigmoid')(intermediate)\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_0/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n model.compile(\n optimizer='adam', \n loss='binary_crossentropy',\n metrics=[keras.metrics.AUC()]\n )\n \n #history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=10, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:07.160067Z","iopub.execute_input":"2021-11-08T12:53:07.160435Z","iopub.status.idle":"2021-11-08T12:53:10.820455Z","shell.execute_reply.started":"2021-11-08T12:53:07.160397Z","shell.execute_reply":"2021-11-08T12:53:10.819618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #2\n\nThe maximum is substituted by an attention layer.","metadata":{}},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n D = 128\n \n def attention(inp):\n out = layers.Dense(D, activation='tanh')(inp)\n out = layers.Dense(1)(out)\n out = tf.reshape(out, shape=(-1,dim[2]))\n return tf.nn.softmax(out, axis=1)\n \n def build_model():\n inp = keras.Input(shape=(*dim,3))\n H = effNet(inp)\n A = attention(H)\n H = tf.reshape(H, shape=(-1,dim[2], H.shape[1]))\n A = tf.expand_dims(A, axis=1)\n intermediate = tf.linalg.matmul(A,H)\n intermediate = tf.squeeze(intermediate, axis=1)\n out = layers.Dense(1, activation='sigmoid')(intermediate)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_1/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_weights_only=True,\n verbose=1)\n\n model.compile(\n optimizer='adam', \n loss='binary_crossentropy',\n metrics=[keras.metrics.AUC()]\n )\n \n # history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=3, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:10.822378Z","iopub.execute_input":"2021-11-08T12:53:10.822638Z","iopub.status.idle":"2021-11-08T12:53:12.707818Z","shell.execute_reply.started":"2021-11-08T12:53:10.822613Z","shell.execute_reply":"2021-11-08T12:53:12.707034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #3\n\nIt's the model #1 plus a gaussian process. I am using the library GPFlux to do so.\n\nNote: Install gpflow before loading any library, otherwise it won't work","metadata":{}},{"cell_type":"markdown","source":"#### Toy Example","metadata":{}},{"cell_type":"code","source":"import gpflow\nimport gpflux\n\nfrom gpflow.config import default_float\n\n# Default float must be same type or it will give errors from time to time\ngpflow.config.set_default_float(\"float32\")\ntf.keras.backend.set_floatx(\"float32\")","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:58:14.843568Z","iopub.execute_input":"2021-11-08T12:58:14.843896Z","iopub.status.idle":"2021-11-08T12:58:14.852207Z","shell.execute_reply.started":"2021-11-08T12:58:14.843862Z","shell.execute_reply":"2021-11-08T12:58:14.851253Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = np.linspace(-10,10,1000)\nY = np.sin(X) + np.random.randn(1000) / 2","metadata":{"execution":{"iopub.status.busy":"2021-10-30T14:14:05.935435Z","iopub.execute_input":"2021-10-30T14:14:05.935794Z","iopub.status.idle":"2021-10-30T14:14:05.941823Z","shell.execute_reply.started":"2021-10-30T14:14:05.935757Z","shell.execute_reply":"2021-10-30T14:14:05.940744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = (X - X.mean()) / X.std()\nnum_data = X.shape\ninput_dim = 1","metadata":{"execution":{"iopub.status.busy":"2021-10-30T14:14:05.943316Z","iopub.execute_input":"2021-10-30T14:14:05.943705Z","iopub.status.idle":"2021-10-30T14:14:05.951043Z","shell.execute_reply.started":"2021-10-30T14:14:05.943668Z","shell.execute_reply":"2021-10-30T14:14:05.950082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_data = len(X)\nnum_inducing = 10\noutput_dim = 1\n\nkernel = gpflow.kernels.SquaredExponential()\ninducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,1))\n)\ngp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data, num_latent_gps=1\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:08:59.747485Z","iopub.execute_input":"2021-10-30T15:08:59.747835Z","iopub.status.idle":"2021-10-30T15:08:59.775127Z","shell.execute_reply.started":"2021-10-30T15:08:59.747805Z","shell.execute_reply":"2021-10-30T15:08:59.774156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"likelihood = gpflow.likelihoods.Gaussian(0.1)\n\n# So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\nlikelihood_container = gpflux.layers.TrackableLayer()\nlikelihood_container.likelihood = likelihood\n\ndef build_model():\n inp = keras.Input(shape=(input_dim))\n x = layers.Dense(100, activation=\"relu\")(inp)\n x = layers.Dense(100, activation=\"relu\")(x)\n x = layers.Dense(1, activation=\"linear\")(x)\n x = gp_layer.predict(x)\n # These two operations are to make the covariance matrix appear in the layers output later on\n x = tf.concat([x[0], x[1]], axis=0)\n x, y = tf.split(x, 2,axis=0)\n out = likelihood_container(x)\n return keras.Model(inputs=inp, outputs=out)\n\nmodel = build_model()\nloss = gpflux.losses.LikelihoodLoss(likelihood)","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:08:59.970291Z","iopub.execute_input":"2021-10-30T15:08:59.970559Z","iopub.status.idle":"2021-10-30T15:09:00.313313Z","shell.execute_reply.started":"2021-10-30T15:08:59.970533Z","shell.execute_reply":"2021-10-30T15:09:00.312552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.compile(loss=loss, optimizer=\"adam\")\nhist = model.fit(X, Y, epochs=100, verbose=0)\nplt.plot(hist.history[\"loss\"])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:09:02.527651Z","iopub.execute_input":"2021-10-30T15:09:02.527993Z","iopub.status.idle":"2021-10-30T15:09:08.361901Z","shell.execute_reply.started":"2021-10-30T15:09:02.527949Z","shell.execute_reply":"2021-10-30T15:09:08.361066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras import backend as K\ndef get_all_outputs(model, input_data, learning_phase=1):\n outputs = [layer.output for layer in model.layers[1:]] # exclude Input\n layers_fn = K.function([model.input], outputs)\n return layers_fn([input_data])\n\ndef get_layer_outputs(model, layer_name, input_data, learning_phase=1):\n outputs = [layer.output for layer in model.layers if layer_name in layer.name]\n layers_fn = K.function([model.input], outputs)\n return layers_fn([input_data])","metadata":{"execution":{"iopub.status.busy":"2021-10-30T14:23:17.36626Z","iopub.execute_input":"2021-10-30T14:23:17.366595Z","iopub.status.idle":"2021-10-30T14:23:17.372655Z","shell.execute_reply.started":"2021-10-30T14:23:17.366566Z","shell.execute_reply":"2021-10-30T14:23:17.37171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot(model, X, Y, ax=None):\n if ax is None:\n fig, ax = plt.subplots()\n\n x_margin = 0.2\n N_test = 100\n X_test = np.linspace(X.min() - x_margin, X.max() + x_margin, N_test).reshape(-1, 1)\n f_distribution = model(X_test)\n\n mean = f_distribution.numpy().squeeze()\n var = np.max(get_layer_outputs(model, 'tf.linalg.adjoint_', X_test)[0],axis=1).squeeze() + model.layers[-1].likelihood.variance.numpy()\n X_test = X_test.squeeze()\n lower = mean - 2 * np.sqrt(var)\n upper = mean + 2 * np.sqrt(var)\n\n ax.set_ylim(Y.min() - 0.5, Y.max() + 0.5)\n ax.plot(X, Y, \"kx\", alpha=0.5)\n ax.plot(X_test, mean, \"C1\")\n\n ax.fill_between(X_test, lower, upper, color=\"C1\", alpha=0.3)\n\n\n# plot(model, X, Y)","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:09:12.757453Z","iopub.execute_input":"2021-10-30T15:09:12.757773Z","iopub.status.idle":"2021-10-30T15:09:12.983286Z","shell.execute_reply.started":"2021-10-30T15:09:12.757745Z","shell.execute_reply":"2021-10-30T15:09:12.982478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Real data","metadata":{}},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n num_data = 4\n num_inducing = 4\n output_dim = 1\n\n kernel = gpflow.kernels.SquaredExponential()\n inducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,1))\n )\n gp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data * dim[2], num_latent_gps=output_dim\n )\n \n likelihood = gpflow.likelihoods.Bernoulli()\n\n # So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\n likelihood_container = gpflux.layers.TrackableLayer()\n likelihood_container.likelihood = likelihood\n\n def build_model():\n inp = keras.Input(shape=(*dim,3))\n latent = effNet(inp)\n latent = layers.Dense(1, activation='sigmoid')(latent)\n out = gp_layer.predict(latent)\n out = out[0]\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n out = likelihood_container(out)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_2/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n loss = gpflux.losses.LikelihoodLoss(likelihood)\n model.compile(\n optimizer='adam', \n loss=loss,\n metrics=[keras.metrics.AUC()]\n )\n \n #history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=10, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:59:29.000518Z","iopub.execute_input":"2021-10-30T15:59:29.000836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #4\n\nGaussian process but using clustering to select the inducing points.","metadata":{}},{"cell_type":"code","source":"class ClusteringLayer(layers.Layer):\n \"\"\"\n Clustering layer converts input sample (feature) to soft label.\n\n # Example\n ```\n model.add(ClusteringLayer(n_clusters=10))\n ```\n # Arguments\n n_clusters: number of clusters.\n weights: list of Numpy array with shape `(n_clusters, n_features)` witch represents the initial cluster centers.\n alpha: degrees of freedom parameter in Student's t-distribution. Default to 1.0.\n # Input shape\n 2D tensor with shape: `(n_samples, n_features)`.\n # Output shape\n 2D tensor with shape: `(n_samples, n_clusters)`.\n \"\"\"\n\n def __init__(self, n_clusters, weights=None, alpha=1.0, **kwargs):\n if 'input_shape' not in kwargs and 'input_dim' in kwargs:\n kwargs['input_shape'] = (kwargs.pop('input_dim'),)\n super(ClusteringLayer, self).__init__(**kwargs)\n self.n_clusters = n_clusters\n self.alpha = alpha\n self.initial_weights = weights\n self.input_spec = layers.InputSpec(ndim=2)\n\n def build(self, input_shape):\n assert len(input_shape) == 2\n input_dim = input_shape[1]\n self.input_spec = layers.InputSpec(dtype=K.floatx(), shape=(None, input_dim))\n self.clusters = self.add_weight(shape=(self.n_clusters, input_dim), initializer='glorot_uniform', name='clusters')\n if self.initial_weights is not None:\n self.set_weights(self.initial_weights)\n del self.initial_weights\n self.built = True\n\n def call(self, inputs, **kwargs):\n \"\"\" student t-distribution, as same as used in t-SNE algorithm. \n q_ij = 1/(1+dist(x_i, µ_j)^2), then normalize it.\n q_ij can be interpreted as the probability of assigning sample i to cluster j.\n (i.e., a soft assignment)\n Arguments:\n inputs: the variable containing data, shape=(n_samples, n_features)\n Return:\n q: student's t-distribution, or soft labels for each sample. shape=(n_samples, n_clusters)\n \"\"\"\n q = 1.0 / (1.0 + (K.sum(K.square(K.expand_dims(inputs, axis=1) - self.clusters), axis=2) / self.alpha))\n q **= (self.alpha + 1.0) / 2.0\n q = K.transpose(K.transpose(q) / K.sum(q, axis=1)) # Make sure each sample's 10 values add up to 1.\n return q\n\n def compute_output_shape(self, input_shape):\n assert input_shape and len(input_shape) == 2\n return input_shape[0], self.n_clusters\n\n def get_config(self):\n config = {'n_clusters': self.n_clusters}\n base_config = super(ClusteringLayer, self).get_config()\n return dict(list(base_config.items()) + list(config.items()))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n num_data = 4\n num_inducing = num_data\n output_dim = 1\n\n kernel = gpflow.kernels.SquaredExponential()\n inducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,1))\n )\n gp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data * dim[2], num_latent_gps=output_dim\n )\n \n likelihood = gpflow.likelihoods.Bernoulli()\n\n # So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\n likelihood_container = gpflux.layers.TrackableLayer()\n likelihood_container.likelihood = likelihood\n\n def build_model():\n inp = keras.Input(shape=(*dim,3))\n latent = effNet(inp)\n latent = ClusteringLayer(2, name='clustering')(latent)\n latent = layers.Dense(1, activation='sigmoid')(latent)\n out = gp_layer.predict(latent)\n out = out[0]\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n out = likelihood_container(out)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_2/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n loss = gpflux.losses.LikelihoodLoss(likelihood)\n model.compile(\n optimizer='adam', \n loss=loss,\n metrics=[keras.metrics.AUC()]\n )\n \n #history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=10, callbacks=[earlyStopping, cp_callback])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Simpler model\n\nLet's try with a simplified version of VGG for the CNN.","metadata":{}},{"cell_type":"code","source":"def conv(inp, n, n_out):\n x = inp\n for i in range(n):\n x = layers.BatchNormalization()(x)\n x = layers.Conv2D(\n n_out, (5,5), padding='same'\n )(x)\n x = layers.Dropout(0.2)(x)\n x = layers.Activation('relu')(x)\n return x\n\ndef pool(inp):\n return layers.MaxPool2D(\n pool_size=(2, 2), strides=None, padding='valid'\n )(inp)\n\ndef dense(inp, n_neur, drop=0.2, act='relu'):\n x = layers.BatchNormalization()(inp)\n x = layers.Dense(n_neur)(x)\n x = layers.Dropout(drop)(x)\n return layers.Activation(act)(x)\n\ndef VGG(inp):\n inp = tf.reshape(tf.transpose(inp, perm=(0,3,1,2,4)), shape=(-1,dim[0], dim[1], 1))\n x = conv(inp, 2, 16)\n x = pool(x)\n \n x = conv(x, 3, 32)\n x = pool(x)\n \n x = layers.Flatten()(x)\n return dense(x, 1, 0.4, act='linear')","metadata":{"execution":{"iopub.status.busy":"2021-11-04T19:19:35.521508Z","iopub.execute_input":"2021-11-04T19:19:35.521863Z","iopub.status.idle":"2021-11-04T19:19:35.534117Z","shell.execute_reply.started":"2021-11-04T19:19:35.521829Z","shell.execute_reply":"2021-11-04T19:19:35.533158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n num_data = batch_size\n num_inducing = 10\n output_dim = 1\n\n kernel = gpflow.kernels.SquaredExponential()\n inducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,output_dim))\n )\n gp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data * dim[2], num_latent_gps=output_dim\n )\n \n likelihood = gpflow.likelihoods.Bernoulli()\n\n # So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\n likelihood_container = gpflux.layers.TrackableLayer()\n likelihood_container.likelihood = likelihood","metadata":{"execution":{"iopub.status.busy":"2021-11-04T19:36:35.362913Z","iopub.execute_input":"2021-11-04T19:36:35.363248Z","iopub.status.idle":"2021-11-04T19:36:35.395375Z","shell.execute_reply.started":"2021-11-04T19:36:35.363219Z","shell.execute_reply":"2021-11-04T19:36:35.394377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n\n def build_model():\n inp = keras.Input(shape=(*dim,3))\n latent = VGG(inp)\n out = gp_layer.predict(latent)\n out = out[0]\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n out = layers.Activation('sigmoid')(out)\n out = likelihood_container(out)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_4/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n loss = gpflux.losses.LikelihoodLoss(likelihood)\n model.compile(\n optimizer='adam', \n loss=loss,\n metrics=[keras.metrics.AUC(), 'accuracy']\n )\n \n history = model.fit(train_dataset, validation_data=val_dataset,\n epochs=30, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-11-04T19:36:37.714636Z","iopub.execute_input":"2021-11-04T19:36:37.71502Z","iopub.status.idle":"2021-11-04T20:28:07.494106Z","shell.execute_reply.started":"2021-11-04T19:36:37.71499Z","shell.execute_reply":"2021-11-04T20:28:07.49322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}
![](\"https://storage.googleapis.com/kaggle-competitions/kaggle/29653/logos/header.png?t=2021-07-07-17-26-56\")
](http://google.com.au/)\n\nIn this notebook I will try to classify the images using differente EfficientNet models. To deal with 3D data I will try several method:\n* Using the maximum as its last layer\n* Adding and attention layer\n* Concatenating the neural network with a gaussian process\n\n## Importing necessary libraries","metadata":{}},{"cell_type":"code","source":"!pip install tensorflow-gpu==2.6 &> /dev/null\n!pip install gpflow &> /dev/null\n# Little hack so that it works on GPU\n!echo \"__version__ = '2.3.0'\" > /opt/conda/lib/python3.7/site-packages/gpflow/versions.py\n# Same but for CPU\n# !cp -r /opt/conda/lib/python3.7/site-packages/tensorboard-2.7.0.dist-info/ /opt/conda/lib/python3.7/site-packages/tensorboard-2.4.1.dist-info/\n!pip install gpflux &> /dev/null\n!echo \"__version__ = '0.2.4'\" > /opt/conda/lib/python3.7/site-packages/gpflux/version.py","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:38:32.261151Z","iopub.execute_input":"2021-11-08T12:38:32.261449Z","iopub.status.idle":"2021-11-08T12:41:25.597401Z","shell.execute_reply.started":"2021-11-08T12:38:32.261387Z","shell.execute_reply":"2021-11-08T12:41:25.596351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Basic libraries\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n# Neural network libraries\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras import layers\nfrom tensorflow.keras.callbacks import EarlyStopping\n\n# Reading images and creating video libraries\nimport cv2\nfrom IPython.display import HTML\nfrom base64 import b64encode\nimport matplotlib.animation as animation\nimport os\n\nimport SimpleITK as sitk","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-11-08T12:43:54.438521Z","iopub.execute_input":"2021-11-08T12:43:54.43889Z","iopub.status.idle":"2021-11-08T12:43:54.908787Z","shell.execute_reply.started":"2021-11-08T12:43:54.438857Z","shell.execute_reply":"2021-11-08T12:43:54.907926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Utility functions to visualize the images\n\nI display a video with a collection of the images of each folder.","metadata":{}},{"cell_type":"code","source":"def play(filename):\n html = ''\n video = open(filename,'rb').read()\n src = 'data:video/mp4;base64,' + b64encode(video).decode()\n html += '' % src \n return HTML(html)\n\ndef create_video(imgs, output='/kaggle/working/predicted.mp4', duration=30, subplot=True, \n frame_delay=200):\n fig, ax = plt.subplots(figsize=(15, 10))\n ims = []\n if not subplot:\n shape = imgs.shape[0]\n for i in range(duration):\n im = ax.imshow(imgs[i % shape], animated=True)\n ims.append([im])\n plt.close(fig)\n else:\n shapes = [imgs[views[0]].shape[0], imgs[views[1]].shape[0], \n imgs[views[2]].shape[0], imgs[views[3]].shape[0]]\n fig, ax = plt.subplots(2,2, figsize=(10,10))\n for k in range(duration):\n im_ = []\n for i in range(2):\n for j in range(2):\n im = ax[i,j].imshow(imgs[views[2*i+j]][k % shapes[2*i+j]], animated=True)\n im_.append(im)\n ax[i,j].set_title(views[2*i+j])\n plt.close()\n ims.append(im_)\n\n ani = animation.ArtistAnimation(fig, ims, interval=frame_delay, blit=True, repeat_delay=1000)\n\n ani.save(output)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:52.547797Z","iopub.execute_input":"2021-11-08T12:52:52.548133Z","iopub.status.idle":"2021-11-08T12:52:52.559335Z","shell.execute_reply.started":"2021-11-08T12:52:52.548102Z","shell.execute_reply":"2021-11-08T12:52:52.558488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Labels","metadata":{}},{"cell_type":"code","source":"target = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv')\npreds = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:53.061888Z","iopub.execute_input":"2021-11-08T12:52:53.062322Z","iopub.status.idle":"2021-11-08T12:52:53.093106Z","shell.execute_reply.started":"2021-11-08T12:52:53.062272Z","shell.execute_reply":"2021-11-08T12:52:53.092364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Read images utility function","metadata":{}},{"cell_type":"code","source":"# specify your image path\nviews = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\ndef load_imgs(idx, ignore_zeros=True, train=True):\n imgs = {}\n for view in views:\n save_ds = []\n if train:\n dir_path = os.walk(os.path.join(\n '../input/rsna-miccai-png/train/', idx, view\n ))\n else:\n dir_path = os.walk(os.path.join(\n '../input/rsna-miccai-png/test/', idx, view\n ))\n for path, subdirs, files in dir_path:\n for name in files:\n image_path = os.path.join(path, name) \n ds = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)\n save_ds.append(np.array(ds))\n if len(save_ds) == 0:\n save_ds = np.zeros((1,256,256))\n imgs[view] = np.array(save_ds)\n return imgs","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:54.165039Z","iopub.execute_input":"2021-11-08T12:52:54.165399Z","iopub.status.idle":"2021-11-08T12:52:54.172901Z","shell.execute_reply.started":"2021-11-08T12:52:54.165368Z","shell.execute_reply":"2021-11-08T12:52:54.172033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we try loading 32 images to see how much it takes. This will be the base to set the batch size later on so that each iteration is less expensive in time.","metadata":{}},{"cell_type":"code","source":"# %%time\n# for i in range(32):\n# idx = str(target.BraTS21ID[i]).zfill(5)\n# imgs = load_imgs(idx)","metadata":{"execution":{"iopub.status.busy":"2021-10-28T14:17:18.168139Z","iopub.execute_input":"2021-10-28T14:17:18.168598Z","iopub.status.idle":"2021-10-28T14:19:31.771205Z","shell.execute_reply.started":"2021-10-28T14:17:18.168553Z","shell.execute_reply":"2021-10-28T14:19:31.770292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Also, there are some folders without images. For those we simply define a zero-valued image so that the models work fine.","metadata":{}},{"cell_type":"code","source":"# Pathological one\nidx = str(109).zfill(5)\nimgs = load_imgs(idx)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Example of Image Visualization","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(2,2, figsize=(10,10))\nfor i in range(2):\n for j in range(2):\n m = ax[i,j].imshow(imgs[views[2*i+j]].mean(axis=0))\n ax[i,j].set_title(views[2*i+j])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-28T14:20:18.461128Z","iopub.execute_input":"2021-10-28T14:20:18.461449Z","iopub.status.idle":"2021-10-28T14:20:19.168249Z","shell.execute_reply.started":"2021-10-28T14:20:18.461421Z","shell.execute_reply":"2021-10-28T14:20:19.166628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Example of Video Visualization","metadata":{}},{"cell_type":"code","source":"create_video(imgs, duration=60, subplot=True, frame_delay=300)\nplay('predicted.mp4')","metadata":{"execution":{"iopub.status.busy":"2021-10-28T14:20:26.918198Z","iopub.execute_input":"2021-10-28T14:20:26.918551Z","iopub.status.idle":"2021-10-28T14:20:46.763541Z","shell.execute_reply.started":"2021-10-28T14:20:26.918522Z","shell.execute_reply":"2021-10-28T14:20:46.762653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### EfficientNet","metadata":{}},{"cell_type":"code","source":"!pip install efficientnet &> /dev/null","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:52:59.077823Z","iopub.execute_input":"2021-11-08T12:52:59.078144Z","iopub.status.idle":"2021-11-08T12:53:06.163896Z","shell.execute_reply.started":"2021-11-08T12:52:59.078112Z","shell.execute_reply":"2021-11-08T12:53:06.16283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## DataGenerator3D\n\nSince the data is massive we need to use a data generator. In the preprocessing we fetch all the images to the same dimensions and reduce bias using N4 Bias Field Correction. Apart from that the values of the images are zero mean and unit variance.\n\nNote: Since we need 3 channels I am just erasing the last one.","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\n\nclass DataGenerator3D(keras.utils.Sequence):\n 'Generates data for Keras'\n def __init__(self, list_IDs, labels=None, batch_size=256, dim=(512,512,512), n_channels=1,\n n_classes=2, shuffle=True, is_train=True):\n 'Initialization'\n self.dim = dim\n self.batch_size = batch_size\n self.labels = labels\n self.is_train = (labels is not None)\n self.list_IDs = list_IDs\n self.n_channels = n_channels\n self.n_classes = n_classes\n self.shuffle = shuffle\n self.on_epoch_end()\n\n def __len__(self):\n 'Denotes the number of batches per epoch'\n return int(np.floor(len(self.list_IDs) / self.batch_size))\n\n def __getitem__(self, index):\n 'Generate one batch of data'\n # Generate indexes of the batch\n list_IDs_temp = self.list_IDs[index*self.batch_size:(index+1)*self.batch_size]\n\n X = self.__data_generation(list_IDs_temp)\n # Generate data\n if self.is_train:\n y = self.labels[index*self.batch_size:(index+1)*self.batch_size]\n return np.array(X), np.array(y, dtype='float64')\n else:\n return np.array(X)\n\n def on_epoch_end(self):\n 'Updates indexes after each epoch'\n self.indexes = np.arange(len(self.list_IDs))\n if self.shuffle == True:\n np.random.shuffle(self.indexes)\n\n def __data_generation(self, list_IDs_temp):\n 'Generates data containing batch_size samples' # X : (n_samples, *dim, n_channels)\n # Initialization\n X = np.empty((self.batch_size, *self.dim, self.n_channels))\n\n # Generate data\n for i, ID in enumerate(list_IDs_temp):\n # Store sample\n idx = str(ID).zfill(5)\n imgs = load_imgs(idx, ignore_zeros=False, train=self.is_train)\n new_imgs = []\n for ii in range(2):\n for jj in range(2):\n if (ii == 1 and jj == 1):\n img_ = imgs[views[2*ii+jj]]\n img_ = np.array([cv2.resize(img_[k], dsize=(self.dim[1],self.dim[0]), interpolation=cv2.INTER_LINEAR) for k in range(img_.shape[0])])\n img_ = np.array([cv2.resize(img_.transpose(1,2,0)[k], dsize=(self.dim[2],self.dim[1]), interpolation=cv2.INTER_LINEAR) for k in range(self.dim[0])])\n\n # Removing radiofrequency inhomogeneity using N4 Bias Field Correction \n for p in range(len(img_)):\n inputImage = sitk.GetImageFromArray(img_[p])\n maskImage = sitk.GetImageFromArray((img_[p] >0.1) * 1)\n inputImage = sitk.Cast(inputImage, sitk.sitkFloat32)\n maskImage = sitk.Cast(maskImage, sitk.sitkUInt8)\n corrector = sitk.N4BiasFieldCorrectionImageFilter()\n numberFittingLevels = 4\n maxIter = 100\n if maxIter is not None:\n corrector.SetMaximumNumberOfIterations([maxIter]\n * numberFittingLevels)\n corrected_image = corrector.Execute(inputImage, maskImage)\n img_[p] = sitk.GetArrayFromImage(corrected_image)\n\n # Normalization\n sc = StandardScaler()\n img_ = np.array([sc.fit_transform(img_[i]) for i in range(img_.shape[0])])\n\n new_imgs.append(img_)\n new_imgs = np.concatenate(new_imgs).transpose(1,2,0).reshape((*self.dim,-1))\n X[i,] = new_imgs\n \n return X","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.16693Z","iopub.execute_input":"2021-11-08T12:53:06.167317Z","iopub.status.idle":"2021-11-08T12:53:06.69007Z","shell.execute_reply.started":"2021-11-08T12:53:06.167277Z","shell.execute_reply":"2021-11-08T12:53:06.689249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Usual train-validation split","metadata":{}},{"cell_type":"code","source":"target_red = target[288:(288+288)]","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.692028Z","iopub.execute_input":"2021-11-08T12:53:06.692431Z","iopub.status.idle":"2021-11-08T12:53:06.697145Z","shell.execute_reply.started":"2021-11-08T12:53:06.692392Z","shell.execute_reply":"2021-11-08T12:53:06.696137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train, X_val, y_train, y_val = train_test_split(target_red.BraTS21ID, target_red.MGMT_value,\n test_size=0.3, random_state=0,\n stratify=target_red.MGMT_value)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.699011Z","iopub.execute_input":"2021-11-08T12:53:06.6997Z","iopub.status.idle":"2021-11-08T12:53:06.760795Z","shell.execute_reply.started":"2021-11-08T12:53:06.699657Z","shell.execute_reply":"2021-11-08T12:53:06.760047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dim=(32,32,34)\nbatch_size=24\ntrain_dataset = DataGenerator3D(X_train, y_train, batch_size=batch_size, dim=dim)\nval_dataset = DataGenerator3D(X_val, y_val, batch_size=batch_size, dim=dim)\ntest_dataset = DataGenerator3D(preds.BraTS21ID, batch_size=batch_size, dim=dim)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.762015Z","iopub.execute_input":"2021-11-08T12:53:06.762373Z","iopub.status.idle":"2021-11-08T12:53:06.768914Z","shell.execute_reply.started":"2021-11-08T12:53:06.762338Z","shell.execute_reply":"2021-11-08T12:53:06.767856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #1\n\nApply efficientNet to each band of the image, then take the maximum.","metadata":{}},{"cell_type":"code","source":"import efficientnet.tfkeras as efn\n\nwith tf.device('/gpu:0'):\n def effNet(inp):\n inp = tf.reshape(tf.transpose(inp, perm=(0,3,1,2,4)), shape=(-1,dim[0], dim[1], 3))\n return efn.EfficientNetB0(include_top=False, pooling='avg')(inp)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:06.770586Z","iopub.execute_input":"2021-11-08T12:53:06.771049Z","iopub.status.idle":"2021-11-08T12:53:07.158719Z","shell.execute_reply.started":"2021-11-08T12:53:06.77101Z","shell.execute_reply":"2021-11-08T12:53:07.157838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n \n def build_model():\n inp = keras.Input(shape=(*dim,3))\n intermediate = effNet(inp)\n out = layers.Dense(1, activation='sigmoid')(intermediate)\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_0/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n model.compile(\n optimizer='adam', \n loss='binary_crossentropy',\n metrics=[keras.metrics.AUC()]\n )\n \n #history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=10, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:07.160067Z","iopub.execute_input":"2021-11-08T12:53:07.160435Z","iopub.status.idle":"2021-11-08T12:53:10.820455Z","shell.execute_reply.started":"2021-11-08T12:53:07.160397Z","shell.execute_reply":"2021-11-08T12:53:10.819618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #2\n\nThe maximum is substituted by an attention layer.","metadata":{}},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n D = 128\n \n def attention(inp):\n out = layers.Dense(D, activation='tanh')(inp)\n out = layers.Dense(1)(out)\n out = tf.reshape(out, shape=(-1,dim[2]))\n return tf.nn.softmax(out, axis=1)\n \n def build_model():\n inp = keras.Input(shape=(*dim,3))\n H = effNet(inp)\n A = attention(H)\n H = tf.reshape(H, shape=(-1,dim[2], H.shape[1]))\n A = tf.expand_dims(A, axis=1)\n intermediate = tf.linalg.matmul(A,H)\n intermediate = tf.squeeze(intermediate, axis=1)\n out = layers.Dense(1, activation='sigmoid')(intermediate)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_1/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_weights_only=True,\n verbose=1)\n\n model.compile(\n optimizer='adam', \n loss='binary_crossentropy',\n metrics=[keras.metrics.AUC()]\n )\n \n # history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=3, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:53:10.822378Z","iopub.execute_input":"2021-11-08T12:53:10.822638Z","iopub.status.idle":"2021-11-08T12:53:12.707818Z","shell.execute_reply.started":"2021-11-08T12:53:10.822613Z","shell.execute_reply":"2021-11-08T12:53:12.707034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #3\n\nIt's the model #1 plus a gaussian process. I am using the library GPFlux to do so.\n\nNote: Install gpflow before loading any library, otherwise it won't work","metadata":{}},{"cell_type":"markdown","source":"#### Toy Example","metadata":{}},{"cell_type":"code","source":"import gpflow\nimport gpflux\n\nfrom gpflow.config import default_float\n\n# Default float must be same type or it will give errors from time to time\ngpflow.config.set_default_float(\"float32\")\ntf.keras.backend.set_floatx(\"float32\")","metadata":{"execution":{"iopub.status.busy":"2021-11-08T12:58:14.843568Z","iopub.execute_input":"2021-11-08T12:58:14.843896Z","iopub.status.idle":"2021-11-08T12:58:14.852207Z","shell.execute_reply.started":"2021-11-08T12:58:14.843862Z","shell.execute_reply":"2021-11-08T12:58:14.851253Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = np.linspace(-10,10,1000)\nY = np.sin(X) + np.random.randn(1000) / 2","metadata":{"execution":{"iopub.status.busy":"2021-10-30T14:14:05.935435Z","iopub.execute_input":"2021-10-30T14:14:05.935794Z","iopub.status.idle":"2021-10-30T14:14:05.941823Z","shell.execute_reply.started":"2021-10-30T14:14:05.935757Z","shell.execute_reply":"2021-10-30T14:14:05.940744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = (X - X.mean()) / X.std()\nnum_data = X.shape\ninput_dim = 1","metadata":{"execution":{"iopub.status.busy":"2021-10-30T14:14:05.943316Z","iopub.execute_input":"2021-10-30T14:14:05.943705Z","iopub.status.idle":"2021-10-30T14:14:05.951043Z","shell.execute_reply.started":"2021-10-30T14:14:05.943668Z","shell.execute_reply":"2021-10-30T14:14:05.950082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_data = len(X)\nnum_inducing = 10\noutput_dim = 1\n\nkernel = gpflow.kernels.SquaredExponential()\ninducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,1))\n)\ngp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data, num_latent_gps=1\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:08:59.747485Z","iopub.execute_input":"2021-10-30T15:08:59.747835Z","iopub.status.idle":"2021-10-30T15:08:59.775127Z","shell.execute_reply.started":"2021-10-30T15:08:59.747805Z","shell.execute_reply":"2021-10-30T15:08:59.774156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"likelihood = gpflow.likelihoods.Gaussian(0.1)\n\n# So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\nlikelihood_container = gpflux.layers.TrackableLayer()\nlikelihood_container.likelihood = likelihood\n\ndef build_model():\n inp = keras.Input(shape=(input_dim))\n x = layers.Dense(100, activation=\"relu\")(inp)\n x = layers.Dense(100, activation=\"relu\")(x)\n x = layers.Dense(1, activation=\"linear\")(x)\n x = gp_layer.predict(x)\n # These two operations are to make the covariance matrix appear in the layers output later on\n x = tf.concat([x[0], x[1]], axis=0)\n x, y = tf.split(x, 2,axis=0)\n out = likelihood_container(x)\n return keras.Model(inputs=inp, outputs=out)\n\nmodel = build_model()\nloss = gpflux.losses.LikelihoodLoss(likelihood)","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:08:59.970291Z","iopub.execute_input":"2021-10-30T15:08:59.970559Z","iopub.status.idle":"2021-10-30T15:09:00.313313Z","shell.execute_reply.started":"2021-10-30T15:08:59.970533Z","shell.execute_reply":"2021-10-30T15:09:00.312552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.compile(loss=loss, optimizer=\"adam\")\nhist = model.fit(X, Y, epochs=100, verbose=0)\nplt.plot(hist.history[\"loss\"])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:09:02.527651Z","iopub.execute_input":"2021-10-30T15:09:02.527993Z","iopub.status.idle":"2021-10-30T15:09:08.361901Z","shell.execute_reply.started":"2021-10-30T15:09:02.527949Z","shell.execute_reply":"2021-10-30T15:09:08.361066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras import backend as K\ndef get_all_outputs(model, input_data, learning_phase=1):\n outputs = [layer.output for layer in model.layers[1:]] # exclude Input\n layers_fn = K.function([model.input], outputs)\n return layers_fn([input_data])\n\ndef get_layer_outputs(model, layer_name, input_data, learning_phase=1):\n outputs = [layer.output for layer in model.layers if layer_name in layer.name]\n layers_fn = K.function([model.input], outputs)\n return layers_fn([input_data])","metadata":{"execution":{"iopub.status.busy":"2021-10-30T14:23:17.36626Z","iopub.execute_input":"2021-10-30T14:23:17.366595Z","iopub.status.idle":"2021-10-30T14:23:17.372655Z","shell.execute_reply.started":"2021-10-30T14:23:17.366566Z","shell.execute_reply":"2021-10-30T14:23:17.37171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot(model, X, Y, ax=None):\n if ax is None:\n fig, ax = plt.subplots()\n\n x_margin = 0.2\n N_test = 100\n X_test = np.linspace(X.min() - x_margin, X.max() + x_margin, N_test).reshape(-1, 1)\n f_distribution = model(X_test)\n\n mean = f_distribution.numpy().squeeze()\n var = np.max(get_layer_outputs(model, 'tf.linalg.adjoint_', X_test)[0],axis=1).squeeze() + model.layers[-1].likelihood.variance.numpy()\n X_test = X_test.squeeze()\n lower = mean - 2 * np.sqrt(var)\n upper = mean + 2 * np.sqrt(var)\n\n ax.set_ylim(Y.min() - 0.5, Y.max() + 0.5)\n ax.plot(X, Y, \"kx\", alpha=0.5)\n ax.plot(X_test, mean, \"C1\")\n\n ax.fill_between(X_test, lower, upper, color=\"C1\", alpha=0.3)\n\n\n# plot(model, X, Y)","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:09:12.757453Z","iopub.execute_input":"2021-10-30T15:09:12.757773Z","iopub.status.idle":"2021-10-30T15:09:12.983286Z","shell.execute_reply.started":"2021-10-30T15:09:12.757745Z","shell.execute_reply":"2021-10-30T15:09:12.982478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Real data","metadata":{}},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n num_data = 4\n num_inducing = 4\n output_dim = 1\n\n kernel = gpflow.kernels.SquaredExponential()\n inducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,1))\n )\n gp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data * dim[2], num_latent_gps=output_dim\n )\n \n likelihood = gpflow.likelihoods.Bernoulli()\n\n # So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\n likelihood_container = gpflux.layers.TrackableLayer()\n likelihood_container.likelihood = likelihood\n\n def build_model():\n inp = keras.Input(shape=(*dim,3))\n latent = effNet(inp)\n latent = layers.Dense(1, activation='sigmoid')(latent)\n out = gp_layer.predict(latent)\n out = out[0]\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n out = likelihood_container(out)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_2/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n loss = gpflux.losses.LikelihoodLoss(likelihood)\n model.compile(\n optimizer='adam', \n loss=loss,\n metrics=[keras.metrics.AUC()]\n )\n \n #history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=10, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-10-30T15:59:29.000518Z","iopub.execute_input":"2021-10-30T15:59:29.000836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model #4\n\nGaussian process but using clustering to select the inducing points.","metadata":{}},{"cell_type":"code","source":"class ClusteringLayer(layers.Layer):\n \"\"\"\n Clustering layer converts input sample (feature) to soft label.\n\n # Example\n ```\n model.add(ClusteringLayer(n_clusters=10))\n ```\n # Arguments\n n_clusters: number of clusters.\n weights: list of Numpy array with shape `(n_clusters, n_features)` witch represents the initial cluster centers.\n alpha: degrees of freedom parameter in Student's t-distribution. Default to 1.0.\n # Input shape\n 2D tensor with shape: `(n_samples, n_features)`.\n # Output shape\n 2D tensor with shape: `(n_samples, n_clusters)`.\n \"\"\"\n\n def __init__(self, n_clusters, weights=None, alpha=1.0, **kwargs):\n if 'input_shape' not in kwargs and 'input_dim' in kwargs:\n kwargs['input_shape'] = (kwargs.pop('input_dim'),)\n super(ClusteringLayer, self).__init__(**kwargs)\n self.n_clusters = n_clusters\n self.alpha = alpha\n self.initial_weights = weights\n self.input_spec = layers.InputSpec(ndim=2)\n\n def build(self, input_shape):\n assert len(input_shape) == 2\n input_dim = input_shape[1]\n self.input_spec = layers.InputSpec(dtype=K.floatx(), shape=(None, input_dim))\n self.clusters = self.add_weight(shape=(self.n_clusters, input_dim), initializer='glorot_uniform', name='clusters')\n if self.initial_weights is not None:\n self.set_weights(self.initial_weights)\n del self.initial_weights\n self.built = True\n\n def call(self, inputs, **kwargs):\n \"\"\" student t-distribution, as same as used in t-SNE algorithm. \n q_ij = 1/(1+dist(x_i, µ_j)^2), then normalize it.\n q_ij can be interpreted as the probability of assigning sample i to cluster j.\n (i.e., a soft assignment)\n Arguments:\n inputs: the variable containing data, shape=(n_samples, n_features)\n Return:\n q: student's t-distribution, or soft labels for each sample. shape=(n_samples, n_clusters)\n \"\"\"\n q = 1.0 / (1.0 + (K.sum(K.square(K.expand_dims(inputs, axis=1) - self.clusters), axis=2) / self.alpha))\n q **= (self.alpha + 1.0) / 2.0\n q = K.transpose(K.transpose(q) / K.sum(q, axis=1)) # Make sure each sample's 10 values add up to 1.\n return q\n\n def compute_output_shape(self, input_shape):\n assert input_shape and len(input_shape) == 2\n return input_shape[0], self.n_clusters\n\n def get_config(self):\n config = {'n_clusters': self.n_clusters}\n base_config = super(ClusteringLayer, self).get_config()\n return dict(list(base_config.items()) + list(config.items()))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n num_data = 4\n num_inducing = num_data\n output_dim = 1\n\n kernel = gpflow.kernels.SquaredExponential()\n inducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,1))\n )\n gp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data * dim[2], num_latent_gps=output_dim\n )\n \n likelihood = gpflow.likelihoods.Bernoulli()\n\n # So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\n likelihood_container = gpflux.layers.TrackableLayer()\n likelihood_container.likelihood = likelihood\n\n def build_model():\n inp = keras.Input(shape=(*dim,3))\n latent = effNet(inp)\n latent = ClusteringLayer(2, name='clustering')(latent)\n latent = layers.Dense(1, activation='sigmoid')(latent)\n out = gp_layer.predict(latent)\n out = out[0]\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n out = likelihood_container(out)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_2/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n loss = gpflux.losses.LikelihoodLoss(likelihood)\n model.compile(\n optimizer='adam', \n loss=loss,\n metrics=[keras.metrics.AUC()]\n )\n \n #history = model.fit(train_dataset, validation_data=val_dataset,\n # epochs=10, callbacks=[earlyStopping, cp_callback])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Simpler model\n\nLet's try with a simplified version of VGG for the CNN.","metadata":{}},{"cell_type":"code","source":"def conv(inp, n, n_out):\n x = inp\n for i in range(n):\n x = layers.BatchNormalization()(x)\n x = layers.Conv2D(\n n_out, (5,5), padding='same'\n )(x)\n x = layers.Dropout(0.2)(x)\n x = layers.Activation('relu')(x)\n return x\n\ndef pool(inp):\n return layers.MaxPool2D(\n pool_size=(2, 2), strides=None, padding='valid'\n )(inp)\n\ndef dense(inp, n_neur, drop=0.2, act='relu'):\n x = layers.BatchNormalization()(inp)\n x = layers.Dense(n_neur)(x)\n x = layers.Dropout(drop)(x)\n return layers.Activation(act)(x)\n\ndef VGG(inp):\n inp = tf.reshape(tf.transpose(inp, perm=(0,3,1,2,4)), shape=(-1,dim[0], dim[1], 1))\n x = conv(inp, 2, 16)\n x = pool(x)\n \n x = conv(x, 3, 32)\n x = pool(x)\n \n x = layers.Flatten()(x)\n return dense(x, 1, 0.4, act='linear')","metadata":{"execution":{"iopub.status.busy":"2021-11-04T19:19:35.521508Z","iopub.execute_input":"2021-11-04T19:19:35.521863Z","iopub.status.idle":"2021-11-04T19:19:35.534117Z","shell.execute_reply.started":"2021-11-04T19:19:35.521829Z","shell.execute_reply":"2021-11-04T19:19:35.533158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n num_data = batch_size\n num_inducing = 10\n output_dim = 1\n\n kernel = gpflow.kernels.SquaredExponential()\n inducing_variable = gpflow.inducing_variables.InducingPoints(\n np.ones((num_inducing,output_dim))\n )\n gp_layer = gpflux.layers.GPLayer(\n kernel, inducing_variable, num_data=num_data * dim[2], num_latent_gps=output_dim\n )\n \n likelihood = gpflow.likelihoods.Bernoulli()\n\n # So that Keras can track the likelihood variance, we need to provide the likelihood as part of a \"dummy\" layer:\n likelihood_container = gpflux.layers.TrackableLayer()\n likelihood_container.likelihood = likelihood","metadata":{"execution":{"iopub.status.busy":"2021-11-04T19:36:35.362913Z","iopub.execute_input":"2021-11-04T19:36:35.363248Z","iopub.status.idle":"2021-11-04T19:36:35.395375Z","shell.execute_reply.started":"2021-11-04T19:36:35.363219Z","shell.execute_reply":"2021-11-04T19:36:35.394377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tf.device('/gpu:0'):\n\n def build_model():\n inp = keras.Input(shape=(*dim,3))\n latent = VGG(inp)\n out = gp_layer.predict(latent)\n out = out[0]\n out = tf.reshape(out, shape=(-1,dim[2],out.shape[1]))\n out = layers.maximum(tf.unstack(out, axis=1))\n out = layers.Activation('sigmoid')(out)\n out = likelihood_container(out)\n return keras.Model(inputs=inp, outputs=out)\n \n model = build_model()\n \n earlyStopping = EarlyStopping(patience=2, min_delta=0.001, verbose=1)\n \n checkpoint_path = \"training_4/cp.ckpt\"\n checkpoint_dir = os.path.dirname(checkpoint_path)\n # model.load_weights(checkpoint_path)\n\n # Create a callback that saves the model's weights\n cp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n save_best_only=True,\n save_weights_only=True,\n verbose=1)\n \n loss = gpflux.losses.LikelihoodLoss(likelihood)\n model.compile(\n optimizer='adam', \n loss=loss,\n metrics=[keras.metrics.AUC(), 'accuracy']\n )\n \n history = model.fit(train_dataset, validation_data=val_dataset,\n epochs=30, callbacks=[earlyStopping, cp_callback])","metadata":{"execution":{"iopub.status.busy":"2021-11-04T19:36:37.714636Z","iopub.execute_input":"2021-11-04T19:36:37.71502Z","iopub.status.idle":"2021-11-04T20:28:07.494106Z","shell.execute_reply.started":"2021-11-04T19:36:37.71499Z","shell.execute_reply":"2021-11-04T20:28:07.49322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}