notebooks/Train Patch Classifier Weak Labeling - Ensemble, SVM, 1px, LARGE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f1c69676",
   "metadata": {},
   "source": [
    "# Train a Patch Classifier via Weak Labeling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de8c4dcd",
   "metadata": {},
   "source": [
    "This notebook takes strongly labeled data, generates an ensemble of patch classifiers, trains an SVM, then uses all of this to generate weak labels on a ton of unlabeled data. This unlabeled data with weak labels is then combined with the strongly labeled data to train a new model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d803c9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#system and IO\n",
    "import sys\n",
    "import os\n",
    "import pickle\n",
    "\n",
    "#Math and display tools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "from tqdm import tqdm\n",
    "\n",
    "#Machine learning\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "from tensorflow.keras.preprocessing.image import ImageDataGenerator\n",
    "from sklearn import svm\n",
    "from sklearn.metrics import classification_report\n",
    "\n",
    "ENSEMBLE_SIZE = 32\n",
    "BATCH_SIZE = 128\n",
    "EPOCHS = 256\n",
    "\n",
    "data_root = \"spectrogram_patches_3mo-mosaics_2x-int\"\n",
    "data_paths = list()\n",
    "#find all data files in the directory\n",
    "for f in os.listdir(data_root):\n",
    "    if \"arrays.pkl\" in f:\n",
    "        data_paths.append(f[:-17])\n",
    "\n",
    "#Means and standard deviations used for normalization\n",
    "MEANS = [1367.8407, 1104.4116, 1026.8099, 856.1295, \n",
    "         1072.1476, 1880.3287, 2288.875, 2104.5999, \n",
    "         2508.7764, 305.3795, 1686.0194, 946.1319]\n",
    "STDEVS = [249.14418, 317.69983, 340.8048, 467.8019, \n",
    "          390.11594, 529.972, 699.90826, 680.56006, \n",
    "          798.34937, 108.10846, 651.8683, 568.5347]\n",
    "\n",
    "#class to one-hot (keras has a function for this, but it has issues on my machine)\n",
    "def to_onehot(array, num_classes):\n",
    "    array = list(array) #ensure the array is a list\n",
    "    for i in range(len(array)):\n",
    "        zeros = [0]*num_classes\n",
    "        zeros[int(array[i])] = 1\n",
    "        array[i] = zeros\n",
    "    return np.float32(array)\n",
    "\n",
    "#takes a random crop of an image\n",
    "def random_crop2D(img, out_dim):\n",
    "    in_dim = (img.shape[0], img.shape[1])\n",
    "    top_left = (random.randrange(0, in_dim[0] - out_dim[0]), random.randrange(0, in_dim[1] - out_dim[1]))\n",
    "    return img[top_left[0]:(top_left[0]+out_dim[0]), top_left[1]:(top_left[1]+out_dim[1]),:].copy()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ce213b1",
   "metadata": {},
   "source": [
    "# Load the strongly labeled data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaffbc47",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_data = list()\n",
    "labels = list()\n",
    "\n",
    "#open each data file and read\n",
    "for path in data_paths:\n",
    "    with open(os.path.join(data_root, path + \"_patch_arrays.pkl\"), \"rb\") as f:\n",
    "        patch = pickle.load(f) #load set of patches\n",
    "        for img in patch:\n",
    "            raw_data.append(np.asarray(img)[:,:48,:48,:]) #naive upper left crop\n",
    "    with open(os.path.join(data_root, path + \"_patch_array_labels.pkl\"), \"rb\") as f:\n",
    "        label = pickle.load(f)\n",
    "        labels += label\n",
    "raw_data = np.float32(raw_data) #convert to np array\n",
    "labels = np.float32(to_onehot(labels, 2)) #convert labels to one-hot, then to np array\n",
    "\n",
    "#Normalize training data\n",
    "normalized_data = raw_data.copy()\n",
    "for i in range(0, 12):\n",
    "    #normalize each channel to global unit norm\n",
    "    normalized_data[:,:,:,:,i] = (raw_data[:,:,:,:,i] - MEANS[i])/STDEVS[i]\n",
    "\n",
    "#produce 48x48x24 dataset\n",
    "flattened_data = np.concatenate((normalized_data[:, 0, :, :, :], normalized_data[:, 1, :, :, :]), 3)\n",
    "\n",
    "#produce a (28*28*24,) for the SVM training dataset via randomly cropping the 48x48x24 dataset\n",
    "X_cropped = list()\n",
    "Y_cropped = list()\n",
    "for i in range(0, flattened_data.shape[0]):\n",
    "    for _ in range(0, 4):\n",
    "        X_cropped.append(random_crop2D(flattened_data[i], (28, 28)))\n",
    "        Y_cropped.append(labels[i])\n",
    "\n",
    "#shuffle cropped data\n",
    "zipped = list(zip(X_cropped, Y_cropped)) #inefficient, but it works\n",
    "random.shuffle(zipped)\n",
    "X_cropped, Y_cropped = zip(*zipped)\n",
    "X_cropped = np.float32(X_cropped)\n",
    "Y_cropped = np.float32(Y_cropped)\n",
    "\n",
    "#Create SVM data\n",
    "X_flat = np.asarray(X_cropped)\n",
    "X_flat = X_flat.reshape(X_flat.shape[0],-1) #flatten data to one dimension for SVM\n",
    "Y_flat = np.asarray(Y_cropped)\n",
    "\n",
    "#shuffle 2D-data\n",
    "zipped = list(zip(flattened_data, labels))\n",
    "random.shuffle(zipped)\n",
    "flattened_data, labels = zip(*zipped)\n",
    "flattened_data = np.float32(flattened_data)\n",
    "labels = np.float32(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86c896b8",
   "metadata": {},
   "source": [
    "# Train ensemble of patch classifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd1b0d7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def crop_generator(batch_size):\n",
    "    augmentation_parameters = {\n",
    "        'featurewise_center': False,\n",
    "        'rotation_range': 360,\n",
    "        'vertical_flip': True,\n",
    "        'horizontal_flip': True,\n",
    "        'fill_mode': 'reflect'\n",
    "    }\n",
    "\n",
    "    datagen = ImageDataGenerator(**augmentation_parameters) #create the datagenerator for augmentation\n",
    "    generator = datagen.flow((flattened_data, labels), batch_size=batch_size)\n",
    "    \n",
    "    while True:\n",
    "        augmented_X, augmented_Y = generator.next()\n",
    "        \n",
    "        cropped_X = list()\n",
    "        \n",
    "        for i in range(0, augmented_X.shape[0]):\n",
    "            cropped_X.append(random_crop2D(augmented_X[i], (28, 28))) #random crop\n",
    "\n",
    "        yield (np.asarray(cropped_X), augmented_Y)\n",
    "\n",
    "for i in range(0, ENSEMBLE_SIZE):\n",
    "    crop_gen = crop_generator(BATCH_SIZE)\n",
    "\n",
    "    #define model as determined from the hyperparameter optimization experiment\n",
    "    input_img = keras.Input(shape=(28, 28, 24))\n",
    "    x = layers.Conv2D(32, (5, 5), padding='same')(input_img)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.Conv2D(32, (5, 5), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.Conv2D(32, (5, 5), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.MaxPooling2D(2)(x)\n",
    "\n",
    "    x = layers.Conv2D(32, (4, 4), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.Conv2D(32, (4, 4), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.Conv2D(32, (4, 4), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.MaxPooling2D(2)(x)\n",
    "\n",
    "    x = layers.Conv2D(32, (3, 3), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.Conv2D(32, (3, 3), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.Conv2D(32, (3, 3), padding='same')(x)\n",
    "    x = layers.LeakyReLU()(x)\n",
    "    x = layers.MaxPooling2D(2)(x)\n",
    "\n",
    "    x = layers.Flatten()(x)\n",
    "    x = layers.Dropout(0.5)(x)\n",
    "    x = layers.Dense(512, activation=layers.ELU())(x)\n",
    "    x = layers.Dropout(0.5)(x)\n",
    "    x = layers.Dense(512, activation=layers.ELU())(x)\n",
    "    x = layers.Dropout(0.5)(x)\n",
    "    out = layers.Dense(2, activation=\"softmax\")(x)\n",
    "\n",
    "    model = keras.Model(input_img, out, name=\"classifier\")\n",
    "\n",
    "    model.compile(loss=\"binary_crossentropy\", #binary_crossentropy is default for two classes\n",
    "                  optimizer=keras.optimizers.Adam(epsilon=0.001), #haven't tried other optimizers, might be a good idea, but Adam is quite good usually\n",
    "                  metrics=[\"accuracy\"])\n",
    "    \n",
    "    mcp_save = keras.callbacks.ModelCheckpoint(\"patch_ensemble/model_\" + str(i) + \".hdf5\", \n",
    "                                               save_best_only=True, monitor='accuracy', verbose=1, mode='max')\n",
    "    reduce_lr_loss = keras.callbacks.ReduceLROnPlateau(monitor='loss', factor=0.1, patience=7, verbose=1, min_delta=1e-4, mode='min')\n",
    "\n",
    "\n",
    "    #train the model\n",
    "    model.fit(crop_gen,\n",
    "              epochs=EPOCHS,\n",
    "              steps_per_epoch=flattened_data.shape[0]//BATCH_SIZE,\n",
    "              callbacks=[mcp_save, reduce_lr_loss],\n",
    "              verbose = 1\n",
    "             )\n",
    "    #model.save(\"patch_ensemble/model_\" + str(i) + \".hdf5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cfb1d0f",
   "metadata": {},
   "source": [
    "# Train SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5b453e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Train SVM\n",
    "svm_classifier = svm.SVC()\n",
    "svm_classifier.fit(X_flat, Y_flat[:,1]) #this will take a loooong time!\n",
    "\n",
    "pickle.dump(svm_classifier, open(\"SVM.pkl\", \"wb\")) #save SVM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fa32df3",
   "metadata": {},
   "source": [
    "# Generate Weak Labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f4b835c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate weak labels\n",
    "def predict_ensemble(X, models, return_disagreement=False):\n",
    "    preds = list() #list of model predictions, will be (len(models), len(X), 2) in dimension\n",
    "    \n",
    "    #run ensemble\n",
    "    for model in models:\n",
    "        preds.append(model.predict(X))\n",
    "    \n",
    "    #convert soft predictions to absolute\n",
    "    for pred in preds:\n",
    "        for i in range(0, pred.shape[0]):\n",
    "            y = [0, 0]\n",
    "            y[np.argmax(pred[i])] = 1.0\n",
    "            pred[i] = y\n",
    "    \n",
    "    #aggregate predictions in votes\n",
    "    pred_sum = preds[0]\n",
    "    for i in range(1, len(preds)):\n",
    "        pred_sum += preds[i]\n",
    "    \n",
    "    #pick the prediction with the highest votes\n",
    "    for i in range(0, len(pred_sum)):\n",
    "        y = [0, 0]\n",
    "        y[np.argmax(pred_sum[i])] = 1.0\n",
    "        pred_sum[i] = y\n",
    "    \n",
    "    #compute the std of the votes (a measure of disagreement)\n",
    "    if return_disagreement:\n",
    "        return np.asarray(pred_sum), np.std(np.asarray(preds)[:,:,0],0)\n",
    "    \n",
    "    return np.asarray(pred_sum)\n",
    "\n",
    "def create_weak_labels(X, patch_ensemble, pixel_model, SVM, weights, return_all=False):\n",
    "    print(\"Running ensemble\")\n",
    "    patch_preds, patch_disagreement = predict_ensemble(X, patch_ensemble, True) #patch predictions\n",
    "    SVM_preds = list()\n",
    "    print(\"Running pixel classifier\")\n",
    "    \n",
    "    #we need to do some pretty massive adjustments to the data to work with the px classifier\n",
    "    px_2D = np.reshape(X, (-1,24))\n",
    "    px_2D = np.asarray([px_2D[:,:12], px_2D[:,12:]]) #this essentially doubles the RAM requirement, not optimal\n",
    "    px_2D = np.transpose(px_2D, (1,2,0))\n",
    "    px_2D.shape = (px_2D.shape[0], px_2D.shape[1], px_2D.shape[2], 1)\n",
    "    for i in range(0, 12):\n",
    "        px_2D[:,i,:,:] = (px_2D[:,i,:,:]*STDEVS[i]+MEANS[i])/3000.0 #undo the original normalization\n",
    "        \n",
    "    #Run the 1 px classifier\n",
    "    pixel_out = pixel_model.predict(px_2D)\n",
    "    pixel_out.shape = (-1, 28*28, 2)\n",
    "    per_image_out = np.mean(pixel_out, axis=1)\n",
    "    per_image_out = np.int32(per_image_out[:,1] > 0.02) #threshold at 0.02, a predetermined value\n",
    "    per_image_out = to_onehot(per_image_out, 2)\n",
    "    \n",
    "    print(\"Running SVM\")\n",
    "    for i in tqdm(range(0, X.shape[0])):\n",
    "        SVM_preds.append(to_onehot(list(SVM.predict(X[i].reshape((1,-1)))), 2)[0])\n",
    "    SVM_preds = np.asarray(SVM_preds)\n",
    "    \n",
    "    weak_labels = list()\n",
    "    for i in range(0, X.shape[0]):\n",
    "        weak_label = patch_preds[i]*(1-2*patch_disagreement[i])\n",
    "        \n",
    "        #ensure weak label total adds to 1.0\n",
    "        if weak_label[0] == 0:\n",
    "            weak_label[0] = 1 - weak_label[1]\n",
    "        else:\n",
    "            weak_label[1] = 1 - weak_label[0]\n",
    "            \n",
    "        weak_label = weak_label*weights[0] + per_image_out[i]*weights[1] + SVM_preds[i]*weights[2]\n",
    "        \n",
    "        total_prob = np.sum(weak_label)\n",
    "        if total_prob != 1:\n",
    "            weak_label /= total_prob\n",
    "        weak_labels.append(weak_label)\n",
    "    if return_all:\n",
    "        return np.asarray(weak_labels), patch_preds, patch_disagreement, per_image_out, SVM_preds\n",
    "    return np.asarray(weak_labels)\n",
    "\n",
    "def create_weak_labels_bayesian(X, patch_ensemble, pixel_model, SVM, return_all=False):\n",
    "    print(\"Running ensemble\")\n",
    "    patch_preds, patch_disagreement = predict_ensemble(X, patch_ensemble, True) #patch predictions\n",
    "    SVM_preds = list()\n",
    "    print(\"Running pixel classifier\")\n",
    "    \n",
    "    #we need to do some pretty massive adjustments to the data to work with the px classifier\n",
    "    px_2D = np.reshape(X, (-1,24))\n",
    "    px_2D = np.asarray([px_2D[:,:12], px_2D[:,12:]]) #this essentially doubles the RAM requirement, not optimal\n",
    "    px_2D = np.transpose(px_2D, (1,2,0))\n",
    "    px_2D.shape = (px_2D.shape[0], px_2D.shape[1], px_2D.shape[2], 1)\n",
    "    for i in range(0, 12):\n",
    "        px_2D[:,i,:,:] = (px_2D[:,i,:,:]*STDEVS[i]+MEANS[i])/3000.0 #undo the original normalization\n",
    "        \n",
    "    #Run the 1 px classifier\n",
    "    pixel_out = pixel_model.predict(px_2D)\n",
    "    pixel_out.shape = (-1, 28*28, 2)\n",
    "    per_image_out = np.mean(pixel_out, axis=1)\n",
    "    per_image_out = np.int32(per_image_out[:,1] > 0.02) #threshold at 0.02, a predetermined value\n",
    "    \n",
    "    print(\"Running SVM\")\n",
    "    for i in tqdm(range(0, X.shape[0])):\n",
    "        SVM_preds.append(SVM.predict(X[i].reshape((1,-1))))\n",
    "    SVM_preds = np.asarray(SVM_preds)\n",
    "    \n",
    "    weak_labels = list()\n",
    "    for i in range(0, X.shape[0]):\n",
    "        weak_label = patch_preds[i]*(1-2*patch_disagreement[i])\n",
    "        \n",
    "        #ensure weak label total adds to 1.0\n",
    "        if weak_label[0] == 0:\n",
    "            weak_label[0] = 1 - weak_label[1]\n",
    "        else:\n",
    "            weak_label[1] = 1 - weak_label[0]\n",
    "        \n",
    "        '''\n",
    "            ===============1 PX CLASSIFIER STATS===============\n",
    "            Sensitivity of 1px (recall of positive class): 0.83\n",
    "            Specificity of 1px (recall of negative class): 0.89\n",
    "            True positive rate: 0.83\n",
    "            False negative rate: 0.17\n",
    "            True negative rate: 0.89\n",
    "            False positive rate: 0.11\n",
    "        '''\n",
    "        #update hypothesis for 1px classifier\n",
    "        if per_image_out[i] == 0: #No plastic detected\n",
    "            weak_label[0] = weak_label[0]*0.89/(0.89*weak_label[0] + (1-weak_label[0])*0.17)\n",
    "            weak_label[1] = 1 - weak_label[0]\n",
    "        elif per_image_out[i] == 1: #plastic detected\n",
    "            weak_label[1] = weak_label[1]*0.83/(weak_label[1]*0.83 + (1-weak_label[1])*0.11)\n",
    "            weak_label[0] = 1 - weak_label[1]\n",
    "        else:\n",
    "            print(\"Error in 1px Bayesian update, this should not be reached!\")\n",
    "        \n",
    "        '''\n",
    "            ===============SVM CLASSIFIER STATS===============\n",
    "            Sensitivity of 1px (recall of positive class): 0.78\n",
    "            Specificity of 1px (recall of negative class): 0.90\n",
    "            True positive rate: 0.78\n",
    "            False negative rate: 0.22\n",
    "            True negative rate: 0.90\n",
    "            False positive rate: 0.10\n",
    "        '''\n",
    "        if SVM_preds[i] == 0: #No plastic detected\n",
    "            weak_label[0] = weak_label[0]*0.9/(weak_label[0]*0.9 + (1-weak_label[0])*0.22)\n",
    "            weak_label[1] = 1 - weak_label[0]\n",
    "        elif SVM_preds[i] == 1: #plastic detected\n",
    "            weak_label[1] = weak_label[1]*0.78/(weak_label[1]*0.78 + (1-weak_label[1])*0.10)\n",
    "            weak_label[0] = 1 - weak_label[1]\n",
    "        else:\n",
    "            print(\"Error in 1px Bayesian update, this should not be reached!\")\n",
    "        \n",
    "        weak_labels.append(weak_label)\n",
    "    if return_all:\n",
    "        return np.asarray(weak_labels), patch_preds, patch_disagreement, per_image_out, SVM_preds\n",
    "    return np.asarray(weak_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99082981",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load unlabeled data\n",
    "unlabeled_norm = np.load(\"unlabeled_28x28x24_99.npy\")\n",
    "\n",
    "#Load ensemble models\n",
    "models = list()\n",
    "for i in range(0, ENSEMBLE_SIZE):\n",
    "    models.append(keras.models.load_model(\"patch_ensemble_v1.0/model_\" + str(i) + \".hdf5\"))\n",
    "\n",
    "#Load 1px model\n",
    "pixel_model = keras.models.load_model(\"spectrogram_v0.0.11_2021-07-13.h5\")\n",
    "\n",
    "#Load SVM\n",
    "SVM = pickle.load( open( \"SVM.pkl\", \"rb\" ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f261d8c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create weak labels\n",
    "weak_labels = create_weak_labels_bayesian(unlabeled_norm, models, pixel_model, SVM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "051afd31",
   "metadata": {},
   "outputs": [],
   "source": [
    "#save new data\n",
    "with open('weak_data.npy', 'wb') as f:\n",
    "    np.save(f, np.concatenate((X_cropped,unlabeled_norm)))\n",
    "with open('weak_labels.npy', 'wb') as f:\n",
    "    np.save(f, np.concatenate((Y_cropped,weak_labels)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b31ab100",
   "metadata": {},
   "source": [
    "# Train new model on weakly labeled data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc9019c1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1487eb4b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#train new model on weak data\n",
    "X_weak = np.load(\"../Weak Data/weak_data_v2.0.npy\")\n",
    "Y_weak = np.load(\"../Weak Data/weak_labels_v2.0.npy\")\n",
    "\n",
    "#shuffle data\n",
    "zipped = list(zip(X_weak, Y_weak)) #inefficient, but it works\n",
    "random.shuffle(zipped)\n",
    "X_weak, Y_weak = zip(*zipped)\n",
    "X_weak = np.float32(X_weak)\n",
    "Y_weak = np.float32(Y_weak)\n",
    "\n",
    "#define data augmentation parameters\n",
    "augmentation_parameters = {\n",
    "    'featurewise_center': False,\n",
    "    'rotation_range': 360,\n",
    "    'vertical_flip': True,\n",
    "    'horizontal_flip': True,\n",
    "    'fill_mode': 'reflect'\n",
    "}\n",
    "\n",
    "datagen = ImageDataGenerator(**augmentation_parameters) #create the datagenerator for augmentation\n",
    "#define model as determined from the hyperparameter optimization experiment\n",
    "\n",
    "input_img = keras.Input(shape=(28, 28, 24))\n",
    "x = layers.Conv2D(32, (5, 5), padding='same')(input_img)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.Conv2D(32, (5, 5), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.Conv2D(32, (5, 5), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.MaxPooling2D(2)(x)\n",
    "\n",
    "x = layers.Conv2D(32, (4, 4), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.Conv2D(32, (4, 4), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.Conv2D(32, (4, 4), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.MaxPooling2D(2)(x)\n",
    "\n",
    "x = layers.Conv2D(32, (3, 3), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.Conv2D(32, (3, 3), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.Conv2D(32, (3, 3), padding='same')(x)\n",
    "x = layers.LeakyReLU()(x)\n",
    "x = layers.MaxPooling2D(2)(x)\n",
    "\n",
    "x = layers.Flatten()(x)\n",
    "x = layers.Dropout(0.5)(x)\n",
    "x = layers.Dense(512, activation=layers.ELU())(x)\n",
    "x = layers.Dropout(0.5)(x)\n",
    "x = layers.Dense(512, activation=layers.ELU())(x)\n",
    "x = layers.Dropout(0.5)(x)\n",
    "out = layers.Dense(2, activation=\"softmax\")(x)\n",
    "\n",
    "model = keras.Model(input_img, out, name=\"classifier\")\n",
    "\n",
    "model.compile(loss=\"binary_crossentropy\", #binary_crossentropy is default for two classes\n",
    "              optimizer=keras.optimizers.Adam(epsilon=0.001), #haven't tried other optimizers, might be a good idea, but Adam is quite good usually\n",
    "              metrics=[\"accuracy\"])\n",
    "\n",
    "mcp_save = keras.callbacks.ModelCheckpoint('deep_plastic_model.hdf5', save_best_only=True, \n",
    "                                           monitor='loss', verbose=1, mode='min')\n",
    "reduce_lr_loss = keras.callbacks.ReduceLROnPlateau(monitor='loss', factor=0.1, patience=7, verbose=1, min_delta=1e-4, mode='min')\n",
    "\n",
    "#train the model\n",
    "model.fit(datagen.flow(X_weak, Y_weak, batch_size=BATCH_SIZE), \n",
    "          epochs=EPOCHS,\n",
    "          callbacks=[mcp_save, reduce_lr_loss],\n",
    "          verbose = 1\n",
    "         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6f20231",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "#train new model on weak data\n",
    "X_weak = np.load(\"weak_data.npy\")\n",
    "Y_weak = np.load(\"weak_labels.npy\")\n",
    "\n",
    "#shuffle data\n",
    "zipped = list(zip(X_weak, Y_weak)) #inefficient, but it works\n",
    "random.shuffle(zipped)\n",
    "X_weak, Y_weak = zip(*zipped)\n",
    "X_weak = np.float32(X_weak)\n",
    "Y_weak = np.float32(Y_weak)\n",
    "\n",
    "#define data augmentation parameters\n",
    "augmentation_parameters = {\n",
    "    'featurewise_center': False,\n",
    "    'rotation_range': 360,\n",
    "    'vertical_flip': True,\n",
    "    'horizontal_flip': True,\n",
    "    'fill_mode': 'reflect'\n",
    "}\n",
    "\n",
    "datagen = ImageDataGenerator(**augmentation_parameters) #create the datagenerator for augmentation\n",
    "#define model as determined from the hyperparameter optimization experiment\n",
    "model = keras.Sequential([\n",
    "            keras.Input(shape=(28,28,24)),\n",
    "            layers.Conv2D(32, kernel_size=(3), padding='same'),\n",
    "            layers.LeakyReLU(),\n",
    "            layers.Conv2D(32, kernel_size=(3), padding='same'),\n",
    "            layers.LeakyReLU(),\n",
    "            layers.MaxPooling2D(2),\n",
    "        \n",
    "            layers.Conv2D(32, kernel_size=(4), padding='same'),\n",
    "            layers.LeakyReLU(),\n",
    "            layers.Conv2D(32, kernel_size=(4), padding='same'),\n",
    "            layers.LeakyReLU(),\n",
    "            layers.MaxPooling2D(2),\n",
    "        \n",
    "            layers.Conv2D(32, kernel_size=(2), padding='same'),\n",
    "            layers.LeakyReLU(),\n",
    "            layers.Conv2D(32, kernel_size=(2), padding='same'),\n",
    "            layers.LeakyReLU(),\n",
    "            layers.MaxPooling2D(2),\n",
    "        \n",
    "            layers.Flatten(),\n",
    "            layers.Dropout(0.5),\n",
    "            layers.Dense(512, activation=layers.ELU()),\n",
    "            layers.Dropout(0.5),\n",
    "            layers.Dense(512, activation=layers.ELU()),\n",
    "            layers.Dropout(0.5),\n",
    "            layers.Dense(2, activation=\"softmax\")])\n",
    "model.summary()\n",
    "\n",
    "model.compile(loss=\"binary_crossentropy\", #binary_crossentropy is default for two classes\n",
    "              optimizer=\"adam\", #haven't tried other optimizers, might be a good idea, but Adam is quite good usually\n",
    "              metrics=[\"accuracy\"])\n",
    "\n",
    "mcp_save = keras.callbacks.ModelCheckpoint('model.hdf5', save_best_only=True, \n",
    "                                           monitor='loss', verbose=1, mode='min')\n",
    "reduce_lr_loss = keras.callbacks.ReduceLROnPlateau(monitor='loss', factor=0.1, patience=7, verbose=1, min_delta=1e-4, mode='min')\n",
    "\n",
    "#train the model\n",
    "model.fit(datagen.flow(X_weak, Y_weak, batch_size=BATCH_SIZE), \n",
    "          epochs=EPOCHS,\n",
    "          callbacks=[mcp_save, reduce_lr_loss],\n",
    "          verbose = 1\n",
    "         )\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a868ed5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test ensemble performance\n",
    "def predict_ensemble(X, models, return_disagreement=False):\n",
    "    preds = list() #list of model predictions, will be (len(models), len(X), 2) in dimension\n",
    "    \n",
    "    #run ensemble\n",
    "    for model in models:\n",
    "        preds.append(model.predict(X))\n",
    "    \n",
    "    #convert soft predictions to absolute\n",
    "    for pred in preds:\n",
    "        for i in range(0, pred.shape[0]):\n",
    "            y = [0, 0]\n",
    "            y[np.argmax(pred[i])] = 1.0\n",
    "            pred[i] = y\n",
    "    \n",
    "    #aggregate predictions in votes\n",
    "    pred_sum = preds[0]\n",
    "    for i in range(1, len(preds)):\n",
    "        pred_sum += preds[i]\n",
    "    \n",
    "    #pick the prediction with the highest votes\n",
    "    for i in range(0, len(pred_sum)):\n",
    "        y = [0, 0]\n",
    "        y[np.argmax(pred_sum[i])] = 1.0\n",
    "        pred_sum[i] = y\n",
    "    \n",
    "    #compute the std of the votes (a measure of disagreement)\n",
    "    if return_disagreement:\n",
    "        return np.asarray(pred_sum), np.std(np.asarray(preds)[:,:,0],0)\n",
    "    \n",
    "    return np.asarray(pred_sum)\n",
    "\n",
    "#Load ensemble models\n",
    "models = list()\n",
    "for i in range(0, ENSEMBLE_SIZE):\n",
    "    models.append(keras.models.load_model(\"patch_ensemble_large/model_\" + str(i) + \".hdf5\"))\n",
    "\n",
    "print(classification_report(Y_cropped[:,1], predict_ensemble(X_cropped, models)[:,1] > 0.5, target_names=['No TPA', 'TPA']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b4ee037",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test weak labeled model performance\n",
    "model = keras.models.load_model(\"deep_plastic_model.hdf5\")\n",
    "print(classification_report(Y_cropped[:,1], model.predict(X_cropped)[:,1] > 0.5, target_names=['No TPA', 'TPA']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8857da80",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test SVM Performance\n",
    "SVM = pickle.load( open( \"SVM.pkl\", \"rb\" ) )\n",
    "SVM_preds = list()\n",
    "for i in tqdm(range(0, X_flat.shape[0])):\n",
    "    SVM_preds.append(SVM.predict([X_flat[i]])[0])\n",
    "print(classification_report(Y_flat[:,1], SVM_preds, target_names=['No TPA', 'TPA']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e2e008d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Load ensemble models\n",
    "models = list()\n",
    "for i in range(0, ENSEMBLE_SIZE):\n",
    "    models.append(keras.models.load_model(\"patch_ensemble/model_\" + str(i) + \".hdf5\"))\n",
    "\n",
    "for i in range(0, len(models)):\n",
    "    print(\"=======================MODEL \" + str(i) + \"===========================\")\n",
    "    print(classification_report(Y_cropped[:,1], models[i].predict(X_cropped)[:,1] > 0.5, target_names=['No TPA', 'TPA']))\n",
    "    print(\"=========================================================\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60b7c1ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test weak labeled model performance\n",
    "print(classification_report(Y_cropped[:,1], model.predict(X_cropped)[:,1] > 0.5, target_names=['No TPA', 'TPA']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6c319fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load ensemble models\n",
    "models = list()\n",
    "for i in range(0, ENSEMBLE_SIZE):\n",
    "    models.append(keras.models.load_model(\"patch_ensemble_large/model_\" + str(i) + \".hdf5\"))\n",
    "\n",
    "#Load 1px model\n",
    "pixel_model = keras.models.load_model(\"spectrogram_v0.0.11_2021-07-13.h5\")\n",
    "\n",
    "#Load SVM\n",
    "SVM = pickle.load( open( \"SVM.pkl\", \"rb\" ) )\n",
    "asdf = create_weak_labels_bayesian(X_cropped, models, pixel_model, SVM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d42f23c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(model.history.history['loss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "187255ab",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14ac9e1a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}