<jupyter_start><jupyter_text>Bayesian Classification for ECG Time-Series> Copyright 2019 <NAME>. All Rights Reserved.> > Licensed under the Apache License, Version 2.0 (the "License");> you may not use this file except in compliance with the License.> You may obtain a copy of the License at>> http://www.apache.org/licenses/LICENSE-2.0>  > Unless required by applicable law or agreed to in writing, software> distributed under the License is distributed on an "AS IS" BASIS,> WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.> See the License for the specific language governing permissions and> limitations under the License.  OverviewThis 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.- Data for this analysis should be prepared using the `PreprocessECG.ipynb` notebook from this project.- Original data is from: https://www.kaggle.com/shayanfazeli/heartbeat<jupyter_code>import numpy as np
import tensorflow as tf
import tensorflow.keras.layers as keras
import tensorflow_probability as tfp
from tensorflow_probability import distributions as tfd
import matplotlib.pyplot as plt
import pickle

from google.colab import drive
drive.mount('/content/drive')

TRAIN_SET = '/content/drive/My Drive/Colab Notebooks/Data/train_set.pickle'
TEST_SET = '/content/drive/My Drive/Colab Notebooks/Data/test_set.pickle'

BATCH_SIZE = 125

with open(TEST_SET, 'rb') as file:
    test_set = pickle.load(file)
    x_test = test_set['x'].astype('float32')
    y_test = test_set['y'].astype('int32')

with open(TRAIN_SET, 'rb') as file:
    train_set = pickle.load(file)
    x_train = train_set['x'].astype('float32')
    y_train = train_set['y'].astype('int32')

TRAIN_COUNT = len(y_train)
TEST_COUNT = len(y_test)
BATCHES_PER_EPOCH = TRAIN_COUNT // BATCH_SIZE
TEST_BATCHES_PER_EPOCH = TEST_COUNT // BATCH_SIZE
INPUT_SIZE = np.shape(x_train)[1]

print('Train count:', TRAIN_COUNT, 'Test count:', TEST_COUNT)<jupyter_output><empty_output><jupyter_text>Input Datasets<jupyter_code>def combined_dataset(features, labels):
    assert features.shape[0] == labels.shape[0]
    dataset = tf.data.Dataset.from_tensor_slices((np.expand_dims(features, axis=-1), labels))
    return dataset

# For training
def train_input_fn():
    dataset = combined_dataset(x_train, y_train)
    return dataset.shuffle(TRAIN_COUNT, reshuffle_each_iteration=True).repeat().batch(BATCH_SIZE, drop_remainder=True).prefetch(1)

# For evaluation and metrics
def eval_input_fn():
    dataset = combined_dataset(x_test, y_test)
    return dataset.repeat().batch(BATCH_SIZE).prefetch(1)

training_batches = train_input_fn()
training_iterator = tf.compat.v1.data.make_one_shot_iterator(training_batches)
heldout_iterator = tf.compat.v1.data.make_one_shot_iterator(eval_input_fn())

# Combine these into a feedable iterator that can switch between training
# and validation inputs.
handle = tf.compat.v1.placeholder(tf.string, shape=[])
feedable_iterator = tf.compat.v1.data.Iterator.from_string_handle(handle, training_batches.output_types, training_batches.output_shapes)
series, labels = feedable_iterator.get_next()<jupyter_output><empty_output><jupyter_text>Define the model Bayesian CNN Model* The convolutional model is taken from: https://arxiv.org/pdf/1805.00794.pdfModel consists of:* An initial 1-D convolutional layer* 5 repeated residual blocks (`same` padding)* A fully-connected layer* A linear layer with softmax output* Flipout layers are used instead of standard layers<jupyter_code>KL_ANNEALING = 30

MODEL_PATH = '/content/drive/My Drive/Colab Notebooks/Models/BayesianCNN/BNN.tfmodel'

def kernel_prior(dtype, shape, name, trainable, add_variable_fn):
    return tfp.layers.default_multivariate_normal_fn(dtype, shape, name, trainable, add_variable_fn)
#     return tfd.Horseshoe(scale=5.0)

#     mix = 0.75
#     mixture = tfd.Mixture(name=name,
#         cat=tfd.Deterministic([mix, 1. - mix]),
#         components=[tfd.Normal(loc=0., scale=1.), tfd.Normal(loc=0., scale=7.)])
#     return mixture

def conv_unit(unit, input_layer):
    s = '_' + str(unit)
    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)
    layer = tfp.layers.Convolution1DFlipout(name='Conv2' + s, filters=32, kernel_size=5, strides=1, padding='same', activation=None, kernel_prior_fn=kernel_prior)(layer )
    layer = keras.Add(name='ResidualSum' + s)([layer, input_layer])
    layer = keras.Activation("relu", name='Act' + s)(layer)
    layer = keras.MaxPooling1D(name='MaxPool' + s, pool_size=5, strides=2)(layer)
    return layer

def model_fn(input_shape):
    time_series = tf.keras.layers.Input(shape=input_shape, dtype='float32')
    
    current_layer = tfp.layers.Convolution1DFlipout(filters=32, kernel_size=5, strides=1, kernel_prior_fn=kernel_prior)(time_series)

    for i in range(5):
        current_layer = conv_unit(i + 1, current_layer)

    current_layer = keras.Flatten()(current_layer)
    current_layer = tfp.layers.DenseFlipout(32, name='FC1', activation='relu', kernel_prior_fn=kernel_prior)(current_layer)
    logits = tfp.layers.DenseFlipout(5, name='Output', kernel_prior_fn=kernel_prior)(current_layer)
    
    model = tf.keras.Model(inputs=time_series, outputs=logits, name='bayes_cnn_model')
    return model
  
# Compute the negative Evidence Lower Bound (ELBO) loss
t = tf.compat.v1.Variable(0.0)

def loss_fn(labels, logits):
    labels_distribution = tfd.Categorical(logits=logits)

    # Perform KL annealing. The optimal number of annealing steps
    # depends on the dataset and architecture.
    kl_regularizer = t / (KL_ANNEALING * BATCHES_PER_EPOCH)

    # Compute the -ELBO as the loss. The kl term is annealed from 0 to 1 over
    # the epochs specified by the kl_annealing flag.
    log_likelihood = labels_distribution.log_prob(labels)
    neg_log_likelihood = -tf.reduce_mean(input_tensor=log_likelihood)
    kl = sum(model.losses) / len(x_train) * tf.minimum(1.0, kl_regularizer)
    return neg_log_likelihood + kl, kl, kl_regularizer, labels_distribution

model = model_fn([INPUT_SIZE, 1])
model.summary()<jupyter_output><empty_output><jupyter_text>Train model<jupyter_code>INITIAL_LEARNING_RATE = 0.0001
EPOCHS = 100

assert (EPOCHS > 0)

logits = model(series)
loss, kl, kl_reg, labels_distribution = loss_fn(labels, logits)

# Build metrics for evaluation. Predictions are formed from a single forward
# pass of the probabilistic layers. They are cheap but noisy
# predictions.
predictions = tf.argmax(input=logits, axis=1)
with tf.compat.v1.name_scope("train"):
    train_accuracy, train_accuracy_update_op = tf.compat.v1.metrics.accuracy(labels=labels, predictions=predictions)
    opt = tf.compat.v1.train.AdamOptimizer(INITIAL_LEARNING_RATE)
    train_op = opt.minimize(loss)
    update_step_op = tf.compat.v1.assign(t, t + 1)

with tf.compat.v1.name_scope("valid"):
    valid_accuracy, valid_accuracy_update_op = tf.compat.v1.metrics.accuracy(labels=labels, predictions=predictions)

init_op = tf.group(tf.compat.v1.global_variables_initializer(), tf.compat.v1.local_variables_initializer())

stream_vars_valid = [
    v for v in tf.compat.v1.local_variables() if "valid/" in v.name
]
reset_valid_op = tf.compat.v1.variables_initializer(stream_vars_valid)

with tf.compat.v1.Session() as sess:
    sess.run(init_op)

    # Run the training loop
    train_handle = sess.run(training_iterator.string_handle())
    heldout_handle = sess.run(heldout_iterator.string_handle())
    training_steps = EPOCHS * BATCHES_PER_EPOCH
    
    for step in range(training_steps):
        _ = sess.run([train_op, train_accuracy_update_op, update_step_op], feed_dict={handle: train_handle})

        # Manually print the frequency
        if step % (BATCHES_PER_EPOCH // 5) == 0:
            loss_value, accuracy_value, kl_value, kl_reg_value = sess.run([loss, train_accuracy, kl, kl_reg], feed_dict={handle: train_handle})
            print("   Loss: {:.3f} Accuracy: {:.3f} KL: {:.3f} KL-reg: {:.3f}".format(loss_value, accuracy_value, kl_value, kl_reg_value))

        if (step + 1) % BATCHES_PER_EPOCH == 0:
            # Calculate validation accuracy
            for _ in range(TEST_BATCHES_PER_EPOCH):
                sess.run(valid_accuracy_update_op, feed_dict={handle: heldout_handle})
            
            valid_value = sess.run(valid_accuracy, feed_dict={handle: heldout_handle})
            print("Epoch: {:>3d} Validation Accuracy: {:.3f}".format((step + 1) // BATCHES_PER_EPOCH, valid_value))

            sess.run(reset_valid_op)
            
    model.save_weights(MODEL_PATH)<jupyter_output><empty_output><jupyter_text>Evaluate model<jupyter_code>NUM_MONTE_CARLO = 1000

model.load_weights(MODEL_PATH)

mc_counts = np.zeros((TEST_COUNT, 5))
x = np.expand_dims(x_test, -1)
sample_index = np.arange(TEST_COUNT)

for i in range(NUM_MONTE_CARLO):
    y_pred = np.argmax(model.predict(x), axis=1)
    mc_counts[sample_index, y_pred] += 1
    
y_pred = np.argmax(mc_counts, axis=1)
y_prob = mc_counts[sample_index, y_pred] / NUM_MONTE_CARLO

y_prob_correct = y_prob[y_pred == y_test]
y_prob_mis = y_prob[y_pred != y_test]<jupyter_output><empty_output><jupyter_text>Check probability estimates<jupyter_code>from astropy.stats import binom_conf_interval

_, _, _ = plt.hist(y_prob, 10, (0, 1))
plt.xlabel('Belief')
plt.ylabel('Count')
plt.title('All Predictions')
plt.show();

n_all, bins = np.histogram(y_prob, 10, (0, 1))
n_correct, bins = np.histogram(y_prob_correct, 10, (0, 1))

f_correct = n_correct / np.clip(n_all, 1, None)
f_bins = 0.5 * (bins[:-1] + bins[1:])

n_correct = n_correct[n_all > 0]
n_total = n_all[n_all > 0]
f_correct = n_correct / n_total
f_bins = f_bins[n_all > 0]

lower_bound, upper_bound = binom_conf_interval(n_correct, n_total)
error_bars = np.array([f_correct - lower_bound, upper_bound - f_correct])

plt.plot([0., 1.], [0., 1.])
plt.errorbar(f_bins, f_correct, yerr=error_bars, fmt='o')
plt.xlabel('Monte Carlo Probability')
plt.ylabel('Frequency')
plt.title('Correct Predictions')
plt.show();<jupyter_output><empty_output><jupyter_text>Compute metrics<jupyter_code>import sklearn.metrics as skm
import seaborn

# Classification report
report = skm.classification_report(y_test, y_pred)
print(report)

# Confusion matrix
cm = skm.confusion_matrix(y_test, y_pred)
seaborn.heatmap(cm, annot=True,annot_kws={"size": 16})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11508199829367713, 'fraction_numerical': 0.014788131576452744, 'mean_word_length': 5.180433509080258, 'pattern_counts': {'":': 1, '<': 30, '<?xml version=': 0, '>': 48, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '526571', 'license': 'Apache-2.0', 'n_tokens_mistral': 3844, 'n_tokens_neox': 3636, 'n_words': 787, 'path': 'BayesClassifierECG.ipynb', 'repo_name': 'mostly-sunny/ECGClassifier', 'size': '10557'}