--- a
+++ b/BayesClassifierECG.ipynb
@@ -0,0 +1,448 @@
+{
+  "nbformat": 4,
+  "nbformat_minor": 0,
+  "metadata": {
+    "colab": {
+      "name": "BayesClassifierECG-v1.ipynb",
+      "version": "0.3.2",
+      "provenance": [],
+      "collapsed_sections": [
+        "x2Q8Cy8HZ8dB"
+      ]
+    },
+    "kernelspec": {
+      "name": "python3",
+      "display_name": "Python 3"
+    },
+    "accelerator": "GPU"
+  },
+  "cells": [
+    {
+      "metadata": {
+        "id": "DfByttJIWbgQ",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "# Bayesian Classification for ECG Time-Series\n",
+        "\n",
+        "> Copyright 2019 Dave Fernandes. All Rights Reserved.\n",
+        "> \n",
+        "> Licensed under the Apache License, Version 2.0 (the \"License\");\n",
+        "> you may not use this file except in compliance with the License.\n",
+        "> You may obtain a copy of the License at\n",
+        ">\n",
+        "> http://www.apache.org/licenses/LICENSE-2.0\n",
+        ">  \n",
+        "> Unless required by applicable law or agreed to in writing, software\n",
+        "> distributed under the License is distributed on an \"AS IS\" BASIS,\n",
+        "> WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
+        "> See the License for the specific language governing permissions and\n",
+        "> limitations under the License."
+      ]
+    },
+    {
+      "metadata": {
+        "id": "MTd4aHhYWhdN",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Overview\n",
+        "This notebook classifies time-series for segmented heartbeats from ECG lead II recordings. Either of two models (CNN or RNN) can be selected from training and evaluation.\n",
+        "- Data for this analysis should be prepared using the `PreprocessECG.ipynb` notebook from this project.\n",
+        "- Original data is from: https://www.kaggle.com/shayanfazeli/heartbeat"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "hjmdX-HbWdeI",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "import numpy as np\n",
+        "import tensorflow as tf\n",
+        "import tensorflow.keras.layers as keras\n",
+        "import tensorflow_probability as tfp\n",
+        "from tensorflow_probability import distributions as tfd\n",
+        "import matplotlib.pyplot as plt\n",
+        "import pickle\n",
+        "\n",
+        "from google.colab import drive\n",
+        "drive.mount('/content/drive')\n",
+        "\n",
+        "TRAIN_SET = '/content/drive/My Drive/Colab Notebooks/Data/train_set.pickle'\n",
+        "TEST_SET = '/content/drive/My Drive/Colab Notebooks/Data/test_set.pickle'\n",
+        "\n",
+        "BATCH_SIZE = 125\n",
+        "\n",
+        "with open(TEST_SET, 'rb') as file:\n",
+        "    test_set = pickle.load(file)\n",
+        "    x_test = test_set['x'].astype('float32')\n",
+        "    y_test = test_set['y'].astype('int32')\n",
+        "\n",
+        "with open(TRAIN_SET, 'rb') as file:\n",
+        "    train_set = pickle.load(file)\n",
+        "    x_train = train_set['x'].astype('float32')\n",
+        "    y_train = train_set['y'].astype('int32')\n",
+        "\n",
+        "TRAIN_COUNT = len(y_train)\n",
+        "TEST_COUNT = len(y_test)\n",
+        "BATCHES_PER_EPOCH = TRAIN_COUNT // BATCH_SIZE\n",
+        "TEST_BATCHES_PER_EPOCH = TEST_COUNT // BATCH_SIZE\n",
+        "INPUT_SIZE = np.shape(x_train)[1]\n",
+        "\n",
+        "print('Train count:', TRAIN_COUNT, 'Test count:', TEST_COUNT)"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "x2Q8Cy8HZ8dB",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Input Datasets"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "uiTIHzr5Wsmn",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "def combined_dataset(features, labels):\n",
+        "    assert features.shape[0] == labels.shape[0]\n",
+        "    dataset = tf.data.Dataset.from_tensor_slices((np.expand_dims(features, axis=-1), labels))\n",
+        "    return dataset\n",
+        "\n",
+        "# For training\n",
+        "def train_input_fn():\n",
+        "    dataset = combined_dataset(x_train, y_train)\n",
+        "    return dataset.shuffle(TRAIN_COUNT, reshuffle_each_iteration=True).repeat().batch(BATCH_SIZE, drop_remainder=True).prefetch(1)\n",
+        "\n",
+        "# For evaluation and metrics\n",
+        "def eval_input_fn():\n",
+        "    dataset = combined_dataset(x_test, y_test)\n",
+        "    return dataset.repeat().batch(BATCH_SIZE).prefetch(1)\n",
+        "\n",
+        "training_batches = train_input_fn()\n",
+        "training_iterator = tf.compat.v1.data.make_one_shot_iterator(training_batches)\n",
+        "heldout_iterator = tf.compat.v1.data.make_one_shot_iterator(eval_input_fn())\n",
+        "\n",
+        "# Combine these into a feedable iterator that can switch between training\n",
+        "# and validation inputs.\n",
+        "handle = tf.compat.v1.placeholder(tf.string, shape=[])\n",
+        "feedable_iterator = tf.compat.v1.data.Iterator.from_string_handle(handle, training_batches.output_types, training_batches.output_shapes)\n",
+        "series, labels = feedable_iterator.get_next()"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "PnGsC48GaGTk",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Define the model\n",
+        "#### Bayesian CNN Model\n",
+        "* The convolutional model is taken from: https://arxiv.org/pdf/1805.00794.pdf\n",
+        "\n",
+        "Model consists of:\n",
+        "* An initial 1-D convolutional layer\n",
+        "* 5 repeated residual blocks (`same` padding)\n",
+        "* A fully-connected layer\n",
+        "* A linear layer with softmax output\n",
+        "* Flipout layers are used instead of standard layers"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "Q8XdxTVYaO1q",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "KL_ANNEALING = 30\n",
+        "\n",
+        "MODEL_PATH = '/content/drive/My Drive/Colab Notebooks/Models/BayesianCNN/BNN.tfmodel'\n",
+        "\n",
+        "def kernel_prior(dtype, shape, name, trainable, add_variable_fn):\n",
+        "    return tfp.layers.default_multivariate_normal_fn(dtype, shape, name, trainable, add_variable_fn)\n",
+        "#     return tfd.Horseshoe(scale=5.0)\n",
+        "\n",
+        "#     mix = 0.75\n",
+        "#     mixture = tfd.Mixture(name=name,\n",
+        "#         cat=tfd.Deterministic([mix, 1. - mix]),\n",
+        "#         components=[tfd.Normal(loc=0., scale=1.), tfd.Normal(loc=0., scale=7.)])\n",
+        "#     return mixture\n",
+        "\n",
+        "def conv_unit(unit, input_layer):\n",
+        "    s = '_' + str(unit)\n",
+        "    layer = tfp.layers.Convolution1DFlipout(name='Conv1' + s, filters=32, kernel_size=5, strides=1, padding='same', activation='relu', kernel_prior_fn=kernel_prior)(input_layer)\n",
+        "    layer = tfp.layers.Convolution1DFlipout(name='Conv2' + s, filters=32, kernel_size=5, strides=1, padding='same', activation=None, kernel_prior_fn=kernel_prior)(layer )\n",
+        "    layer = keras.Add(name='ResidualSum' + s)([layer, input_layer])\n",
+        "    layer = keras.Activation(\"relu\", name='Act' + s)(layer)\n",
+        "    layer = keras.MaxPooling1D(name='MaxPool' + s, pool_size=5, strides=2)(layer)\n",
+        "    return layer\n",
+        "\n",
+        "def model_fn(input_shape):\n",
+        "    time_series = tf.keras.layers.Input(shape=input_shape, dtype='float32')\n",
+        "    \n",
+        "    current_layer = tfp.layers.Convolution1DFlipout(filters=32, kernel_size=5, strides=1, kernel_prior_fn=kernel_prior)(time_series)\n",
+        "\n",
+        "    for i in range(5):\n",
+        "        current_layer = conv_unit(i + 1, current_layer)\n",
+        "\n",
+        "    current_layer = keras.Flatten()(current_layer)\n",
+        "    current_layer = tfp.layers.DenseFlipout(32, name='FC1', activation='relu', kernel_prior_fn=kernel_prior)(current_layer)\n",
+        "    logits = tfp.layers.DenseFlipout(5, name='Output', kernel_prior_fn=kernel_prior)(current_layer)\n",
+        "    \n",
+        "    model = tf.keras.Model(inputs=time_series, outputs=logits, name='bayes_cnn_model')\n",
+        "    return model\n",
+        "  \n",
+        "# Compute the negative Evidence Lower Bound (ELBO) loss\n",
+        "t = tf.compat.v1.Variable(0.0)\n",
+        "\n",
+        "def loss_fn(labels, logits):\n",
+        "    labels_distribution = tfd.Categorical(logits=logits)\n",
+        "\n",
+        "    # Perform KL annealing. The optimal number of annealing steps\n",
+        "    # depends on the dataset and architecture.\n",
+        "    kl_regularizer = t / (KL_ANNEALING * BATCHES_PER_EPOCH)\n",
+        "\n",
+        "    # Compute the -ELBO as the loss. The kl term is annealed from 0 to 1 over\n",
+        "    # the epochs specified by the kl_annealing flag.\n",
+        "    log_likelihood = labels_distribution.log_prob(labels)\n",
+        "    neg_log_likelihood = -tf.reduce_mean(input_tensor=log_likelihood)\n",
+        "    kl = sum(model.losses) / len(x_train) * tf.minimum(1.0, kl_regularizer)\n",
+        "    return neg_log_likelihood + kl, kl, kl_regularizer, labels_distribution\n",
+        "\n",
+        "model = model_fn([INPUT_SIZE, 1])\n",
+        "model.summary()"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "1h96cjFraUo_",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Train model"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "9uHtgxoLynkU",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "INITIAL_LEARNING_RATE = 0.0001\n",
+        "EPOCHS = 100\n",
+        "\n",
+        "assert (EPOCHS > 0)\n",
+        "\n",
+        "logits = model(series)\n",
+        "loss, kl, kl_reg, labels_distribution = loss_fn(labels, logits)\n",
+        "\n",
+        "# Build metrics for evaluation. Predictions are formed from a single forward\n",
+        "# pass of the probabilistic layers. They are cheap but noisy\n",
+        "# predictions.\n",
+        "predictions = tf.argmax(input=logits, axis=1)\n",
+        "with tf.compat.v1.name_scope(\"train\"):\n",
+        "    train_accuracy, train_accuracy_update_op = tf.compat.v1.metrics.accuracy(labels=labels, predictions=predictions)\n",
+        "    opt = tf.compat.v1.train.AdamOptimizer(INITIAL_LEARNING_RATE)\n",
+        "    train_op = opt.minimize(loss)\n",
+        "    update_step_op = tf.compat.v1.assign(t, t + 1)\n",
+        "\n",
+        "with tf.compat.v1.name_scope(\"valid\"):\n",
+        "    valid_accuracy, valid_accuracy_update_op = tf.compat.v1.metrics.accuracy(labels=labels, predictions=predictions)\n",
+        "\n",
+        "init_op = tf.group(tf.compat.v1.global_variables_initializer(), tf.compat.v1.local_variables_initializer())\n",
+        "\n",
+        "stream_vars_valid = [\n",
+        "    v for v in tf.compat.v1.local_variables() if \"valid/\" in v.name\n",
+        "]\n",
+        "reset_valid_op = tf.compat.v1.variables_initializer(stream_vars_valid)\n",
+        "\n",
+        "with tf.compat.v1.Session() as sess:\n",
+        "    sess.run(init_op)\n",
+        "\n",
+        "    # Run the training loop\n",
+        "    train_handle = sess.run(training_iterator.string_handle())\n",
+        "    heldout_handle = sess.run(heldout_iterator.string_handle())\n",
+        "    training_steps = EPOCHS * BATCHES_PER_EPOCH\n",
+        "    \n",
+        "    for step in range(training_steps):\n",
+        "        _ = sess.run([train_op, train_accuracy_update_op, update_step_op], feed_dict={handle: train_handle})\n",
+        "\n",
+        "        # Manually print the frequency\n",
+        "        if step % (BATCHES_PER_EPOCH // 5) == 0:\n",
+        "            loss_value, accuracy_value, kl_value, kl_reg_value = sess.run([loss, train_accuracy, kl, kl_reg], feed_dict={handle: train_handle})\n",
+        "            print(\"   Loss: {:.3f} Accuracy: {:.3f} KL: {:.3f} KL-reg: {:.3f}\".format(loss_value, accuracy_value, kl_value, kl_reg_value))\n",
+        "\n",
+        "        if (step + 1) % BATCHES_PER_EPOCH == 0:\n",
+        "            # Calculate validation accuracy\n",
+        "            for _ in range(TEST_BATCHES_PER_EPOCH):\n",
+        "                sess.run(valid_accuracy_update_op, feed_dict={handle: heldout_handle})\n",
+        "            \n",
+        "            valid_value = sess.run(valid_accuracy, feed_dict={handle: heldout_handle})\n",
+        "            print(\"Epoch: {:>3d} Validation Accuracy: {:.3f}\".format((step + 1) // BATCHES_PER_EPOCH, valid_value))\n",
+        "\n",
+        "            sess.run(reset_valid_op)\n",
+        "            \n",
+        "    model.save_weights(MODEL_PATH)"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "Ml9cWK8j0vEU",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Evaluate model"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "p40qKgiIIVNd",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "NUM_MONTE_CARLO = 1000\n",
+        "\n",
+        "model.load_weights(MODEL_PATH)\n",
+        "\n",
+        "mc_counts = np.zeros((TEST_COUNT, 5))\n",
+        "x = np.expand_dims(x_test, -1)\n",
+        "sample_index = np.arange(TEST_COUNT)\n",
+        "\n",
+        "for i in range(NUM_MONTE_CARLO):\n",
+        "    y_pred = np.argmax(model.predict(x), axis=1)\n",
+        "    mc_counts[sample_index, y_pred] += 1\n",
+        "    \n",
+        "y_pred = np.argmax(mc_counts, axis=1)\n",
+        "y_prob = mc_counts[sample_index, y_pred] / NUM_MONTE_CARLO\n",
+        "\n",
+        "y_prob_correct = y_prob[y_pred == y_test]\n",
+        "y_prob_mis = y_prob[y_pred != y_test]"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "dWh_nMZN9J0N",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Check probability estimates"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "K0lRvUS9gQIk",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "from astropy.stats import binom_conf_interval\n",
+        "\n",
+        "_, _, _ = plt.hist(y_prob, 10, (0, 1))\n",
+        "plt.xlabel('Belief')\n",
+        "plt.ylabel('Count')\n",
+        "plt.title('All Predictions')\n",
+        "plt.show();\n",
+        "\n",
+        "n_all, bins = np.histogram(y_prob, 10, (0, 1))\n",
+        "n_correct, bins = np.histogram(y_prob_correct, 10, (0, 1))\n",
+        "\n",
+        "f_correct = n_correct / np.clip(n_all, 1, None)\n",
+        "f_bins = 0.5 * (bins[:-1] + bins[1:])\n",
+        "\n",
+        "n_correct = n_correct[n_all > 0]\n",
+        "n_total = n_all[n_all > 0]\n",
+        "f_correct = n_correct / n_total\n",
+        "f_bins = f_bins[n_all > 0]\n",
+        "\n",
+        "lower_bound, upper_bound = binom_conf_interval(n_correct, n_total)\n",
+        "error_bars = np.array([f_correct - lower_bound, upper_bound - f_correct])\n",
+        "\n",
+        "plt.plot([0., 1.], [0., 1.])\n",
+        "plt.errorbar(f_bins, f_correct, yerr=error_bars, fmt='o')\n",
+        "plt.xlabel('Monte Carlo Probability')\n",
+        "plt.ylabel('Frequency')\n",
+        "plt.title('Correct Predictions')\n",
+        "plt.show();"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "ggvEGluM9TEw",
+        "colab_type": "text"
+      },
+      "cell_type": "markdown",
+      "source": [
+        "### Compute metrics"
+      ]
+    },
+    {
+      "metadata": {
+        "id": "TeV6NtwnP0Lm",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        "import sklearn.metrics as skm\n",
+        "import seaborn\n",
+        "\n",
+        "# Classification report\n",
+        "report = skm.classification_report(y_test, y_pred)\n",
+        "print(report)\n",
+        "\n",
+        "# Confusion matrix\n",
+        "cm = skm.confusion_matrix(y_test, y_pred)\n",
+        "seaborn.heatmap(cm, annot=True,annot_kws={\"size\": 16})"
+      ],
+      "execution_count": 0,
+      "outputs": []
+    },
+    {
+      "metadata": {
+        "id": "MnQeNtCkV3Uv",
+        "colab_type": "code",
+        "colab": {}
+      },
+      "cell_type": "code",
+      "source": [
+        ""
+      ],
+      "execution_count": 0,
+      "outputs": []
+    }
+  ]
+}
\ No newline at end of file